當特徵值很小時，獲得準確的特徵向量
我有一個對稱矩陣
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
問題是由於大特徵向量的“洩漏”。 我將介紹一個簡短的分析，然後提供一個帶有代碼的解決方案，然後對解決方案的性質和局限性進行一些評論。
分析
讓實對稱矩陣的特徵值 $ \mathbb M $ 是 $ \lambda_1,\lambda_2,\ldots,\lambda_d $ 按絕對值遞減排序，使得 $ |\lambda_1|\ge |\lambda_2|\ge \cdots \ge |\lambda_d|. $ 如對該問題的評論中所述，目標是找到向量 $ x $ 以便 $ \mathbb{M}x $ 比較小。因為映射 $ x\to\mathbb{M}x $ 是線性的，我們可以集中研究這樣的 $ x $ 關於單位長度的。
譜定理保證正交基的存在 $ {e_1,e_2,\ldots, e_d} $ 對應的特徵向量，其中 $ x $ 可以表示為線性組合
$$ x = \xi_1 e_1 + \xi_2 e_2 + \cdots + \xi_d e_d. $$
正交性向我們保證
$$ 1 = |x|^2 = \xi_1^2 + \xi_2^2 + \cdots + \xi_d^2. $$
申請 $ \mathbb M $ 產量
$$ \mathbb{M}x = \xi_1 \mathbb{M} e_1 + \xi_2 \mathbb{M}e_2 + \cdots \xi_d \mathbb{M}e_d = \xi_1\lambda_1 e_1 + \xi_2\lambda_2 e_2 + \cdots + \xi_d\lambda_d e_d. $$
假設在特徵向量的計算中存在一些輕微的錯誤 $ e_i, $ 所以我們認為我們正在使用具有特徵值的特徵向量 $ \lambda_i, $ 我們確實在處理擾動 $ e_i^\prime $ 在哪裡
$$ e_i^\prime - e_i = \alpha_{i1}e_1 + \alpha_{i2}e_2 + \cdots + \alpha_{id}e_d. $$
這些不再是特徵向量了，因為我們可以通過應用來檢查 $ \mathbb M $ 給他們：
$$ \mathbb{M}e_i^\prime = \mathbb{M}e_i + \alpha_{i1}\mathbb{M}e_1 + \cdots + \alpha_{id}\mathbb{M}e_d = \lambda_i e_i + \alpha_{i1}\lambda_1 e_1 + \cdots + \alpha_{id}\lambda_d e_d. $$
特別注意，當 $ \lambda_i $ 相比起來很小 $ \lambda_1, $ 複數的出現 $ \alpha_{i1}\lambda_1 $ 的 $ e_1 $ 可以極大地改變結果，即使 $ \alpha_{i1} $ 很小，因為乘以 $ \lambda_1. $ 這是問題的癥結所在。
申請時問題蔓延 $ \mathbb{M} $ 到線性組合
$$ x^\prime = \xi_1 e_1^\prime + \xi_2 e_2^\prime + \cdots + \xi_d e_d^\prime, $$
我們認為是 $ x $ （因為我們認為 $ e_i^\prime $ 是 $ e_i $ ）。申請 $ \mathbb{M} $ 分別對每個 $ e_i^\prime $ 在右邊“洩漏”潛在的大倍數 $ e_1 $ 進入結果。
解決方案
幸運的是，有一個簡單的解決方案：從結果中刪除意外的特徵向量。 當（說）第一個 $ k $ 的係數 $ x $ 為零， $ \xi_1=\xi_2=\cdots=\xi_k, $ 那麼*不應該有任何倍數 $ e_1 $ 通過 $ e_k $ 在 $ \mathbb{M} x. $ 我們可以通過投影結果來刪除它們 $ \mathbb{M}x $ 到由 $ e_{k+1},\ldots, e_d. $
這並不能讓我們得到正確的值 $ \mathbb{M}x, $ 因為較小的特徵向量之間存在“洩漏”。例如，申請 $ \mathbb{M} $ 到 $ x=e_d $ 產生所有的線性組合 $ e_i. $ 投射出第一個 $ k $ 特徵向量給我們留下了線性組合
$$ \lambda_{k+1}\alpha_{d,k+1} e_{k+1} + \lambda_{k+2}\alpha_{d,k+2} e_{k+2} + \cdots + \lambda_{d}(1+\alpha_{d,d}) e_{d}. $$
但是，如果 (a) $ |\lambda_{k+1}| \approx |\lambda_{k+2}| \approx \cdots \approx |\lambda_d| $ 具有可比較的數量級，並且 (b) $ \alpha_{d,*} $ 係數都比較小，結果還是比較準確的。
例子
讓我們看看它是如何與矩陣一起工作的 $ \mathbb M $ 的問題。這裡 $ k=8 $ 和 $ \lambda_1\approx 1111 $ 和 $ \lambda_2 \approx 5\times 10^{-8} $ 遠大於任何其他特徵向量。在預期（假設）非負值的設置中出現微小的負特徵值進一步證明了最後六個（實際上是七個）特徵值可能只是零的噪聲版本。
下面的
R代碼創建隨機向量 $ x^\prime $ 在最後生成的空間中 $ 8-2=6 $ 特徵向量。它適用 $ \mathbb M $ 有兩種方式：直接作為 $ y=\mathbb{M}x^\prime $ ; 並使用“投影法”，其中 $ y $ 針對前兩個特徵向量進行回歸。它返回這兩個結果的規範。*我們希望規範真的很小，*因為 $ \mathbb M $ 不應該改變規範超過 $ \max{|\lambda_3|,\ldots, |\lambda_8|}\approx 6.25\times 10^{-8}. $ 為了了解它的真正作用，我繪製了兩組結果（左側和中間）的直方圖。為了比較這兩種方法，我還在右側繪製了它們的散點圖。
**直接法（左）很糟糕：**相當多的 $ e_1 $ 已洩漏到特徵向量中，導致規範膨脹了大約十（！）個數量級。在將前兩個特徵向量投影出來之後，範數總是具有預期的大小。散點圖顯示結果之間沒有關係。（在某些情況下，對於其他矩陣 $ \mathbb M, $ 誤差係數的大小引起了一些有趣的關係 $ \alpha_{ij}. $ )
評論
該方案保證了計算的結果 $ \mathbb{M}x $ 將與最大特徵空間正交，並以相對於最大特徵值大小的精度反映多少 $ \mathbb M $ 影響長度。這可能就是預期應用程序所需的全部內容。
也沒有做得更好的前景。 矩陣中的許多條目如 $ \mathbb M $ 已經使用浮點算術計算，因此不能認為比這些操作的結果更精確。矩陣的大條件數，如 $ \mathbb M $ 保證即使使用無限精確的算術，（a）更改任何條目的效果 $ \mathbb M $ 在其最低有效位或（b）改變任何係數 $ x $ 在其最低有效位將極大地擾亂結果。
（
svd奇異值分解）函數在R計算特徵值和向量方面做得更好。 這不會改善直接計算 $ \mathbb 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))
