2012-01-09 209 views
3

有點背景優先。我找到了一個實對稱矩陣的特徵值和特徵向量,其中行總和爲0.更具體地說,一旦我找到一個特徵向量,我使用$ argsort $來查找排列其中一個特徵值的排列,並將排列應用於原始矩陣。numpy中的精確性:比較數字時的問題

現在,我使用numpy包在python中實現了代碼。代碼本身是遞歸的,如果它在特徵向量中找到一組相等的值,則它提取對應於我們具有相等值的索引的對稱子矩陣,並將算法重新應用於此矩陣。雖然這一切都很好,而且大部分都是咕嚕咕嚕的工作,但是當一堆應該與特徵向量中的相等條目相對應的索引沒有被識別爲具有相等的值時,我感到驚訝。問題是這些值是通過一些算法計算出來的(可能是Lanczos,但我對numpy並不完全熟悉)。這是一個示例輸出,其中我明確檢查兩個條目之間的差在本徵向量:

>>> T=spectral.seriation(A,index) 

    columns [ 0 1 2 3 4 5 6 7 8 9 10 11] 

    [ 3.30289130e-01 -2.75240941e-01 -2.75240941e-01 3.30289130e-01 
    -2.75240941e-01 3.30289130e-01 -2.75240941e-01 3.30289130e-01 
    3.30289130e-01 -2.75240941e-01 -1.69794463e-16 -2.75240941e-01] 

    [ 4 6 9 1 2 11 10 0 5 7 8 3] 

    difference -5.55111512313e-17 

例程系列化()是一個遞歸函數。浮點數組是正在考慮的特徵向量,下面的數組給出了列的排序順序。請注意,列[4,6,9,1,2,11]具有相同的值。但是,特徵向量和特徵值計算總是近似值,實際上,當我輸出第9列和第2列中的入口之間的差值時,它不是零。如果算法應該將[4,6,9,1,2,11]組合在一起,那麼它只將[4,6,9]分組,並將其餘的組合放在另一個組中,向作品中投擲扳手。

所以問題是這樣的:有沒有一種方法來執行numpy的任意精度計算?如果做不到這一點,這個問題會是一個「好」的解決方法嗎?

另外,我應該提到,它可以在數學上證明這些條目必須相同。這是矩陣的一個屬性,但希望與這個問題沒有密切關係。

+0

正如@amit已經說過的,不要檢查與浮點數的相等性。檢查它們是否在容錯範圍內。如果你想要一個'numpy'函數,使用'numpy.allclose(a,b)'而不是'a == b'。爲了更直接地回答你的問題,雖然,numpy不支持任意精度。 (你可以使用'decimal'對象數組來僞造它,但是1)會破壞numpy數組的全部目的,2)'numpy.linalg'使用lapack,它不支持任意精度,所以你不會在這種情況下獲得任何東西。) – 2012-01-09 17:37:56

回答

4

雙打併不是真正的數字[甚至沒有理性]。在每個範圍內都有無數的有理數[好,每個範圍至少有兩個元素,準確地說],但只有有限的位數來表示它們。
因此,您應該期待「精確」計算出現一些舍入誤差。

更多inforamtion,你可能需要閱讀what every computer scientist should know about floating-point arithmetic

+0

確實,如果僅僅是找到特徵值和特徵向量,我可以接受精度的損失。畢竟,這些不一定是合理的,甚至不是終止小數。然而,對我來說重要的是,程序應該忽略由這些舍入錯誤引起的小差異。否則,就像我們在上面的例子中看到的那樣,列的分區會出錯。 – user1137683 2012-01-09 17:37:09

2

當執行同等規模的兩個浮點數的減法,精度不應該是一個問題,即,如果[2]和[9]真的是同樣的差異將是零。

我懷疑實際情況是,默認情況下,輸出將數字顯示爲8位小數,但除此之外,數字會有所不同,通常雙精度大約有16個小數位數(要找出運行numpy.finfo(numpy.float).eps以獲得該機器epsilon給出了以上最小可能的數字1)

嘗試使用輸出格式"%.16f\n%.16f" % myarray[[2, 9]]檢查數字。

如果他們確實不同,但您對7d.p的相似性感到滿意,那麼您可以使用像numpy.around(differences, 7)這樣的東西截短結果。另外,如果你想預處理數據,那麼你可以使用類似下面的東西(儘管可能有更好的方法來做到這一點)。

sigcnd, expn = numpy.frexp(myarray) 
sigcnd = numpy.around(sigcnd, 7) 
truncated_myarray = numpy.ldexp(sigcnd, expn)
1

如果你想直逼元素的給定容差,你可以做類似的指標:

def almost_matches(x, array, rtol=1e-05, atol=1e-08): 
    answer = [] 
    for y in xrange(len(array)): 
     if abs(x-array[y]) <= (atol + rtol * abs(array[y])): 
      answer.append(y) 
    return answer 

(使用相同近似的比較,因爲numpy.allclose()使用)

>>> a = [3.30289130e-01, -2.75240941e-01, -2.75240941e-01, 3.30289130e-01, -2.75240941e-01, 3.30289130e-01, -2.75240941e-01, 3.30289130e-01, 3.30289130e-01, -2.75240941e-01, -1.69794463e-16, -2.75240941e-01] 
>>> almost_matches(min(a), a) 
[1, 2, 4, 6, 9, 11] 
2

檢查numpy.allclosenumpy.isclose用於在公差範圍內測試相等性的函數。