隨機數據的 SVD 結果中的奇怪相關性；他們有數學解釋還是 LAPACK 錯誤？
我在隨機數據的 SVD 結果中觀察到一個非常奇怪的行為，我可以在 Matlab 和 R 中重現這種行為。這看起來像是 LAPACK 庫中的一些數值問題；是嗎？
我畫樣本來自具有零均值和恆等協方差的維高斯：. 我把它們組裝成一個數據矩陣. （我可以選擇居中與否，它不影響以下內容。）然後我執行奇異值分解（SVD）得到. 讓我們採取一些兩個特定的元素，例如和，並詢問它們在不同的抽獎中的相關性是什麼. 我希望如果數字of draws 相當大，那麼所有這些相關性應該在零附近（即總體相關性應該為零，樣本相關性會很小）。
然而，我觀察到一些奇怪的強相關性（大約） 之間,,， 和，並且僅在這些元素之間。正如預期的那樣，所有其他元素對的相關性大約為零。以下是相關矩陣的方式“上”元素看起來像（首先第一列的元素，然後是第一列第二列的元素）：
請注意每個像限左上角的異常高值。
正是這個@whuber 的評論引起了我的注意。@whuber 認為 PC1 和 PC2 不是獨立的，並提出了這種強相關性作為證據。然而，我的印像是他無意中發現了 LAPACK 庫中的一個數字錯誤。這裡發生了什麼？
這是@whuber 的 R 代碼：
stat <- function(x) {u <- svd(x)$u; c(u[1,1], u[2, 2])}; Sigma <- matrix(c(1,0,0,1), 2); sim <- t(replicate(1e3, stat(MASS::mvrnorm(10, c(0,0), Sigma)))); cor.test(sim[,1], sim[,2]);這是我的 Matlab 代碼：
clear all rng(7) n = 1000; %// Number of variables k = 2; %// Number of observations Nrep = 1000; %// Number of iterations (draws) for rep = 1:Nrep X = randn(n,k); %// X = bsxfun(@minus, X, mean(X)); [U,S,V] = svd(X,0); t(rep,:) = [U(1:10,1)' U(1:10,2)']; end figure imagesc(corr(t), [-.5 .5]) axis square hold on plot(xlim, [10.5 10.5], 'k') plot([10.5 10.5], ylim, 'k')
這不是錯誤。
正如我們在評論中（廣泛地）探索的那樣，發生了兩件事。第一個是列 $ U $ 被約束以滿足 SVD 要求：每個都必須具有單位長度並且與所有其他的正交。查看 $ U $ 作為從隨機矩陣創建的隨機變量 $ X $ 通過特定的 SVD 算法，我們因此註意到這些 $ k(k+1)/2 $ 功能獨立的約束在 $ U $ .
通過研究組件之間的相關性，可以或多或少地揭示這些依賴關係。 $ U $ ，但出現了第二個現象：SVD 解決方案不是唯一的。至少，每一列 $ U $ 可以獨立否定，至少給出 $ 2^k $ 不同的解決方案 $ k $ 列。強相關性（超過 $ 1/2 $ ) 可以通過適當地改變列的符號來誘導。（在此線程中我對 Amoeba 的回答的第一條評論中給出了一種方法：我強制所有 $ u_{ii},i=1,\ldots, k $ 具有相同的符號，使它們全部為負或以相等的概率全部為正。）另一方面，可以通過以相等的概率隨機、獨立地選擇符號來使所有相關性消失。（我在下面的“編輯”部分給出了一個例子。）
仔細閱讀時，我們可以部分辨別出這兩種現象 $ U $ . 某些特徵——例如在明確定義的圓形區域內幾乎均勻分佈的點的出現——表明缺乏獨立性。其他的，例如顯示明顯非零相關性的散點圖，顯然取決於算法中所做的選擇——但這種選擇之所以可能，是因為首先缺乏獨立性。
對像 SVD（或 Cholesky、LR、LU 等）這樣的分解算法的最終測試是它是否符合它的要求。在這種情況下，檢查當 SVD 返回矩陣的三元組時就足夠了 $ (U, D, V) $ ， 那 $ X $ 由產品恢復到預期的浮點錯誤 $ UDV^\prime $ ; 的列 $ U $ 和 $ V $ 是正交的；然後 $ D $ 是對角線，它的對角線元素是非負的，並且按降序排列。
svd我已經在算法中應用了這樣的測試，R並且從未發現它是錯誤的。儘管這不能保證它是完全正確的，但這種經驗——我相信很多人都分享過——表明任何錯誤都需要某種非凡的輸入才能顯現出來。以下是對問題中提出的具體觀點的更詳細分析。
使用
R’ssvd程序，首先您可以將其檢查為 $ k $ 增加，係數之間的相關性 $ U $ 變弱，但它們仍然非零。 如果您只是進行更大的模擬，您會發現它們很重要。（什麼時候 $ k=3 $ , 50000 次迭代應該足夠了。）與問題中的斷言相反，相關性並沒有“完全消失”。其次，研究這種現象的更好方法是回到係數*獨立性的基本問題。*儘管在大多數情況下相關性往往接近於零，但缺乏獨立性是顯而易見的。 通過研究係數的完整多元分佈，這一點最為明顯 $ U $ . 即使在非零相關性（還）無法檢測到的小型模擬中，分佈的性質也會出現。例如，檢查係數的散點圖矩陣。為了使這可行，我將每個模擬數據集的大小設置為 $ 4 $ 並保持 $ k=2 $ ，從而繪製 $ 1000 $ 的實現 $ 4\times 2 $ 矩陣 $ U $ , 創建一個 $ 1000\times 8 $ 矩陣。這是其完整的散點圖矩陣，其中的變量按它們的位置列出 $ U $ ：
向下掃描第一列揭示了一個有趣的缺乏獨立性 $ u_{11} $ 和另外一個 $ u_{ij} $ ：看看散點圖的上象限如何 $ u_{21} $ 例如，幾乎是空的；或檢查描述 $ (u_{11}, u_{22}) $ 關係和向下傾斜的雲 $ (u_{21}, u_{12}) $ 一對。 **仔細觀察會發現幾乎所有這些係數都明顯缺乏獨立性：**儘管它們中的大多數表現出接近於零的相關性，但它們中的極少數看起來是完全獨立的。
（注意：大多數圓形雲是由歸一化條件創建的超球面的投影，該條件迫使每列的所有分量的平方和為一。）
散點圖矩陣 $ k=3 $ 和 $ k=4 $ 表現出相似的模式：這些現象並不局限於 $ k=2 $ ，它們也不依賴於每個模擬數據集的大小：它們只是變得更難以生成和檢查。
這些模式的解釋轉到用於獲得的算法 $ U $ 在奇異值分解中，但我們知道這種非獨立模式必須存在於 $ U $ ：由於每個連續列（在幾何上）與前面的列正交，因此這些正交性條件在係數之間施加了函數依賴關係，從而轉化為相應隨機變量之間的統計依賴關係。
編輯
作為對評論的回應，可能值得一提的是，這些依賴現像在多大程度上反映了底層算法（計算 SVD），以及它們在多大程度上是過程本質所固有的。
**係數之間相關性的特定模式很大程度上取決於 SVD 算法做出的任意選擇，**因為解決方案不是唯一的： $ U $ 總是可以獨立地乘以 $ -1 $ 或者 $ 1 $ . 沒有內在的方式來選擇標誌。因此，當兩個 SVD 算法做出不同的（任意甚至隨機的）符號選擇時，它們可能會導致不同的散點圖模式 $ (u_{ij}, u_{i^\prime j^\prime}) $ 價值觀。如果您想看到這個，請將
stat下面代碼中的函數替換為stat <- function(x) { i <- sample.int(dim(x)[1]) # Make a random permutation of the rows of x u <- svd(x[i, ])$u # Perform SVD as.vector(u[order(i), ]) # Unpermute the rows of u }這首先隨機重新
x排序觀察，執行 SVD，然後應用逆序u來匹配原始觀察序列。因為效果是形成原始散點圖的反射和旋轉版本的混合，所以矩陣中的散點圖看起來更加均勻。所有樣本相關性都將非常接近於零（通過構造：基礎相關性正好為零）。儘管如此，缺乏獨立性仍然很明顯（在出現的均勻圓形中，特別是在 $ u_{i,j} $ 和 $ u_{i,j^\prime} $ ）。一些原始散點圖（如上圖所示）的某些象限中缺少數據是由於
RSVD 算法如何為列選擇符號。結論沒有任何變化。 因為第二列 $ U $ 與第一個正交，它（被認為是多元隨機變量）依賴於第一個（也被認為是多元隨機變量）。您不能讓一列的所有組件獨立於另一列的所有組件；您所能做的就是以模糊依賴關係的方式查看數據——但依賴關係將持續存在。
R這是處理案例的更新代碼 $ k\gt 2 $ 並繪製散點圖矩陣的一部分。k <- 2 # Number of variables p <- 4 # Number of observations n <- 1e3 # Number of iterations stat <- function(x) as.vector(svd(x)$u) Sigma <- diag(1, k, k); Mu <- rep(0, k) set.seed(17) sim <- t(replicate(n, stat(MASS::mvrnorm(p, Mu, Sigma)))) colnames(sim) <- as.vector(outer(1:p, 1:k, function(i,j) paste0(i,",",j))) pairs(sim[, 1:min(11, p*k)], pch=".")
