Why is cholesky decomposition not giving me the same result as simply inverting the matrix?

Question

Setting to reproduce my Minimal Working Example

I have the following matrix

K <- matrix(c(1.250000e+00, 3.366892e-07, 4.641930e-10, 1.455736e-08, 1.049863e-06, 
              3.366892e-07, 1.250000e+00, 5.482775e-01, 8.606555e-01, 9.776887e-01,
              4.641930e-10, 5.482775e-01, 1.250000e+00, 8.603413e-01, 4.246732e-01,
              1.455736e-08, 8.606555e-01, 8.603413e-01, 1.250000e+00, 7.490100e-01,
              1.049863e-06, 9.776887e-01, 4.246732e-01, 7.490100e-01, 1.250000e+00), nrow=5)

and the following vector

y <- matrix(c(39.13892, 12.34428, 12.38426, 14.71951, 10.52160), nrow=5)

Problem

I want to calculate the product between the inverse of K and the vector y.

Naive Method - It works

The naive method works (I have a way of checking, but doesn't matter here)

solve(K) %*% y

           [,1]
[1,] 31.3111308
[2,]  3.0620869
[3,]  3.7383357
[4,]  6.6257060
[5,]  0.7820081

Cholesky Decomposition - Fails

However the "clever" method fails. I use cholesky decomposition which gives me an upper triangular matrix. Then I solve the system L z = y by backward substitution and the system L^T x = L^{-1} y by forward substitution.

L <- chol(K)  ## upper triangular
forwardsolve(t(L), backsolve(L, y))

          [,1]
[1,]  31.31112
[2,] -14.16259
[3,]   9.84534
[4,]  39.67900
[5,]  33.54842

What is happening? This matrix K and this vector 'y` are just an example. It happens with many other similar vectors and matrices. Why?

@Matt incredibly, this gives an answer that is different from both! `[,1] [1,] 39.138913 [2,] -14.385526 [3,] -5.993738 [4,] 27.188652 [5,] 6.766934` — Euler_Salter, Dec 23 '19 at 18:35
@Matt I think someone had a similar problem [here](https://stats.stackexchange.com/questions/147654/cholesky-factorization-and-forward-substitution-less-accurate-than-inversion) but none really managed to answer — Euler_Salter, Dec 23 '19 at 18:36
`chol2inv(chol(K)) %*% y` gives the results as `solve(K) %*% y` for me. See [this](https://stackoverflow.com/questions/26718334/comparing-matrix-inversions-in-r-what-is-wrong-with-the-cholesky-method) — Matt, Dec 23 '19 at 18:39
@Matt that’s right, same result! do you know if chol2inv should be preferred to solving the system using forwardsolve and backsolve? — Euler_Salter, Dec 23 '19 at 19:42
@Matt I think I managed to find a way around it, not sure why it works though. I found out that doing `backsolve(L, forwardsolve(t(myL), ytrain))` actually works.. — Euler_Salter, Dec 23 '19 at 20:08
I meant `backsolve(L, forwardsolve(t(L), ytrain))` of course — Euler_Salter, Dec 23 '19 at 20:18

G. Grothendieck · Accepted Answer · 2019-12-23T21:08:40.807

The key point is that when taking the inverse of a product one must reverse the product of the inverses:

solve(A %*% B) = solve(B) %*% solve(A)

In the question the order was not reversed.

If R = chol(K) where we use R to emphasize that it is upper right triangular then:

solve(K, y)
= solve(t(R) %*% R, y)   since K = t(R) %*% R
= solve(t(R) %*% R) %*% y
= solve(R) %*% solve(t(R)) %*% y  note that we have reversed the order
= solve(R) %*% solve(t(R), y)
= backsolve(R, forwardsolve(t(R), y))

In the last line we used the fact that the transpose of R is lower left triangular and forwardsolve works on such matrices whereas backsolve works on upper right triangular matrices.

We can check that this does give the same answer as using solve direclty:

R = chol(K)
all.equal(backsolve(R, forwardsolve(t(R), y)), solve(K, y))
# [1] TRUE

score 1 · Answer 2 · answered Dec 23 '19 at 20:16

This is just a partial answer, but it's better than nothing. Basically finding K^{-1}y is equivalent to solving the following system

Writing this using the cholesky decomposition we have

Basically we now consider Lz to be our variable first, call it x

and since L is upper-triangular, it's transpose is lower-triangular, so that we can use forwardsolve to find x

forwardsolve(t(L), y)

Once we've found x, we can just remember what it means, i.e.

so that we can use 'backsolveto findz`

backsolve(L, forwardsolve(t(L), y))

and this gives the correct answer. Not sure why the other way around does't work as well.

Why is cholesky decomposition not giving me the same result as simply inverting the matrix?

Setting to reproduce my Minimal Working Example

Problem

Naive Method - It works

Cholesky Decomposition - Fails

2 Answers2