How to deal with the Hermitian matrix that is meant to be positive semidefinite but turned out to have negative eigenvalue due to numerical issues?
14 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
奥 刘
el 29 de Oct. de 2023
So I am using CVX to solve a Boolean least square problem, the code is shown below
cvx_begin
variable F(J,J) hermitian semidefinite
variable v(J,1) complex
minimize( real(trace(Z * F)) - 2 * real(J * as_ele_Tx.' * v) + J^2 )
subject to
for i = 1:(K_num)
trace(Hm(i,:)' * Hm(i,:) * F) <= gamma_seq(i2) ./ P;
end
diag(F) == 1;
[F v;v' 1] == hermitian_semidefinite(J+1);
cvx_end
% extract the weight solution by randomization
mu = v; % expectation
cov = F - v * v'; % covariance
L = 5000; % number of samples
% randomization procedure
z = 1 ./ sqrt(2) .* (randn(J, L) + 1j .* randn(J, L));
f_L = mu .* ones(J, L) + chol(cov)' * z ;
The problem can be well solved by CVX with solution of
and
. But, when I try to use the Gaussian randomization to extract a vector solution from the positive semidefinite matrix
and vector
, although I have constraint
, the matrix cov = F - v * v' will have a trivial negative eigenvalue, which leads to the error of chol(cov) operation, where for Cholesky factorization requires cov to be positive definite.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1524026/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1524031/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1524026/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1524031/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1524036/image.png)
So, my first problem is that why this negative eigenvalue happen? As the CVX solved this problem with the
constraint, it doesn't seem to be a CVX problem to me but because of the numerical issues here.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1524036/image.png)
Secondly, how could I resonably resolve this problem to achieve a cov without negative eigenvalue? Maybe somehow set a reasonable precision?
0 comentarios
Respuesta aceptada
John D'Errico
el 29 de Oct. de 2023
Editada: John D'Errico
el 29 de Oct. de 2023
This is not a question of a precision you can set.
You can download my nearestSPD function from the file exchange. It finds the smallest perturbation to a matrix such that the matrix will be SPD, where chol will now succeed. (That is the common test in MATLAB for a matrix to be positive definite.)
The code is based on a Nick Higham algorithm.
For example:
A = rand(4,3);
A = A*A.';
chol(A)
Error using chol
Matrix must be positive definite.
So a fails the chol test.
A
A =
1.524 0.87422 1.6569 1.1529
0.87422 1.0754 1.1442 0.98086
1.6569 1.1442 1.9891 1.6667
1.1529 0.98086 1.6667 1.8131
In fact, eig even thinks it may be just barely positive definite, but chol fails.
eig(A)
ans =
3.4315e-16
0.38764
0.53059
5.4834
However, nearestSPD tweaks the matrix just slightly, enough for chol to now be happy.
chol(nearestSPD(A))
ans =
1.2345 0.70815 1.3422 0.93392
0 0.7576 0.25576 0.42174
0 0 0.34964 0.87352
0 0 0 2.5564e-08
The difference is tiny of course, down in the least significant bits of A.
A - nearestSPD(A)
ans =
0 0 0 0
0 2.2204e-16 -2.2204e-16 -1.1102e-16
0 -2.2204e-16 -2.2204e-16 0
0 -1.1102e-16 0 0
3 comentarios
John D'Errico
el 29 de Oct. de 2023
Editada: John D'Errico
el 29 de Oct. de 2023
Can you use only MATLAB functions? nearestSPD only uses MATLAB functions! It is written in MATLAB. So I don't really see the point of your question. They are just collected into one function call. It is designed to be a minimal perturbation so the result will be SPD, and where you don't need to tweak anything.
Can you add a multiple of the identity matrix to your matrix? Well, yes. If the matrix is already symmetric. But even then, you need to guess how large of a multiple to add! Just using eig may not be sufficient here.
For example, in the case I used before, we had:
eig(A)
ans =
3.4315e-16
0.38764
0.53059
5.4834
So as far as eig is concerned, A already is SPD, all positive eigenvalues. But due to roundoff errors, chol still failed.
fac = 2;
chol(A + eye(size(A))*abs(fac*min(eig(A))))
ans =
1.2345 0.70815 1.3422 0.93392
0 0.7576 0.25576 0.42174
0 0 0.34964 0.87352
0 0 0 9.8292e-08
Something like that MAY work, as long as A is already symmetric. But the choice of fac==2 is arbitrary, and it may also easily fail on some other matrices. It may need to be larger or smaller, and that forces you to fool around to find the right multiple, as small as possible, yet just large enough.
Más respuestas (1)
Matt J
el 29 de Oct. de 2023
You can add a small multiple of eye(N) to raise the eigenvalues above zero
F=F+small_number*eye(size(F))
1 comentario
Ver también
Categorías
Más información sobre Logical en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!