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?

Question

奥刘 2023-10-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2040101-how-to-deal-with-the-hermitian-matrix-that-is-meant-to-be-positive-semidefinite-but-turned-out-to-ha

编辑：奥刘 2023-10-30

采纳的回答： John D'Errico

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 corresponding code is at https://github.com/la0719/temp/blob/main/main

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.

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.

Secondly, how could I resonably resolve this problem to achieve a cov without negative eigenvalue? Maybe somehow set a reasonable precision?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2023-10-29

编辑：John D'Errico 2023-10-29

在 MATLAB Online 中打开

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.)

https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd?s_tid=srchtitle

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 个评论
显示 1更早的评论隐藏 1更早的评论

John D'Errico 2023-10-29

编辑：John D'Errico 2023-10-29

在 MATLAB Online 中打开

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.

奥刘 2023-10-30