Home > Software engineering >  Why is numpy saying that my totally positive matrix is not positive semi-definite?
Why is numpy saying that my totally positive matrix is not positive semi-definite?

Time:09-21

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

CodePudding user response:

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.

  • Related