Error in R function `is.positive.definite`

Question

I want to test whether a matrix is positive definite or not in R. I used the R function is.positive.definite but kept getting the following error message although my matrix is symmetric as function isSymmetric shows. Is this because of rounding errors, please?

Error in is.positive.definite(S) : argument x is not a symmetric matrix

My working code is attached below. Could anyone help me with this, please? Thanks.

library(Matrix) # isSymmetric
library(matrixcalc) # is.positive.definite
library(expm) # sqrtm

###################################################################################################

theta0 <- c(0.2, 10)

###################################################################################################

OS.mean <- function(shape, rank, n=10){
  term1 <- factorial(n)/(factorial(rank-1)*factorial(n-rank))
  term2 <- beta(n-rank+1, rank) - beta(n-rank+shape+1, rank)
  term1*term2/shape
}

OS.mean.theta0.10 <- as.matrix(OS.mean(theta0[1], rank=seq(1, 10, by=1)))

###################################################################################################

OSsq.mean <- function(shape, rank, n=10){
  term1 <- factorial(n)/(factorial(rank-1)*factorial(n-rank))
  term2 <- beta(n-rank+1, rank) - 2*beta(n-rank+shape+1, rank) + beta(n-rank+2*shape+1, rank)
  term1*term2/(shape*shape)
}

OSsq.mean.theta0.10 <- as.matrix(OSsq.mean(theta0[1], rank=seq(1, 10, by=1)))

###################################################################################################

OSprod.mean <- function(shape, rank1, rank2, n=10){
  term1 <- factorial(n)/(factorial(rank1-1)*factorial(rank2-rank1-1)*factorial(n-rank2))
  term2 <- beta(n-rank1+1, rank1) - beta(n-rank1+shape+1, rank1)
  term3 <- beta(n-rank2+1, rank2-rank1)
  term4 <- beta(n-rank1+shape+1, rank1) - beta(n-rank1+2*shape+1, rank1)
  term5 <- beta(n-rank2+shape+1, rank2-rank1)
  term1*(term2*term3-term4*term5)/(shape*shape)
}

OS.cov <- function(shape, rank1, rank2){
  OSprod.mean(shape, rank1, rank2) - OS.mean(shape, rank1)*OS.mean(shape, rank2)  
}

###################################################################################################

spacing <- seq(1, 10, by=1)

OS.varcov.10 <- function(shape, n=10){
  V.diag <- diag(c(OSsq.mean.theta0.10 - OS.mean.theta0.10^2))

  V.upper <- matrix(0, nrow=10, ncol=10)
  for(i in 1:9){
    for(j in (i+1):10){
      V.upper[i, j] <- OS.cov(shape, spacing[i], spacing[j])
    }
  }

  V.upper + V.diag + t(V.upper)
}

###################################################################################################

V.theta0.10 <- OS.varcov.10(theta0[1])

kappa(V.theta0.10)

isSymmetric(V.theta0.10)
is.positive.definite(V.theta0.10)

S <- sqrtm(V.theta0.10)
isSymmetric(S)
is.positive.definite(S)

How does this differ from [your previous question](http://stackoverflow.com/questions/31643031/covariance-matrix-is-no-longer-positive-definite)? — josliber, Jul 27 '15 at 05:32
`is.positive.definite` function uses `is.symmetric.matrix` function, not `isSymmetric` function to test whether a matrix is a numeric symmetric square matrix. — , Jul 27 '15 at 05:37

cyberj0g · Accepted Answer · 2015-08-04T11:39:03.717

7

Your S matrix is not symmetric because of loss of significance, but default number of output decimal places hides it. See for yourself:

> options(digits=20)
> S[1,2]
[1] 0.033457660484940172
> S[2,1]
[1] 0.033457660484940213

The thing is, is.symmetric.matrix from matrixcalc package doesn't account for small differences (i.e. it just compares matrix elements with strict == instead of all.equal method) while isSymmetric from Matrix package does. If you round the matrix, everything will be fine:

> S=round(S,10)
> is.symmetric.matrix(S)
[1] TRUE
> is.positive.definite(S)
[1] TRUE

edited Aug 04 '15 at 11:39

answered Jul 27 '15 at 05:38

cyberj0g

3,707
1
19
34

It is known that a positive definite matrix has a Unique Positive Definite square root. This is calculated by `sqrtm` function. That is, `S` is supposed to be positive definite in theory. However, it is not here. Is it because of rounding error, please? – LaTeXFan Jul 27 '15 at 05:42
Yes, it's because of [loss of significance](https://en.wikipedia.org/wiki/Loss_of_significance). `is.symmetric.matrix` from `matrixcalc` package doesn't account for small differences while `isSymmetric` from `Matrix` package does. All you need to do is round the matrix. – cyberj0g Jul 27 '15 at 05:48
Can arbitrary precision type program solve this rounding problem, please? – LaTeXFan Jul 27 '15 at 05:53
Why just not add `S=round(S,10)`? – cyberj0g Jul 27 '15 at 05:58

Error in R function `is.positive.definite`

1 Answers1