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

Question

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.

Answer 1

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*} \begin{equation} \left\{ \begin{array}{@{}ll@{}} i)\ LY=I & \rightarrow Y\ ✔ \\ ii)\ L^{T}X = Y & \rightarrow X=A^{-1}\ ✔ \end{array}\right. \end{equation}

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.}}

Answer 2

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)$.