# sqrtm vs chol to produce draws the multivariate normal distribution

9 visualizaciones (últimos 30 días)
zym el 3 de Mzo. de 2023
Comentada: zym el 25 de Mzo. de 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?
##### 1 comentarioMostrar NingunoOcultar Ninguno
zym el 3 de Mzo. de 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.

Iniciar sesión para comentar.

John D'Errico el 25 de Mzo. de 2023
Editada: John D'Errico el 25 de Mzo. de 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 comentarioMostrar NingunoOcultar Ninguno
zym el 25 de Mzo. de 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

Iniciar sesión para comentar.

### Más respuestas (1)

the cyclist el 25 de Mzo. de 2023
Editada: the cyclist el 25 de Mzo. de 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 comentariosMostrar 1 comentario más antiguoOcultar 1 comentario más antiguo
zym el 25 de Mzo. de 2023
Thanks a lot for the answer! It helped me to understand this better

Iniciar sesión para comentar.

### Categorías

Más información sobre Discriminant Analysis en Help Center y File Exchange.

R2022b

### Community Treasure Hunt

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

Start Hunting!