sqrtm vs chol to produce draws the multivariate normal distribution

30 views (last 30 days)
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!
  1 Comment
zym
zym on 3 Mar 2023
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.

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 25 Mar 2023
Edited: John D'Errico on 25 Mar 2023
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.
  1 Comment
zym
zym on 25 Mar 2023
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

Sign in to comment.

More Answers (1)

the cyclist
the cyclist on 25 Mar 2023
Edited: the cyclist on 25 Mar 2023
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
  2 Comments
the cyclist
the cyclist on 25 Mar 2023
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
zym on 25 Mar 2023
Thanks a lot for the answer! It helped me to understand this better

Sign in to comment.

Categories

Find more on Linear Algebra in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!