當特徵值很小時，獲得準確的特徵向量
我有一個對稱矩陣
A。我無法準確計算所有特徵向量，我相信這是由於最後幾個特徵值A非常非常小。是否有一種方法/庫可以讓我獲得非常準確的特徵值/向量，即使它非常慢？（在 R 或 Python 中）
（一個建議是使用“投影儀”迭代地消除最大的特徵值/向量，但我不知道如何去做/是否可行）
編輯：
一些背景：我正在考慮的矩陣基本上是通過對多參數模型的輸出進行 Hessian 矩陣生成的 Fisher 信息矩陣，每個模型都有顯著的內部參數“冗餘”，您可以從特徵分解中發現。（在這些情況下獲得非常高的條件數並不罕見，如下所述：https://www.lassp.cornell.edu/sethna/Sloppy/WhatAreSloppyModels.html。）
這是一些 R 代碼，其中包含我感興趣的矩陣示例。
#A matrix we're interested in: m <- matrix(c( 0.002307414, -1.318435973, 8.63e-11, -3.33e-11, 0.486238053, 0.41182779, -0.169525454, 0.624275558, -1.318435973, 753.3425186, -4.93e-08, 1.90e-08, -277.8320732, -235.3147146, 96.86532753, -356.7054684, 8.63e-11, -4.93e-08, 3.23e-18, -1.24e-18, 1.82e-08, 1.54e-08, -6.34e-09, 2.34e-08, -3.33e-11, 1.90e-08, -1.24e-18, 4.79e-19, -7.01e-09, -5.94e-09, 2.44e-09, -9.00e-09, 0.486238053, -277.8320732, 1.82e-08, -7.01e-09, 102.4642297, 86.78386443, -35.72384951, 131.5526701, 0.41182779, -235.3147146, 1.54e-08, -5.94e-09, 86.78386443, 73.50310589, -30.25693671, 111.4208258, -0.169525454, 96.86532753, -6.34e-09, 2.44e-09, -35.72384951, -30.25693671, 12.45501408, -45.8654479, 0.624275558, -356.7054684, 2.34e-08, -9.00e-09, 131.5526701, 111.4208258, -45.8654479, 168.8989909 ), nrow=8) print('matrix m:') print(m) #Get eigenvalues/vectors eig <- eigen(m) print('Eigenvalues/vectors') print(eig) #Check - does Matrix*eigenvector == eigenvalue*eigenvector? (expected result is c(1,1,1,1,1,1,1,1) # print('Eigenvector 1') print((m %*% eig$vectors[,1]) / (eig$values[1] * eig$vectors[,1])) #Fine - yields c(1,1,1,1,1,1,1,1) print('Eigenvector 2') print((m %*% eig$vectors[,2]) / (eig$values[2] * eig$vectors[,2])) #Mostly ok - all approximately 1 print('Eigenvector 3 (bad!)') print((m %*% eig$vectors[,3]) / (eig$values[3] * eig$vectors[,3])) #---Bad (elements vary a lot)! print('Eigenvector 4 (really bad!)') print((m %*% eig$vectors[,4]) / (eig$values[4] * eig$vectors[,4])) #---Really bad! print('Eigenvector 5') print((m %*% eig$vectors[,5]) / (eig$values[5] * eig$vectors[,5])) #Mostly ok print('Eigenvector 6') print((m %*% eig$vectors[,6]) / (eig$values[6] * eig$vectors[,6])) #Mostly ok print('Eigenvector 7') print((m %*% eig$vectors[,7]) / (eig$values[7] * eig$vectors[,7])) #Mostly ok print('Eigenvector 8') print((m %*% eig$vectors[,8]) / (eig$values[8] * eig$vectors[,8])) #Mostly ok
問題是由於大特徵向量的“洩漏”。 我將介紹一個簡短的分析，然後提供一個帶有代碼的解決方案，然後對解決方案的性質和局限性進行一些評論。
分析
讓實對稱矩陣的特徵值 M 是 λ1,λ2,…,λd 按絕對值遞減排序，使得 |λ1|≥|λ2|≥⋯≥|λd|. 如對該問題的評論中所述，目標是找到向量 x 以便 Mx 比較小。因為映射 x→Mx 是線性的，我們可以集中研究這樣的 x 關於單位長度的。
譜定理保證正交基的存在 e1,e2,…,ed 對應的特徵向量，其中 x 可以表示為線性組合
x=ξ1e1+ξ2e2+⋯+ξded.
正交性向我們保證
1=|x|2=ξ21+ξ22+⋯+ξ2d.
申請 M 產量
Mx=ξ1Me1+ξ2Me2+⋯ξdMed=ξ1λ1e1+ξ2λ2e2+⋯+ξdλded.
假設在特徵向量的計算中存在一些輕微的錯誤 ei, 所以我們認為我們正在使用具有特徵值的特徵向量 λi, 我們確實在處理擾動 e′i 在哪裡
e′i−ei=αi1e1+αi2e2+⋯+αided.
這些不再是特徵向量了，因為我們可以通過應用來檢查 M 給他們：
Me′i=Mei+αi1Me1+⋯+αidMed=λiei+αi1λ1e1+⋯+αidλded.
特別注意，當 λi 相比起來很小 λ1, 複數的出現 αi1λ1 的 e1 可以極大地改變結果，即使 αi1 很小，因為乘以 λ1. 這是問題的癥結所在。
申請時問題蔓延 M 到線性組合
x′=ξ1e′1+ξ2e′2+⋯+ξde′d,
我們認為是 x （因為我們認為 e′i 是 ei ）。申請 M 分別對每個 e′i 在右邊“洩漏”潛在的大倍數 e1 進入結果。
解決方案
幸運的是，有一個簡單的解決方案：從結果中刪除意外的特徵向量。 當（說）第一個 k 的係數 x 為零， ξ1=ξ2=⋯=ξk, 那麼*不應該有任何倍數 e1 通過 ek 在 Mx. 我們可以通過投影結果來刪除它們 Mx 到由 ek+1,…,ed.
這並不能讓我們得到正確的值 Mx, 因為較小的特徵向量之間存在“洩漏”。例如，申請 M 到 x=ed 產生所有的線性組合 ei. 投射出第一個 k 特徵向量給我們留下了線性組合
λk+1αd,k+1ek+1+λk+2αd,k+2ek+2+⋯+λd(1+αd,d)ed.
但是，如果 (a) |λk+1|≈|λk+2|≈⋯≈|λd| 具有可比較的數量級，並且 (b) αd,∗ 係數都比較小，結果還是比較準確的。
例子
讓我們看看它是如何與矩陣一起工作的 M 的問題。這裡 k=8 和 λ1≈1111 和 λ2≈5×10−8 遠大於任何其他特徵向量。在預期（假設）非負值的設置中出現微小的負特徵值進一步證明了最後六個（實際上是七個）特徵值可能只是零的噪聲版本。
下面的
R代碼創建隨機向量 x′ 在最後生成的空間中 8−2=6 特徵向量。它適用 M 有兩種方式：直接作為 y=Mx′ ; 並使用“投影法”，其中 y 針對前兩個特徵向量進行回歸。它返回這兩個結果的規範。*我們希望規範真的很小，*因為 M 不應該改變規範超過 max|λ3|,…,|λ8|≈6.25×10−8. 為了了解它的真正作用，我繪製了兩組結果（左側和中間）的直方圖。為了比較這兩種方法，我還在右側繪製了它們的散點圖。
**直接法（左）很糟糕：**相當多的 e1 已洩漏到特徵向量中，導致規範膨脹了大約十（！）個數量級。在將前兩個特徵向量投影出來之後，範數總是具有預期的大小。散點圖顯示結果之間沒有關係。（在某些情況下，對於其他矩陣 M, 誤差係數的大小引起了一些有趣的關係 αij. )
評論
該方案保證了計算的結果 Mx 將與最大特徵空間正交，並以相對於最大特徵值大小的精度反映多少 M 影響長度。這可能就是預期應用程序所需的全部內容。
也沒有做得更好的前景。 矩陣中的許多條目如 M 已經使用浮點算術計算，因此不能認為比這些操作的結果更精確。矩陣的大條件數，如 M 保證即使使用無限精確的算術，（a）更改任何條目的效果 M 在其最低有效位或（b）改變任何係數 x 在其最低有效位將極大地擾亂結果。
（
svd奇異值分解）函數在R計算特徵值和向量方面做得更好。 這不會改善直接計算 M 可想而知，雖然。如果您想比較這些方法，請將行替換eig <- eigen(m)為eig <- with(svd(m), list(values = d, vectors = u)).代碼
這
m在問題中的定義和使用 計算其特徵向量之後得到了提升eigen。lambda <- eig$values i <- lambda >= max(lambda) * 1e-12 # Indexes of large eigenvalues. # # Generate unit movements in the small directions. # nu <- function(x) { # Standardize to unit length if possible a <- sum(x^2) if(a==0) x else x/sqrt(a) } set.seed(17) sim <- replicate(500, { # Naive (uncorrected) algorithm. dx <- eig$vectors %*% rnorm(length(lambda)) * ifelse(i, 0, 1) dy <- m %*% nu(dx) # Corrected algorithm with projection. dx <- residuals(lm(dx ~ eig$vectors[,i])) dz <- m %*% nu(dx) c(sqrt(sum(dy^2)), sqrt(sum(dz^2))) }) # # Compare. # par(mfrow=c(1,3)) hist(sim[1,], main="Magnitudes of resultants: direct method", cex.main=1) hist(sim[2,], main="Magnitudes of resultants: projection method", cex.main=1) plot(t(sim), xlab="Direct method", ylab="Projection method", main="Comparison", cex.main=1) par(mfrow=c(1,1))
