Why is numpy saying that my totally positive matrix is not positive semi-definite?

Question

I generate a correlation matrix by drawing from a uniform distribution:

corr = np.random.uniform(low=0.1, high=0.3, size=[20, 20])

and set the main diagonal elements to one

corr[np.diag_indices_from(corr)] = 1

Then, I make the correlation matrix symmetric

corr[np.triu_indices(n=20, k=1)] = corr.T[np.triu_indices(n=20, k=1)]

which yields a totally positive matrix ,i.e., all values of the matrix are strictly positive.

According to numpy, however, the matrix is not positive (semi-) definit.

np.all(np.linalg.eigvals(corr) >= 0)
False

You're not setting the main diagonal elements to zero, but to one — Riccardo Bucco, Sep 20 '22 at 17:14
Because positive definite matrices are not defined that way. — Mechanic Pig, Sep 20 '22 at 17:14
With the updated code (set diagonal to 1), I can't reproduce this. I've run the code several times, and the result of `np.all(np.linalg.eigvals(corr) >= 0)` has been `True` each time. — Warren Weckesser, Sep 20 '22 at 18:26

Bob · Accepted Answer · 2022-09-20T18:47:06.067

That's still not guaranteed to be PSD

I will give you two easy ways:

Sny square non-singular matrix can be used to create a PSD matrix with

A = A @ A.T

Any matrix can be used to produce a PSD matrix with

A = (A + A.T)/2
A = A - np.eye(len(A)) * (np.min(np.linalg.eigvalsh(A)) - 0.001)

If you want the minimum perturbation to a symmetric matrix (the least squares projection to the positive semidefinite cone)

A_ = (v * np.maximum(w, 0.01)) @ v.T
print(np.linalg.eigvalsh(A_))

Notice that I am giving a margin of 0.01, if I used strictly zero your test could fail due to numerical errors.

1 Answers1