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?
15 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
John D'Errico
2023-10-29
编辑:John D'Errico
2023-10-29
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 个评论
John D'Errico
2023-10-29
编辑:John D'Errico
2023-10-29
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.
更多回答(1 个)
Matt J
2023-10-29
You can add a small multiple of eye(N) to raise the eigenvalues above zero
F=F+small_number*eye(size(F))
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!