sqrtm vs chol to produce draws the multivariate normal distribution

Hi,

I am wondering if it safe to use sqrtm to produce draws from the multivariate normal distribution. Usually, one uses cholesky to make the draws from the distribution with the variance-covariance matrix of sigma. However, if the matrix sigma is not positive definite due to some approximation error, I would like to use sqrtm instead of cholesky to take the draws. Do you happen to know if it is safe to do so?

My function looks like this

function out=cholx(x)
[R,p] = chol(x);
if p==0
    out=R;
else
    out=real(sqrtm(x))';
end

And I would take the draws in a usual fashion:

cholx(sigggma)' * randn(N,1)

Or should I just use cholcov function?

Thanks for your help!

[zym]

Based on the simulation for the positive definite varianace-covariance matrix, it seems sqrtm works well (see below).

clear;
clc;
REPS=10000000;
N=10;
A=randn(N,N);
A=A*A';

baba=randn(N,REPS);
a1=chol(A)'*baba;
a2=real(sqrtm(A))*baba;

A
cov(a1')
cov(a2')

But it is unclear how would it perform if it is close to positive definite.

Answer 1

[John D'Errico]

Don't bother. sqrtm also can have a similar problem, but you won't like what it does either! For example:

A = rand(3,2);
A = A*A' - diag(rand(3,1)*1e-15)

A =
0.92657      0.47991      0.30152
0.47991        0.502      0.36037
0.30152      0.36037      0.26265

This insures that A is not positive definite. As I initially built is, A has rank only 2. But then I subtracted a tiny amount off the diagonal. That will kill any chances chol has, even though the matrix will still be symmetric.

chol(A)
Error using chol
Matrix must be positive definite. 

sqrtm(A)
ans =
0.90255 +  1.337e-10i      0.29619 - 1.1723e-09i      0.15571 + 1.4549e-09i
0.29619 - 1.1723e-09i      0.51422 + 1.0278e-08i       0.3871 - 1.2756e-08i
0.15571 + 1.4549e-09i       0.3871 - 1.2756e-08i      0.29759 + 1.5832e-08i

As you can see, sqrtm produces complex results, and you won't be happy with what that does to your random number draws. eig tells the story here.

eig(A)
ans =
-5.3877e-16
0.27483
1.4164

So one tiny negative eigenvalue.

Instead, make the matrix SPD. And that means to just use my nearestSPD function. It insures chol will survive, because it tweaks the matrix minimally in those least significant bits to find a matrix that chol can use. (The final test in nearestSPD is that chol produces a result.)

Ahat = nearestSPD(A)
Ahat =
0.92657      0.47991      0.30152
0.47991        0.502      0.36037
0.30152      0.36037      0.26265

Ahat - A
ans =
-1.1102e-16  -1.1102e-16            0
-1.1102e-16   1.1102e-16  -3.8858e-16
0  -3.8858e-16   3.3307e-16

So the perturbations are infinitessimal. And as you can see, chol now works.

chol(Ahat)
ans =
0.96259      0.49856      0.31324
0      0.50343      0.40562
0            0   5.4505e-09

Oh, by the way, avoid using sqrtm on semi-definte matrices, even if chol survives, at least to generate random numbers. For example:

B = rand(3,2); B = B*B';

sqrtm(B)
ans =
      0.66778 + 1.4162e-09i      0.18439 + 7.2703e-10i      0.42238 - 2.5564e-09i
      0.18439 + 7.2703e-10i      0.59319 + 3.7324e-10i      0.27085 - 1.3124e-09i
      0.42238 - 2.5564e-09i      0.27085 - 1.3124e-09i      0.31102 + 4.6145e-09i
      
chol(B)
ans =
      0.81137      0.42756      0.57109
            0      0.52577      0.14953
            0            0   4.7521e-09

So chol actually survived this time. I did not subtract off that random crap on the diagonal. B is still rank 2 though. And sqrtm is not doing as well. Those complex entries will give you now complex random numbers, if you use sqrtm.

Find nearestSPD on the File Exchange.

https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd?s_tid=ta_fx_results

[zym]

Wow! Thanks a lot for this great answer and explanations.

Exactly, I was doubting if I should use "nearestSPD" or sqrtm() for my problem, where I am doing a lot of drawing. Sometimes I find "nearestSPD" to take a long time to find reasonable matrix, so I thought that sqrtm can be a faster solution. But given your answer, I should stick with "nearestSPD" since it is robust, even though it just take more time.

I tweaked the function to be a bit faster for my application and not to get stuck

function [Xf,R]=nspd1(A,time2run)

B=(A+A')/2;

[U,S,V]=svd(B);

H=V*S*V';

Xf=(B+H)/2;

Xf=(Xf+Xf')/2;

[R,test2]=chol(Xf,'lower');

if test2<1

elseif test2>=1

n=1;

a1=tic;

while test2>=1 && toc(a1)<time2run

mineig=min(eig(Xf));

Xf=Xf+(-mineig*n.^2+eps(mineig)+0.000000000000000001)*eye(size(Xf));

[R,test2]=chol(Xf,'lower');

n=n+1;

end

end

Answer 2

[the cyclist]

A common reason for not having a positive semidefinite correlation matrix is when the pairwise correlations are estimated (or obtained from different external sources), without consideration of whether they are coherent as system.

For example

sigggma = [1.00 0.99 0.75;
           0.75 1.00 0.10;
           0.99 0.10 1.00];

shows 1 & 2 highly correlated, 1 & 3 quite correlated, but 2 & 3 with minimal correlation. This is simply not possible (but pairwise estimates could result in this). Mine is a pretty exaggerated example, where the matrix is quite "far away" from being positive semidefinite. chol() will reject this, as you understand.

I'm not sure what exactly you mean by "safe". Your method will give an answer, but the correlation structure will not be what the input "requested" -- it's not possible! (See how different with my code below, for this case.)

@John D'Errico posted his very nice solution while I was composing and tweaking mine, and injects a notion of "minimally different" perturbation. The important thing is for you to understand how/why your randomly drawn numbers differ from what your input matrix is requesting.

rng default

NTRIALS = 5000;

sigggma = [1.00 0.99 0.75;
           0.75 1.00 0.10;
           0.99 0.10 1.00]; 

x = nan(NTRIALS,3);

for nt = 1:NTRIALS
    
    x(nt,:) = cholx(sigggma)' * randn(3,1);

end

corrcoef(x)

ans = 3×3

1.0000 0.8300 0.6848 0.8300 1.0000 0.1619 0.6848 0.1619 1.0000

function out=cholx(x)
    out=real(sqrtm(x))';
end

[the cyclist]

Just FYI, here are the empirical results I got from your cholx and @John D'Errico's nearestSPD, for the input

sigggma = [1.00 0.99 0.75;
           0.75 1.00 0.10;
           0.99 0.10 1.00];

cholx

    1.0000    0.8254    0.6871
    0.8254    1.0000    0.1569
    0.6871    0.1569    1.0000


nearestSPD

    1.0000    0.8199    0.8178
    0.8199    1.0000    0.5045
    0.8178    0.5045    1.0000

[zym]

Thanks a lot for the answer! It helped me to understand this better