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

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

例如，我想生成一個 4 x 4 隨機對稱正定矩陣 Ω∈R4×4 ，我們知道 Ω1,2=Ω2,1=Ω1,3=Ω3,1=0 . 我怎麼能在 R 中做到這一點？

---

我想到的是類似於 Cholesky 分解 LLT=Ω , 其中行 Li 和行 Lj 是正交的，如果 Ωij=0 . 可能通過拉格朗日乘數求解。但我不確定如何實現這一點。或者，如果這是可能的。

每一個 d×d 對稱正（半）定矩陣 Σ 可以被分解為

Σ=Λ′,Q′,Q,Λ

在哪裡 Q 是一個正交矩陣並且 Λ 是具有非負（正）項的對角矩陣 λ1,…,λd. ( Σ 總是一些的協方差矩陣 d - 變量分佈和 QQ′ 將是它的相關矩陣；這 λi 是邊際分佈的標準差。）

讓我們解釋一下這個公式。 這 (i,j) 入口 Σi,j 是列的點積 i 和 j 的 Q , 乘以 λiλj. 因此，零約束 Σ 是對列的點積的正交性約束 Q.

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

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

創建了一個實例 Q, 隨機生成對角線 Λ 任何你喜歡的方式（如在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

為方便起見，您可以通過對角線 Λ 作為 的第二個參數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