2

If $A^{-1}$ has been precomputed (or to be more precise: the Cholesky decomposition of A has been precomputed and cached), is there an efficient way to compute either $$C = (A+λI)^{-1}$$ or (more valuably)

the Cholesky Decomposition of: $$(A+λI)$$

A few other (possibly unhelpful/unecessary) notes:

  1. A is a square matrix that was constructed by doing this: $B^TB$, where B was some dense, non-square matrix of typically around 30 rows by 20 columns.
  2. Although A is small (eg 20 by 20), I have ~100 different A matrixes and ~$10^7$ different λ. So, precomputing 100 $A^{-1}$ (or its Cholesky decomposition) takes very little time; but having to compute the inverse of $10^7$ different $(A+λI)$ can take time, unless there's a clever way of taking advatange of the fact that $A^{-1}$ has been precomputed, and we're only adding a scaled version of the Identity matrix....
  3. I have already reviewed Inverse of the sum of matrices but no one seemed to have answered the question/comment by Royi, which I'm guessing was asked for the same reason I'm asking this question above.
  4. Ultimately, I will use the Cholesky decomposition of $(A+λI)$ to solve a linear equation, eg find $x$ in: $Cx = y$, where $C = (A+λI)$ and $y$ is known.
Jonathan Sylvester
  • 239
  • Let's see... The "normal" way to compute matrix inverse takes O(n^3). If you can compute each inverse matrix in optimal time O(n^2) then you get a improvement factor of 20, divide by constant factors. Not sure if it worth the effort. / Only single-threaded solutions are useful in our case. - Irrelevant for mathematics, only relevant in computer algorithm. Is this the wrong stackexchange? Should this question be on computer science instead? – user202729 Dec 06 '17 at 14:51
  • In the linked question you can take $B=\lambda I$ in any of the answers. – Ross Millikan Dec 06 '17 at 15:01
  • @user202729 I took out the reference to "single threaded". If the constant factor is small (eg < 5), then it's still worthwile as the difference between an algorithm than runs in say 60 minutes vs 60/4 = 15 minutes is noticeable for an end-user. But, I get your point...it sounds like your hunch is that the constant factor for any solution (if it exists at all) would be >> 5 etc.? – Jonathan Sylvester Dec 06 '17 at 15:03
  • @RossMillikan I did, but in none of the answers is it an improvement in term of complexity. In fact, it might even take longer to compute the alternative representations that just doing the normal inverse of C, right? (Or am I missing something.) – Jonathan Sylvester Dec 06 '17 at 15:05
  • Let's wait for answers to see. – user202729 Dec 06 '17 at 15:07
  • I recall a section in Numerical Recipes that had a way to update an inverse for a small change in the matrix. I think it was in just one element, so you would have to do it $20$ times, but it would work if $\lambda$ is small enough. It is akin to the Taylor series in one of the answers in the linked question. I think that is your best approach, but again it relies on $\lambda$ being small. – Ross Millikan Dec 06 '17 at 15:25

1 Answers1

3

I am not sure how tied you are to the Cholesky factorization, but the following alternative approach may suit your needs. Observe that since $A = B^TB$, it is symmetric and hence orthogonally diagonalizable, i.e., $A = QDQ^T$, where $Q$ is an orthogonal matrix and $D$ is the diagonal matrix of the real eigenvalues of $A$. Now use the fact that $A + \lambda I = Q(D+\lambda I)Q^T$, which implies that $(A+\lambda I)^{-1} = Q(D + \lambda I)^{-1} Q^T$. Therefore, if you precompute $Q$ and $D$ for each $A$, then you can compute $(A+\lambda I)^{-1}$ efficiently for each $\lambda$ via a few matrix vector products. Since your matrices $A$ are quite small and there are only $\approx 100$ of them, it should be relatively cheap to compute and store these factors.

If you have access to the $B$ matrices, then it may be preferable to compute the singular value decomposition $B = U\Sigma V^T$, and then form the eigendecomposition of $A$ by $$ A = B^TB = V(\Sigma^T\Sigma)V^T. $$ Note that you only need the singular values and the right singular vectors.

K. Miller
  • 4,688
  • This is exactly what I was looking for. I'm going to test it out, but based on my theoretical understanding (and reading), this ought to work. – Jonathan Sylvester Dec 07 '17 at 15:33