I am looking for a computationally cheap way to compute $x$ such that $$(L L^T + \mu^2 I)x = y$$ where $L \in \mathbb{R}^{n \times n}$ is a lower triangular definite positive matrix (with some very small eigenvalues), $y \in \mathbb{R}^n$ and $\mu \in \mathbb{R}$ are known. If it is necessary, I can assume that $$\mu \ll 1$$ But $\mu^2$ is larger than the smallest eigenvalue of $LL^T$.

Basically, I would like to make the most of my knowledge of the Cholesky decomposition $L L^T$. Eventually, I hope to be able to compute $x$ in $\mathcal{O}(n^2)$. Approximate approaches are also welcomed.

I have seen here that this does not seem to be doable in a more general situation, but I hope the smallness of $\mu$ may help...

Any idea, reference or warning? Thanks for your help.


Let $A=LL^T,\lambda=\mu^2,f:X\rightarrow X^{-1}$. That follows is an approximate approach in $O(n^2)$ that is valid only if $\lambda$ is small with respect to $\inf(spectrum(A))$. $Df_A(H)=-A^{-1}HA^{-1},D^2f_A(H,K)=A^{-1}KA^{-1}HA^{-1}+A^{-1}HA^{-1}KA^{-1}$. Thus, according to the Taylor formula, $(A+\lambda I)^{-1}\approx A^{-1}-\lambda A^{-2}+\lambda^2 A^{-3}$. Do not calculate $A^{-2},A^{-3}$ but solve $LL^Tx_0=y,LL^TLL^Tx_1=y,LL^TLL^TLL^Tx_2=y$, that is $LL^Tx_0=y,LL^Tx_1=x_0,LL^Tx_2=x_1$.

  • $\begingroup$ Unfortunately, the role of $\lambda = \mu^2$ is to regularize the problem: $\lambda$ is larger (resp. smaller) than the smallest (resp. largest) eigenvalue of $A$. So I don't think we can apply this Taylor expansion... Thanks for your help, I edited the question accordingly. $\endgroup$ – Mathieu Galtier Jan 8 '14 at 11:46
  • $\begingroup$ Hi Mathieu, to regularize the problem, addind a scalar matrix, is a bad method. It is much better to do as follows: use this robust Cholesky factorization: math.berkeley.edu/~cinnawu/hss.pdf $\endgroup$ – loup blanc Jan 8 '14 at 17:21
  • $\begingroup$ Unfortunately (for you), according your comment, your question is an absolute non-sense. Indeed, if your problem is ill-conditioned, then you know the matrix $L$ with a very poor precision ; then what is the interest to regularize the problem in these conditions ? $\endgroup$ – loup blanc Jan 8 '14 at 22:25
  • $\begingroup$ I don't think we mean the same thing by regularization. I am doing this in the framework of an adaptive algo where I can compute exactly $L$ in $\mathcal{O}(n^2)$. The problem is that it has eigenvalues very close to $0$. Therefore, when I want to invert it I get extremely large values for $x$. I am willing to loose a bit of precision in this inversion if it gives me smaller values for $x$. The regularization needs to be done, not for the Cholesky decomposition, but for the inversion. Thanks for your time. $\endgroup$ – Mathieu Galtier Jan 9 '14 at 9:56

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.