使用 R 生成具有零約束的隨機正定矩陣

如何使用 R 生成具有零約束的隨機對稱正定矩陣？

例如，我想生成一個 4 x 4 隨機對稱正定矩陣 $ \Omega\in\mathbb{R}^{4\times4} $ ，我們知道 $ \Omega_{1,2}=\Omega_{2,1}=\Omega_{1,3}=\Omega_{3,1} = 0 $ . 我怎麼能在 R 中做到這一點？

---

我想到的是類似於 Cholesky 分解 $ LL^T=\Omega $ , 其中行 $ L_i $ 和行 $ L_j $ 是正交的，如果 $ \Omega_{ij}=0 $ . 可能通過拉格朗日乘數求解。但我不確定如何實現這一點。或者，如果這是可能的。

每一個 $ d\times d $ 對稱正（半）定矩陣 $ \Sigma $ 可以被分解為

$$ \Sigma = \Lambda^\prime, Q^\prime ,Q,\Lambda $$

在哪裡 $ Q $ 是一個正交矩陣並且 $ \Lambda $ 是具有非負（正）項的對角矩陣 $ \lambda_1, \ldots, \lambda_d. $ ( $ \Sigma $ 總是一些的協方差矩陣 $ d $ - 變量分佈和 $ QQ^\prime $ 將是它的相關矩陣；這 $ \lambda_i $ 是邊際分佈的標準差。）

讓我們解釋一下這個公式。 這 $ (i,j) $ 入口 $ \Sigma_{i,j} $ 是列的點積 $ i $ 和 $ j $ 的 $ Q $ , 乘以 $ \lambda_i\lambda_j. $ 因此，零約束 $ \Sigma $ 是對列的點積的正交性約束 $ Q. $

（請注意，正定矩陣的所有對角元素都必須是非零的，所以我假設零約束都在對角線之外。我還擴展了對 $ (i,j) $ 進入一個約束 $ (j,i) $ 條目，以確保結果的對稱性。）

施加此類約束的一種（完全通用的）方法是生成 $ Q $ 依次。 使用您喜歡的任何方法創建一個 $ d\times d $ 初始值矩陣。在步驟 $ i=1,2,\ldots, d, $ 更改列 $ i $ 在所有列上對其進行回歸 $ 1, 2, \ldots, i-1 $ 的 $ Q $ 這需要與它正交並保留殘差。將這些結果歸一化，使它們的點積（平方和）一致。那是列 $ i $ 的 $ Q. $

創建了一個實例 $ Q, $ 隨機生成對角線 $ \Lambda $ 任何你喜歡的方式（如在https://stats.stackexchange.com/a/215647/919密切相關的答案中所討論的）。

以下R函數默認rQ使用iid標準正態變量作為初始值。我已經用尺寸對其進行了廣泛的測試 $ d=1 $ 通過 $ 200, $ 系統地檢查預期的約束是否成立。我還用泊松測試了它 $ (0.1) $ 變量，因為它們很可能是零，所以會產生非常有問題的初始解決方案。

的主要輸入rQ是一個邏輯矩陣，指示要在何處應用零約束。這是問題中指定的約束的示例。

set.seed(17)
Q <- matrix(c(FALSE, TRUE, TRUE, FALSE,
              TRUE, FALSE, FALSE, FALSE,
              TRUE, FALSE, FALSE, FALSE,
              FALSE, FALSE, FALSE, FALSE), 4)
Lambda <- rexp(4)
zapsmall(rQ(Q, Lambda))

        [,1]      [,2]      [,3]      [,4]
[1,] 2.646156  0.000000  0.000000  2.249189
[2,] 0.000000  0.079933  0.014089 -0.360013
[3,] 0.000000  0.014089  0.006021 -0.055590
[4,] 2.249189 -0.360013 -0.055590  4.167296

為方便起見，您可以通過對角線 $ \Lambda $ 作為 的第二個參數rQ。它的第三個參數f必須是一個隨機數生成器（或任何其他f(n)返回長度為 的數值向量的函數n）。

rQ <- function(Q, Lambda, f=rnorm) {
 normalize <- function(x) {
   v <- zapsmall(c(1, sqrt(sum(x * x))))[2]
   if (v == 0) v <- 1
   x / v
 }
 Q <- Q | t(Q)                    # Force symmetry by applying all constraints
 d <- nrow(Q) 
 if (missing(Lambda)) Lambda <- rep(1, d)
 R <- matrix(f(d^2), d, d)        # An array of column vectors
 for (i in seq_len(d)) {
   j <- which(Q[seq_len(i-1), i]) # Indices of the preceding orthogonal vectors
   R[, i] <- normalize(residuals(.lm.fit(R[, j, drop=FALSE], R[, i])))
 }
 R <- R %*% diag(Lambda)
 crossprod(R)
}

引用自：https://stats.stackexchange.com/questions/547035