I need the following diagonal:
diag(X %*% solve(A) %*% t(X))
where A
is a full rank square matrix, X
is a rectangular matrix. Both A
and X
are sparse.
I know finding the inverse of a matrix is bad unless you really need it. However, I can't see how to rewrite the formula such that solve(A)
is replaced by solve
with two arguments, such that a linear system is solved without explicitly inverting. Is that possible?
Perhaps before I really start, I need to mention that if you do
matrix inverse is avoided.
solve(A, B)
performs LU factorization and uses the resulting triangular matrix factors to solve linear systemA x = B
. If you leaveB
unspecified, it is default to a diagonal matrix hence you will be explicitly computing matrix inverse ofA
.You should read
?solve
carefully, and possibly many times for hints. It says it is based onLAPACK
routineDGESV
, where you can find the numerical linear algebra behind the scene.The difference between
solve(A, t(X))
andsolve(A) %*% t(X)
is a matter of efficiency. The general matrix multiplication%*%
in the latter is much more expensive thansolve
itself.Yet, even if you use
solve(A, t(X))
, you are not on the fastest track, as you have another%*%
.Also, as you only want diagonal elements, you waste considerable time in first getting the full matrix. My answer below will get you on the fastest track.
I have rewritten everything in LaTeX, and the content is greatly expanded, too, including reference to R implementation. Extra resources are given on QR factorization, Singular Value decomposition and Pivoted Cholesky factorization in the end, in case you find them useful.
Extra resources
In case you are interested in pivoted Cholesky factorization, you can refer to Compute projection / hat matrix via QR factorization, SVD (and Cholesky factorization?). There I also discuss QR factorization and singular value decomposition.
The above link is set in ordinary least square regression context. For weighted least square, you can refer to Get hat matrix from QR decomposition for weighted least square regression.
QR factorization also takes compact form. Should you want to know more on how QR factorization is done and stored, you may refer to What is "qraux" returned by QR decomposition.
These questions and answers all focus on numerical matrix computations. The following gives some statistical application:
lm()
: how to get prediction variance of sum of predicted valueslm()
by a factor of 3, using Cholesky method rather than QR method