# Best way to compute $A^{-1}$ when the Cholesky decomposition $A=LL^T$ is known

Suppose $$\mathbf{A}$$ is symmetric positive definite, and that I have available the Cholesky decomposition of $$\mathbf{A}=\mathbf{L}_A\mathbf{L}_A^T$$. I want to know $$\mathbf{A}^{-1}$$. Which of the two methods below are best (fastest, robust, or any other quality indicators)?

1. Compute $$\mathbf{L}_A^{-1}$$ and then $$\mathbf{A}^{-1}=\mathbf{L}_A^{-T}\mathbf{L}_A^{-1}$$. EDIT: Here, I re-use $$\mathbf{L}_A$$ from a previous step in my calculations, so assume no computational cost of obtaining $$\mathbf{L}_A$$.
2. Compute $$\mathbf{A}^{-1}$$ directly. E.g. MATLAB applies first the LU decomposition and then uses the results to form a linear system whose solution is $$\mathbf{A}^{-1}$$ https://se.mathworks.com/help/matlab/ref/inv.html#d123e769694.
• Tangential question, but why would you want to compute $A^{-1}$? May 3, 2022 at 9:25
• I have $L_A$ from a previous section in my algorithm, and in the next part I need $A^{-1}$. Since $L_A$ is lower triangluar, I was wondering if it is possible to exploit this structure when calculating $L_A^{-1}$. May 5, 2022 at 10:26
• Why does the next part of your algorithm need $A^{-1}$ though? A rule of thumb in numerical linear algebra is that one usually should not form $A^{-1}$ explicitly. May 5, 2022 at 13:54
• In my case, $A^{-1}$ is used as a weight in a subsequent optimization problem. I have $J=x^T(B+A^{-1})x -2(C+dA^{-1})x$ and I want to minimize $J$ wrt $x$. Here, $B$ and $C$ are given matrices and $d$ is a given vector. I don't see how I can avoid calculating $A^{-1}$ explicitly in this case, but if you have an idea that would be great! May 10, 2022 at 8:45

Method (1). An efficient method of inverting a lower triangular matrix ($$L$$, in this case) requires $$< n^3/2 + n^2/2$$ operations.[1] Computing the product $$L^{-T}L^{-1}$$ requires $$2n^3/3+n/3$$ operations using an efficient method for multiplying triangular matrices.[2]

Method (2). Since you already have Cholesky factors, you should not perform LU decomposition in order to compute $$A^{-1}$$; because doing so, would effectively mean that you are deriving Cholesky factors again! This only incurs unnecessary cost. In this case, you can solve two triangular systems (each would cost $$n^3/2 +n^2/2$$ arithmetic operations) using the well-known method:

$$\text{Data: } L, L^T. \text{ Output: }X=A^{-1}$$ \begin{align*} AX &= I \\ L(L^TX) &= I \end{align*} $$$$\left\{ \begin{array}{@{}ll@{}} i)\ LY=I & \rightarrow Y\ ✔ \\ ii)\ L^{T}X = Y & \rightarrow X=A^{-1}\ ✔ \end{array}\right.$$$$

Summing the operations count, you'll find that the second method requires less computation if properly implemented. If you are using MATLAB built-in functions, I believe the following command would be fastest way to compute $$A^{-1}$$:

X = L' \ (L \ eye(n));


As a final note, I would like to give a quote from Meyer[3] discussing the same problem:

A tempting alternate solution might be to use the fact $$A^{-1} = (LU)^{−1} = U^{−1}L^{−1}$$. But computing $$U^{−1}$$ and $$L^{−1}$$ explicitly and then multiplying the results is not as computationally efficient as the method just described.

EDIT:

Below, is a naïve performance comparison of the two methods in MATLAB for randomly generated s.p.d matrices. As you can see, on average, you will get slightly faster results with method (2).

[1] Stewart, Gilbert W. Matrix algorithms: volume 1: basic decompositions. Society for Industrial and Applied Mathematics, 1998, chapter 2, section 2, algorithm 2.3.

[2] Lyche, Tom. Numerical linear algebra and matrix factorizations. Vol. 22. Springer Nature, 2020..

[3] Meyer, Carl D. Matrix analysis and applied linear algebra. Vol. 71. Siam, 2000.

From the book of Nick Higham on "Functions of Matrices Theory and Computation", we can find the theoretical number of flops for both LU and Cholesky factorization where $$n$$ is the A matrix size $$n \times n$$.

In general, Cholesky $$\mathcal{O}(n^3/3)$$ is twice as fast as LU decomposition $$\mathcal{O}(2n^3/3)$$ but depending on your matrix, LU can be made faster in the order magnitude of $$\mathcal{O}(n^2)$$.

• Thanks for the reply. However, I am interested in $A^{-1}$ and not necessarily it's decompositions. Is it faster to compute $L_A^{-1}$ than directly $A^{-1}$? I was thinking that the lower triangular structure of $L_A$ would make it easier somehow to calulate $L_A^{-1}$ than directly $A^{-1}$. May 5, 2022 at 10:32
• You get the inverse $A^{-1}$ by the inverse of your decomposed matrices. Well, the given table shows you the complete computation cost of the decomposed matrices of LU and Cholesky decompositions. If you want to compute $L_A^{-1}$ from $L_A$ to get $A^{-1}$, the computation cost would still be the same at the end since the other method also computes the inverse of its decomposed matrices. May 5, 2022 at 10:45
• Yes, but since I already have $L_A$ from a previous step in my calculations then I avoid the cost of doing the decomposition of $A$ again, right? At least in MATLAB, it seems that decomposition would be some kind of LU-decomposition of $A$, and in my specific case that extra decomposition would probably be redundant. So based on what you said, it seems to me that it is cheaper to re-use $L_A$ (assume no cost) to calculate $L_{A}^{-1}$ and then $A^{-1}$ than calculating $A^{-1}$ "directly" (where "directly" means doing first a LU-decomposition and then inverting the decomposed matrices) May 10, 2022 at 8:58
• I edited the question slightly to emphasize that there is no cost in re-using $L_A$. May 10, 2022 at 9:04

I know this bit old, but for future readers, you said you need $$A^{-1}$$ in order to compute $$J=x^T(B+A^{-1})x -2(C+dA^{-1})x$$. What you should be doing is write the parts with inverses as matrix-vector products e.g. $$x^T(A^{-1}x)$$, where you first compute $$A^{-1}x=z$$ by solving the system $$Az=x$$ (in the way @S4JJ4D but with different RHS), which is cheap ($$O(n^2)$$ to solve as you have the cholesky factors given). Then you finally compute $$x^T z$$.