19

I have a positive definite matrix $A$ ($n \times n$ dimension) for which I have the Cholesky decomposition $A=LL^{'}$. I want to use this to compute

a) The Cholesky decomposition of $A+c^2\times I $ where $c$ is a constant and $I$ is the identity matrix

b) The Cholesky decomposition of $A+BB^{'}$ where $B$ is a $n \times n$ sparse matrix with each row having at most $k$ elements for some fixed $k << n$.

Is there any analytical/ computational method/ R-package that can use the already available Cholesky decomposition of $A$ and perform (a) and (b) in a computationally scalable way i.e. ($O(n)$ complexity). Note that (a) is a special case of (b) with $B=cI$. Any references will be appreciated.

Thanks

Rodrigo de Azevedo
  • 1
Abhirup Datta
  • 191
  • could you please tell us, how you intend to use the cholesky factors – andre Feb 21 '16 at 11:31
  • It is easy to see that in general the entries of the $L$ for $A$ are totally different from the entries of the $L$ of $A+BB'$, so a complexity $O(n)$ is impossible, since we have to compute $O(n^2)$ entries for the new $L$. – Giovanni Resta Feb 21 '16 at 23:33
  • I'm very interested in this question as well. Using a solver for $A$ to build a solver for $A + \sigma I$ would allow one to efficiently compute rational functions of the matrix, $f(A) = p(A)/q(A) = p(A) \left(c_1(A + \sigma_1 I)^{-1} + \dots + c_n(A + \sigma_n I)^{-1}\right)$, where $c_i, \sigma_i$ come from the partial fractions expansion. In turn, this can be used to compute any function of a matrix that is well-approximated by rational functions (e.g., exp(A), sqrt(A), etc...). – Nick Alger Feb 22 '16 at 05:07
  • 2
    Related earlier question: http://math.stackexchange.com/q/940342/11268 I agree with Nick Alger that such an algorithm would be tremendously useful. The fact that it does not appear in research literature (esp. related to matrix functions) suggests to me that no such algorithm is currently known. More related questions: http://scicomp.stackexchange.com/q/10278/713, http://scicomp.stackexchange.com/q/8323/713, http://scicomp.stackexchange.com/q/21717/713, all negative. – Kirill Feb 22 '16 at 23:56

1 Answers1

10

I don't see what you can expect.

The direct way is to solve the system in the unknown triangular matrix $X$: $XX^T=A+BB^T,x_{i,i}>0$.

If you want to use your decomposition $A=LL^T$, then $A+BB^T=(L+Y)(L+Y)^T$; the system to solve is now $YY^T+YL^T+LY^T-BB^T=0,y_{i,i}>0$; clearly, this system is more complicated than the above one!

Note that, if $B$ is a small matrix, then we obtain a good approximation by solving the linear equation $YL^T+LY^T-BB^T=0$. This system is stable when the $l_{i,i}$ are not too small.