R Inconsistent product of matrix inversions

Question

I have a matrix emat generated by taking the sum of outer products of vectors. It should be symmetric and positive definite. I am finding that

solve(emat) %*% solve(emat)

generates a different result from

bmat <- solve(emat)
bmat %*% t(bmat)

In fact, the two differ quite substantially (printed output of emat is shortened by R).

> emat
        V1        V2        V3
1 170.2939  15.77391 110.75499
2  15.7739 444.57862   8.87082
3 110.7550   8.87082  72.03669
> solve(emat) %*% solve(emat) - bmat %*% t(bmat)
       1         2 3
V1 -1024  48.00000 0
V2     8  -0.21875 0
V3  2048 -72.00000 0

But this should not be the case.

Is this a bug? Or is it an issue with precision? Or does it have to do with how R handles matrices stored in memory?

https://www.rdocumentation.org/packages/base/versions/3.5.1/topics/solve — GordonShumway, Oct 24 '18 at 03:33

Sam S. · Answer 1 · 2018-10-24T04:08:14.387

1

The problem could be from your emat matrix. I tried it, and it worked for me.

emat <- matrix(c(2,-1,0,-1,2,-1,0,-1,2),3,3)

# check your matrix to see if it is positive definite matrix or not
library(matrixcalc)
is.positive.definite(emat)

bmat <- solve(emat)

# the result of the following is zero matrix
solve(emat) %*% solve(emat) - bmat %*% t(bmat)

If your matrix is not symmetric (so not positive definite), the result would not be zero as emat is not equal to t(emat).

edited Oct 24 '18 at 04:08

answered Oct 24 '18 at 03:45

Sam S.

627
1
7
23

Thanks! Based on your example, as well as the response by Daniel Fischer, it is likely that this is an issue to do with precision. But given that `emat` was generated by taking the sum of outer products of vectors, it is strange that it is not symmetric. That would imply that there is greater precision for the k-th row than for the k-th column, or the other way around, which would be strange. – jkcshea Oct 24 '18 at 14:32
It would be precision problem. You can also check with is.positive.definite(emat) to see if it is positive definite (which should be also symmetric) or not. Please also check the below link which says the product is positive semidefinite. https://math.stackexchange.com/questions/2136242/is-the-outer-product-of-a-column-vector-with-itself-positive-semi-definite – Sam S. Oct 28 '18 at 07:47

score 1 · Answer 2 · answered Oct 24 '18 at 05:18

1

Your emat matrix is not symmetric, see e.g. emat[1,2] != emat[2,1]

Try to run

emat <- matrix(c(170.2939, 15.77391, 110.75499, 15.7739, 444.57862, 8.87082, 110.7550, 8.87082, 72.03669),ncol=3)
emat <- round(emat,4)
bmat <- solve(emat)

solve(emat) %*% solve(emat) - bmat %*% t(bmat)

where the rounding takes care that your data is symmetric (in this particular example, as your data is close enough to be the same after rounding...).

answered Oct 24 '18 at 05:18

Daniel Fischer

3,280
1
18
29

Thanks so much for trying that example. I apologize for not mentioning that the output above was shortened by R (there are many, many more decimal places). `emat` was generated by taking the sum of outer products of vectors, so it should indeed be symmetric and positive definite. I'll update the question with this information. So am I dealing with an issue related to precision? – jkcshea Oct 24 '18 at 09:30

R Inconsistent product of matrix inversions

2 Answers2