16

I have the Cholesky decomposition of a matrix $M$. However, I need the Cholesky decomposition of the inverse of the matrix, $M^{-1}$. Is there a fast way to do this, without first computing $M^{-1}$? In other words, is there a relationship between the Cholesky decompositions of a matrix and of its inverse?

My matrix is a covariance matrix and, hence, positive-definite.

4 Answers4

7

If $L^T L = R$ is the available Cholesky decomposition, then inverting both sides of the equation you get,

$$L^{-1}(L^{T})^{-1} = R^{-1} $$

And since transposition and inverse are interchangeable:

$$L^{-1}(L^{-1})^{T} = R^{-1} $$

So if you define $P = (L^{-1})^T$ this is your desired answer. In other words,

$$P^{T}P=R^{-1}$$

317070
  • 115
nimrodm
  • 236
  • 14
    The decomposition obtained here for $R^{-1}$ is not the Cholesky decomposition as the first factor ,$P^T$, is upper triangular (should be lower triangular). – Juho Kokkala Jan 29 '17 at 12:58
  • 1
    @JuhoKokkala Why does that matter? It is still triangular. – Rodrigo de Azevedo Jul 11 '18 at 17:33
  • 3
    Whether it matters in practice would depend on what you use the Cholesky factor for (not specified in the question). The question asks how to get the Cholesky decomposition of $M^{-1}$, and this answer gives some decomposition that (in the general case) is not the Cholesky decomposition. – Juho Kokkala Jul 11 '18 at 18:30
  • 2
    In general, is it numerically more stable to take inverse, then choleksy decomp; or cholesky decomp, then inverse? – seeker_after_truth Feb 04 '23 at 22:13
3

I do not believe there is any way to do it without some $O(n^3)$ operation. Different choices can have different constants. In particular, assuming $M=LL^T$ and you already have the Cholesky factor $L$, computing $L^{-1}$ costs ~$n^3/3$ flops, and given that, you can form $M^{-1}$ = $L^{-T}L^{-1}$ in ~$n^3/3$ flops and take its Cholesky in another ~$n^3/3$ flops, which is ~$n^3/3$ flops cheaper than pretending you don't have $L$, forming $M^{-1}$ directly and taking its Cholesky. But it's still $O(n^3)$.

  • Isn't inverting a triangular matrix O(n^2) instead of the O(n^3) of a matrix in general? Which would make the operation O(n^2) too, and thus sometimes worthwhile. – ninpnin Jun 19 '24 at 11:05
1

To add to previous answers, if we view $X$ as covariance matrix of data, the relationship between two decompositions reduces to relationship between coefficients of "right-to-left" autoregressive model and "left-to-right" one (details)

In general, switching order of "inverse" and "Cholesky" gives different results. IE

$$ \text{X=}\left( \begin{array}{ccc} 10 & 0 & 4 \\ 0 & 10 & -2 \\ 4 & -2 & 15 \\ \end{array} \right) $$

Using Cholesky and inverting factors, gives this decomposition

$$X^{-1}= \left( \begin{array}{ccc} 1 & 0 & -\frac{2}{5} \\ 0 & 1 & \frac{1}{5} \\ 0 & 0 & 1 \\ \end{array} \right)\left( \begin{array}{ccc} \frac{1}{10} & 0 & 0 \\ 0 & \frac{1}{10} & 0 \\ 0 & 0 & \frac{1}{13} \\ \end{array} \right)\left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -\frac{2}{5} & \frac{1}{5} & 1 \\ \end{array} \right) $$

While Cholesky on the inverse directly gives

$$X^{-1}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ -\frac{4}{73} & 1 & 0 \\ -\frac{20}{73} & \frac{2}{15} & 1 \\ \end{array} \right),\left( \begin{array}{ccc} \frac{73}{650} & 0 & 0 \\ 0 & \frac{15}{146} & 0 \\ 0 & 0 & \frac{1}{15} \\ \end{array} \right),\left( \begin{array}{ccc} 1 & -\frac{4}{73} & -\frac{20}{73} \\ 0 & 1 & \frac{2}{15} \\ 0 & 0 & 1 \\ \end{array} \right) $$

notebook

0

Minor nitpick: a covariance matrix can be positive semidefinite if it is not full rank.

Indeed conventionally a Cholesky factorization is defined in terms of $M = LL^T$ for lower triangular matrix $L$, or $M = U^TU$. I still think it is fair to consider $M=UU^T$ or $M=L^TL$ as Cholesky factorizations too, as the algorithm to calculate them are almost the same, the difference being only the order of the elements being solved. Either way, you would be doing back substitution when using these factors to solve a linear equation.

Back to the question: given a Cholesky factorization $M = LL^T$, we can find $L^{-1}$ which is a cheaper operation compared to general matrix inverse (although still on the order of $O(n^3)$. The problem is, as already discussed in other answers, $L^{-1}$ is also lower triangular, but $M^{-1} = L^{-T}L^{-1}$. The right hand side is of the form $UU^T$, and there is not of the conventional Cholesky factorization form.

How to convert a $U$ to an $L$? The answer is to use Givens rotation. Suppose the rotation to zero out $U_{1,2}$ is named $R_{1,2}$. $$ \hat U=UR_{1,2} $$ Because $R_{12}$ is orthogonal, we have $$ \hat U \hat U^T = UR_{1,2}R_{1,2}^TU^T = UU^T $$ We can consecutively apply other rotations to zero out the remaining off-diagonal upper triangular elements, which will result in an $L$.

$$ L = UR_{1,2}R_{1,3}\dots R_{n-1,n} $$

Each Givens rotation is an $O(n)$ operation, and there are $n(n-1)/2$ off-diagonal upper triangular elements. Therefore, the whole calculation is on the order of $O(n^3)$.


While not explicitly mentioned in the question, the answer by Yaroslav touched on $LDL^T$ factorization, so I will elaborate on that too. For $M=LDL^T$, for which $M^{-1}=L^{-T}D^{-1}L^{-1}$ and hence in the form of $UDU^T$, one cannot simply insert Givens rotation. Let $$ UDU^T = \hat U \hat D \hat U^T $$ From orthogonality of Givens rotation, we get $$ \hat D = R^TDR $$ but $\hat D$ may no longer be diagonal after this transformation. The correct approach would be to multiply the square roots of $D_i$ to $L$ so we can temporarily have $I$ as the diagonal. Then we apply Givens, and divide $L$ by the its resulting diagonal, while the new diagonal matrix will be its square. There is a concern that if $D$ contains negative elements (when M is symmetric indefinite), columns of $L$ will become imaginary. But this is temporary; after the division, $L$ will return back to real.

The following algorithm for converting $U$ to $L$ is found by algebraic manipulation with Sympy, so I don't have a rigorous proof of it. Assuming we process column-by-column from left to right, and for each column, row-by-row from top to bottom, in order to zero out $U_{i,j}$ where $i<j$, the new diagonal values are calculated as $$ \begin{aligned} D'_i = D_i + U^2_{i,j}D_j \\ D'_j = {D_iD_j \over D'_i} \end{aligned} $$ In the new matrix $U'$, the element at $(i,j)$ is transferred to $(j,i)$ as $$ U'_{j,i} = {U_{i,j}D_j\over D'_i} $$ At column $i$, all other nonzero elements below the diagonal are transformed as $$ U'_{k,i} = {U_{k,i}D_i + U_{i,j}U_{k,j}D_j\over D'_i} $$ and at column $j$, all other nonzero elements above the diagonal are transformed as $$ U'_{k,j} = U_{k,j}-U_{i,j}U_{k,i} $$ A caveat: this process will break down if $D_i'=0$.

syockit
  • 899