CODEWITHSUNDEEP
rmatrixoptimizationlinear-algebramatrix-inverse
Euler_Salter
Euler_Salter

Reputation: 3561

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

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?

Upvotes: 4

Views: 931

Answers (2)

G. Grothendieck
G. Grothendieck

Reputation: 269481

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

Upvotes: 3

Euler_Salter
Euler_Salter

Reputation: 3561

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

system

Writing this using the cholesky decomposition we have

cholesky

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

substitution

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.

finalsystem

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.

Upvotes: 1

Related Questions