2
$\begingroup$

Let $D$ be an arbitrary diagonal matrix and let $A$ be an orthogonal matrix ($A'A = AA' = I$). How to compute the following matrix inverse efficiently?

$$(D + ADA^T)^{-1}$$

Hints or references are also welcome. Thanks a lot!

$\endgroup$
  • 3
    $\begingroup$ note that $(D+ADA^T)^{-1}=A(DA+AD)^{-1}$; I would be surprised if $DA+AD$ has any structure that would simplify finding an inverse. $\endgroup$ – Carlo Beenakker Aug 6 '18 at 6:15
  • 2
    $\begingroup$ We get a lot of questions about "update the inverse of a matrix after adding a diagonal correction" at scicomp.se (see e.g. scicomp.stackexchange.com/a/10282/4405), and the answer is almost invariably that there's nothing better than recomputing the inverse from scratch. There is a little more structure here, but my first guess would be no (unless $D=I$). $\endgroup$ – Federico Poloni Aug 6 '18 at 6:36
  • $\begingroup$ Are $A'$ and $A^T$ the same thing? $\endgroup$ – Federico Poloni Aug 8 '18 at 6:56
6
$\begingroup$

To begin with, the matrix in question can well be degenerate, consider for example

\begin{equation*} D=\left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array}\right), A=\left( \begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array}\right). \end{equation*}

I suspect that because of this you should not count for much for a general $D$.

Things are better when $D>0$, because in this case $$D+ADA^T=D^{1/2}(I+BB^T)D^{1/2},\,B=D^{-1/2}AD^{1/2}.$$ Computing $(I+C)^{-1}$ where $C$ is a positive definite matrix with a known inverse $C^{-1}$ is a standard problem (computing it numerically, that is). As for good references, I simply do not remember (I am through with numerical analysis since a decade ago), but the wiki page https://en.wikipedia.org/wiki/Numerical_linear_algebra may be a good place to start.

$\endgroup$
2
$\begingroup$

I would like to add to Alex's idea. It is good to first make note of the fact that the given matrix which we are trying to invert is a symmetric matrix. Additionally, like Alex suggested having D>0 helps a lot, as, it actually makes the whole matrix to be positive definite (not only C (from Alex's answer)) (obviously! because of the additional Identity matrix) (actually, the fact that it is positive definite is more clearly visible if we keep it in the given original form). So, if D>0, the given matrix as whole is Symmetric Positive Definite (SPD) and, having it to be SPD makes things very comfortable (Both, to choose the right iterative method and also to prove the convergence of the method). At this point you have to ask if the matrix is Sparse or Dense (as you are asking about efficiency, I presume you are dealing with big matrix sizes). So, if it is Sparse then, Conjugate Gradient is the best to choose for SPD matrices (https://en.wikipedia.org/wiki/Conjugate_gradient_method) else, if the matrix is Dense, Cholesky decomposition is the way to go for SPD matrices and also for present structure in general as Federico's comment suggests (more efficient than LU decomposition) (https://en.wikipedia.org/wiki/Cholesky_decomposition). Additionally if the matrix is sparse, to make things more efficient in terms of storage while implementation, it is better to use sparse storage data structures (see https://en.wikipedia.org/wiki/Sparse_matrix).

(PS: the suggestions are from the well established methods, you can additionally look for some preconditioning which helps a lot for little effort! (see https://en.wikipedia.org/wiki/Preconditioner), Multigrid could a nice choice I believe for an SPD matrix for preconditioning.)

$\endgroup$
0
$\begingroup$

An answer to your question heavily depends on whether your matrix is sparse or dense. If it is dense (or even sparse, but not larger than a few thousands rows/columns), than any direct method, like based on LU-factorisation, will do. If it is sparse, the problem becomes a lot more complex. I suggest you visit the following blog page and ask any questions there: http://comecau.blogspot.com/2018_09_05_archive.html

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.