r - 如何在不进行矩阵求逆的情况下有效地计算 diag(X %*% solve(A) %*% t(X))？

Question

我需要以下对角线:

diag(X %*% solve(A) %*% t(X))

其中A 是满秩方阵，X 是矩形矩阵。 A 和 X 都是稀疏的。

我知道求逆矩阵很糟糕，除非你真的需要它。但是，我看不出如何重写公式，将 solve(A) 替换为具有两个参数的 solve，这样就可以求解线性系统而无需显式反转.这可能吗？

Answer 1

也许在我真正开始之前，我需要提一下，如果你这样做的话

diag(X %*% solve(A, t(X)))

避免矩阵求逆。 solve(A, B) 执行 LU 分解并使用生成的三角矩阵因子求解线性系统 A x = B。如果您未指定 B，它默认为对角矩阵，因此您将显式计算 A 的逆矩阵。

您应该仔细阅读 ?solve ，并可能多次阅读提示。它说它基于 LAPACK 例程 DGESV ，您可以在其中找到幕后的数值线性代数。

DGESV computes the solution to a real system of linear equations

   A * X = B,

where A is an N-by-N matrix and X and B are N-by-N RHS matrices.

The LU decomposition with partial pivoting and row interchanges is
used to factor A as

  A = P * L * U,

where P is a permutation matrix, L is unit lower triangular, and U is
upper triangular.  The factored form of A is then used to solve the
system of equations A * X = B.

solve(A, t(X)) 和 solve(A) %*% t(X) 之间的区别是效率问题。后者中的一般矩阵乘法 %*% 比 solve 本身要昂贵得多。

然而，即使您使用 solve(A, t(X))，您也不是在最快的轨道上，因为您还有另一个 %*%。

此外，由于您只需要对角线元素，因此您会在第一次获取完整矩阵时浪费大量时间。 我在下面的回答将使您走上最快的轨道。

我用 LaTeX 重写了所有内容，内容也大大扩展了，包括对 R 实现的引用。如果您发现它们有用，最后会提供有关 QR 分解、奇异值分解和 Pivoted Cholesky 分解的额外资源。

---

---

额外资源

如果您对枢轴 Cholesky 分解感兴趣，可以引用 Compute projection / hat matrix via QR factorization, SVD (and Cholesky factorization?) .在那里我还讨论了 QR 分解和奇异值分解。

以上链接设置在普通最小二乘回归上下文中。加权最小二乘可以引用Get hat matrix from QR decomposition for weighted least square regression .

QR 分解也采用紧凑形式。如果您想了解更多关于 QR 分解是如何完成和存储的，您可以引用 What is "qraux" returned by QR decomposition .

这些问题和答案都集中在数值矩阵计算上。下面给出一些统计应用:

• linear model with lm(): how to get prediction variance of sum of predicted values
• Sampling from multivariate normal distribution with rank-deficient covariance matrix via Pivoted Cholesky Factorization
• Boosting lm() by a factor of 3, using Cholesky method rather than QR method