Effacer les filtres
Effacer les filtres

sqrtm vs chol to produce draws the multivariate normal distribution

25 vues (au cours des 30 derniers jours)
zym
zym le 3 Mar 2023
Commenté : zym le 25 Mar 2023
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 commentaire
zym
zym le 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.

Connectez-vous pour commenter.

Réponse acceptée

John D'Errico
John D'Errico le 25 Mar 2023
Modifié(e) : John D'Errico le 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 commentaire
zym
zym le 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

Connectez-vous pour commenter.

Plus de réponses (1)

the cyclist
the cyclist le 25 Mar 2023
Modifié(e) : the cyclist le 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 commentaires
the cyclist
the cyclist le 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 le 25 Mar 2023
Thanks a lot for the answer! It helped me to understand this better

Connectez-vous pour commenter.

Catégories

En savoir plus sur Linear Algebra dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by