Search code examples
rmatrixnormal-distributioncovariance-matrix

The leading minor of order 14 is not positive (rmvnorm)


I have a 30 by 30 covariance matrix that is symmetric (passes with both the isSymmetric function as well as is.symmetric.matrix from matrixcalc library). However, I can't sample from it using rmvnorm(), as I get the following error:

Error in chol.default(Sigma, ...) : the leading minor of order 14 is not positive

Code to recreate my matrix (V) is below. I saw other people having a similar problem, but most of them realized their matrix was just barely asymmetric. Mine doesn't have that problem. Any fixes?

V <- structure(c(29.7886188034286, 0.000135608020306483, 3.31533243601536e-05, 
1.4149260534386e-08, 1.47707078131224e-10, 3.37580764203078e-07, 
5.64221265795885e-13, 2.25017993462975e-15, 2.29195223482087e-16, 
9.10148117311587e-22, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 0.000135608020306483, 29.7886188034286, 
2.2426518115228e-08, 5.5775762335415e-12, 1.77927684962432e-14, 
5.50469954608396e-09, 4.95218818529839e-17, 3.90846912833641e-11, 
1.50295585574595e-17, 1.16206427894487e-20, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.31533243601536e-05, 
2.2426518115228e-08, 29.7886188034286, 0.00685475942462371, 2.11134810790327e-05, 
0.0534401615284902, 5.72799389787849e-08, 4.18624110716291e-15, 
1.11496083204703e-10, 2.41081638713732e-17, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.4149260534386e-08, 
5.5775762335415e-12, 0.00685475942462371, 29.7886188034286, 0.0485441146570834, 
0.00235235999633464, 0.000151744662859306, 1.92911292294946e-17, 
1.2344146127685e-08, 9.27870135119539e-17, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.47707078131224e-10, 
1.77927684962432e-14, 2.11134810790327e-05, 0.0485441146570834, 
29.7886188034286, 4.12378769436744e-06, 0.077752722328188, 3.69961837029023e-20, 
2.24120473537269e-09, 2.78180051438489e-18, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.37580764203078e-07, 
5.50469954608396e-09, 0.0534401615284902, 0.00235235999633464, 
4.12378769436744e-06, 29.7886188034286, 1.53888416023551e-08, 
2.27088584364582e-13, 1.90414928862951e-08, 1.27417002165266e-14, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5.64221265795885e-13, 
4.95218818529839e-17, 5.72799389787849e-08, 0.000151744662859306, 
0.077752722328188, 1.53888416023551e-08, 29.7886188034286, 1.87089710144973e-22, 
6.07333453807195e-10, 2.7834529260172e-19, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.25017993462975e-15, 
3.90846912833641e-11, 4.18624110716291e-15, 1.92911292294946e-17, 
3.69961837029023e-20, 2.27088584364582e-13, 1.87089710144973e-22, 
29.7886188034286, 1.09623746831197e-16, 7.73956371456588e-14, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.29195223482087e-16, 
1.50295585574595e-17, 1.11496083204703e-10, 1.2344146127685e-08, 
2.24120473537269e-09, 1.90414928862951e-08, 6.07333453807195e-10, 
1.09623746831197e-16, 29.7886188034286, 9.14354495463089e-09, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 9.10148117311587e-22, 
1.16206427894487e-20, 2.41081638713732e-17, 9.27870135119539e-17, 
2.78180051438489e-18, 1.27417002165266e-14, 2.7834529260172e-19, 
7.73956371456588e-14, 9.14354495463089e-09, 29.7886188034286, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1.1, 2.23435064143926e-11, 1.33547219558487e-12, 
2.43247400092917e-19, 2.65083889750281e-23, 1.38463754600901e-16, 
3.86793905062364e-28, 6.15198425841414e-33, 6.38251501624773e-35, 
1.00647935339086e-45, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 2.23435064143926e-11, 1.1, 6.1108893902982e-19, 
3.77982549438626e-26, 3.84651806542874e-31, 3.68169628548815e-20, 
2.9797168689254e-36, 1.85606889272712e-24, 2.74456278291983e-37, 
1.64074474950797e-43, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1.33547219558487e-12, 6.1108893902982e-19, 
1.1, 5.70897210323724e-08, 5.41627098503697e-13, 3.46944376178704e-06, 
3.98644533783302e-18, 2.1292624542656e-32, 1.51042770288842e-23, 
7.06169554938132e-37, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 2.43247400092917e-19, 3.77982549438626e-26, 
5.70897210323724e-08, 1.1, 2.86287682950135e-06, 6.72335136471411e-09, 
2.79774037428088e-11, 4.52164040370437e-37, 1.85141095238044e-19, 
1.04605642973501e-35, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 2.65083889750281e-23, 3.84651806542874e-31, 
5.41627098503697e-13, 2.86287682950135e-06, 1.1, 2.06620454498461e-14, 
7.3439529070571e-06, 1.66300890475797e-42, 6.10300580448813e-21, 
9.40226930623985e-39, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1.38463754600901e-16, 3.68169628548815e-20, 
3.46944376178704e-06, 6.72335136471411e-09, 2.06620454498461e-14, 
1.1, 2.87734924851621e-19, 6.26572494547037e-29, 4.4053732445802e-19, 
1.9725839084239e-31, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 3.86793905062364e-28, 2.9797168689254e-36, 
3.98644533783302e-18, 2.79774037428088e-11, 7.3439529070571e-06, 
2.87734924851621e-19, 1.1, 4.25285449747533e-47, 4.48162101889993e-22, 
9.41344266929708e-41, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 6.15198425841414e-33, 1.85606889272712e-24, 
2.1292624542656e-32, 4.52164040370437e-37, 1.66300890475797e-42, 
6.26572494547037e-29, 4.25285449747533e-47, 1.1, 1.46012488822639e-35, 
7.27802729314438e-30, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 6.38251501624773e-35, 2.74456278291983e-37, 
1.51042770288842e-23, 1.85141095238044e-19, 6.10300580448813e-21, 
4.4053732445802e-19, 4.48162101889993e-22, 1.46012488822639e-35, 
1.1, 1.01580402447928e-19, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1.00647935339086e-45, 1.64074474950797e-43, 
7.06169554938132e-37, 1.04605642973501e-35, 9.40226930623985e-39, 
1.9725839084239e-31, 9.41344266929708e-41, 7.27802729314438e-30, 
1.01580402447928e-19, 1.1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.1, 2.23435064143926e-11, 
1.33547219558487e-12, 2.43247400092917e-19, 2.65083889750281e-23, 
1.38463754600901e-16, 3.86793905062364e-28, 6.15198425841414e-33, 
6.38251501624773e-35, 1.00647935339086e-45, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.23435064143926e-11, 
1.1, 6.1108893902982e-19, 3.77982549438626e-26, 3.84651806542874e-31, 
3.68169628548815e-20, 2.9797168689254e-36, 1.85606889272712e-24, 
2.74456278291983e-37, 1.64074474950797e-43, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.33547219558487e-12, 
6.1108893902982e-19, 1.1, 5.70897210323724e-08, 5.41627098503697e-13, 
3.46944376178704e-06, 3.98644533783302e-18, 2.1292624542656e-32, 
1.51042770288842e-23, 7.06169554938132e-37, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.43247400092917e-19, 
3.77982549438626e-26, 5.70897210323724e-08, 1.1, 2.86287682950135e-06, 
6.72335136471411e-09, 2.79774037428088e-11, 4.52164040370437e-37, 
1.85141095238044e-19, 1.04605642973501e-35, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2.65083889750281e-23, 
3.84651806542874e-31, 5.41627098503697e-13, 2.86287682950135e-06, 
1.1, 2.06620454498461e-14, 7.3439529070571e-06, 1.66300890475797e-42, 
6.10300580448813e-21, 9.40226930623985e-39, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.38463754600901e-16, 
3.68169628548815e-20, 3.46944376178704e-06, 6.72335136471411e-09, 
2.06620454498461e-14, 1.1, 2.87734924851621e-19, 6.26572494547037e-29, 
4.4053732445802e-19, 1.9725839084239e-31, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.86793905062364e-28, 
2.9797168689254e-36, 3.98644533783302e-18, 2.79774037428088e-11, 
7.3439529070571e-06, 2.87734924851621e-19, 1.1, 4.25285449747533e-47, 
4.48162101889993e-22, 9.41344266929708e-41, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6.15198425841414e-33, 
1.85606889272712e-24, 2.1292624542656e-32, 4.52164040370437e-37, 
1.66300890475797e-42, 6.26572494547037e-29, 4.25285449747533e-47, 
1.1, 1.46012488822639e-35, 7.27802729314438e-30, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6.38251501624773e-35, 
2.74456278291983e-37, 1.51042770288842e-23, 1.85141095238044e-19, 
6.10300580448813e-21, 4.4053732445802e-19, 4.48162101889993e-22, 
1.46012488822639e-35, 1.1, 1.01580402447928e-19, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.00647935339086e-45, 
1.64074474950797e-43, 7.06169554938132e-37, 1.04605642973501e-35, 
9.40226930623985e-39, 1.9725839084239e-31, 9.41344266929708e-41, 
7.27802729314438e-30, 1.01580402447928e-19, 1.1), dim = c(30L, 
30L))
rmvnorm(1, sigma = V)

Solution

  • Your problem is not with lack of symmetry, but rather with lack of positive definiteness. This also isn't the usual case where the matrix is almost positive definite (i.e. the smallest eigenvalues are near zero, either positive or negative); the smallest eigenvalue is decidedly negative (-8.9, compared to a dominant eigenvalue of 37.4).

    eigen(V)$values
    [1] 37.424265 29.861301 29.823050 29.788700 29.788619 29.788619 29.788483
     [8] 29.788059 29.735071 29.697642  3.502381  1.100006  1.100006  1.100002
    [15]  1.100002  1.100000  1.100000  1.100000  1.100000  1.100000  1.100000
    [22]  1.100000  1.100000  1.100000  1.100000  1.099997  1.099997  1.099992
    [29]  1.099992 -8.899997
    

    On my system (mvtnorm version 1.2.5), mvtnorm::rmvnorm will let me sample MVN values with this covariance matrix, but after warning ("sigma is numerically not positive semidefinite") it truncates all components corresponding to negative eigenvalues to zero. It is up to you to decide whether this makes sense for your application or not (!!) (This is with the default method = "eigen"; it does similar things (SVD, Cholesky with pivoting) to ignore non-pos-def components when using other methods.)

    For what it's worth there are a lot of rmvnorm() functions in the R ecosystem — these are just packages that I have installed on my system that contain an rmvnorm() function. They will all deal with non-positive-definite covariance matrices in different ways ...

    unique(help.search("rmvnorm", agrep = FALSE)$matches$Package)
    ## [1] "CDM"             "dae"             "dqrng"           "geoR"           
    ## [5] "ks"              "mhsmm"           "miceadds"        "mixtools"       
    ## [9] "mvtnorm"         "rethinking"      "Rfast"           "spam"           
    ## [13] "splus2R"         "stockassessment"