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?

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?

Answer 1

[John D'Errico]

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

[奥 刘]

Really thanks for the function and detailed explanation here, John! I will try this method. Btw, is there any method in your mind maybe, that might use only the bulit-in functions of MATLAB to cope with this round off issue?

[John D'Errico]

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.

[奥 刘]

Thanks again! I think now I totaly get it here.

Answer 2

[Matt J]

You can add a small multiple of eye(N) to raise the eigenvalues above zero

F=F+small_number*eye(size(F))

[奥 刘]

Really thanks for the answer. Yes, but then we have to decide the diagnoal value we added here, like how large should it be.