# Check if a matrix is Positive Semidefinite

I have a question concerning the check whether a given matrix is positive semidefinite or not. In mathematica the function PositiveDefiniteMatrixQ[m] tells me whether the matrix m is positive, but not semidefinite. I wanted to ask whether there exists an implementation for that purpose, and if not, how to implement a code that is as fast as possible.

-
Is your matrix symbolic or numeric? –  The Toad Nov 19 '13 at 19:29
It is a 15x15 matrix simply consisting of integers. –  Dan Nov 19 '13 at 19:31
You could compute the eigenvalues numerically to do the check---that is very fast. –  Szabolcs Nov 19 '13 at 19:41
Given its all integers, you can check exact roots of the characteristic polynomial. –  Daniel Lichtblau Nov 19 '13 at 20:43

Computing the eigenvalues is actually overkill. It could be done by row reduction. For convenience, I set empty list to be positive semidefinite.

psdcheck[m_]:=(
If[Length@m== 0,Return@True];
If[Length@m== 1,Return@If[Negative@m[[1,1]],False,True]];
If[Or@@(Negative@Diagonal@m),Return@False];

Module[{mtemp,ind1},
mtemp=If[SymmetricMatrixQ@m,m,(#+Transpose@#)&@m];
ind1=Flatten@Position[Diagonal@mtemp,diag_/;diag== 0.];

If[Length@ind1!=0,
If[(#!=ConstantArray[0.,Dimensions@#])&@(mtemp[[ind1]]),Return@False,
mtemp=mtemp[[#,#]]&@Complement[Range@Length@mtemp,ind1];
If[Length@mtemp<= 1,Return@psdcheck@mtemp]];];

psdcheck[mtemp[[2;;,2;;]]-
(Outer[Times,#,#/mtemp[[1,1]]])&@mtemp[[1,2;;]]]
])

-
Is diag_/;diag== 0. the same as 0.? –  ybeltukov Jan 15 '14 at 9:49
Hi Hsien-Ching, good to see you are here too! –  Sjoerd C. de Vries Jan 15 '14 at 11:30
diag_/;diag==0 and 0. are different. Position uses pattern matching to identify the elements. 0. == 0 is True but 0. and 0 are different patterns. Good to see you here too, Sjoerd. –  Hsien-Ching Kao Jan 15 '14 at 15:37
There are numerical issues with the code during the pivoting step. This could be avoided by choosing the largest diagonal element as pivot instead of blindly choosing the first one. Besides, the code only works on real case but it is easy to generalize it to the complex situation. One could further utilize the symmetric property (after symmetrization) to reduce the memory usage and the number of multiplications. This code is significantly faster than the method based on eigenvalue calculation. –  Hsien-Ching Kao Jan 16 '14 at 5:11

semiPositiveDefiniteMatrixQ[m_, tol_: 10.^-10] :=
PositiveDefiniteMatrixQ[m + tol IdentityMatrix@Dimensions[m]]

PositiveDefiniteMatrixQ@{{1, 0}, {0, 0}}
semiPositiveDefiniteMatrixQ@{{1, 0}, {0, 0}}


False

True

Explanation:

Let the matrix $M$ has $\lambda_1, \lambda_2, \dots \lambda_n$ eigenvalues. Then the matrix $M+tI$ ($I$ is the identity matrix) has $\lambda_1+t, \lambda_2+t, \dots \lambda_n+t$ eigenvalues. If $M+tI$ is positive definite and $t>0$ is small enough (depends on the working precision) then $M$ is positive semidefinite.

-
This should work as long as the entries in the matrix are all small. Otherwise... maybe not. –  Oleksandr R. Jan 15 '14 at 18:10
@OleksandrR. I think it should works for any matrices. See my update, do you agree with it? –  ybeltukov Jan 16 '14 at 0:19
Sorry for not being clear--mathematically I absolutely agree, but I was thinking about the case when the matrix entries are larger in absolute magnitude than tol/\$MachineEpsilon (= 450 000 for the default tolerance). Then you effectively add the zero matrix (or at least not a full rank identity matrix), so the result will not (or may not) change. I am not sure how this situation can be rescued without resorting to calculating the eigenvalues. –  Oleksandr R. Jan 16 '14 at 4:06

If your matrix is numeric, the following should be a reasonably fast way to determine positive semi-definiteness:

Clear@positiveSemiDefiniteQ
positiveSemiDefiniteQ[mat_?MatrixQ] :=
With[{eigs = Chop@Eigenvalues@N@mat},
Quiet[
Check[
False,
GreaterEqual::nord
],
GreaterEqual::nord
]
]

-
Maybe it's a good idea to check if the matrix is symmetric (Hermitian). Usually positive definiteness is only defined for such matrices. –  Szabolcs Nov 19 '13 at 20:06
@Szabolcs Nope. Example: {{1, 0}, {1, 1}} –  The Toad Nov 19 '13 at 20:07

Here's a function that checks to see if the roots of the characteristic polynomial are all non-negative (as suggested by Daniel Lichtblau:

semiDefQ[m_] := Module[{p, roots},
p[x_] := CharacteristicPolynomial[m, x];
roots = Chop[Roots[p[z] == 0, z] // N];
Length[Select[Apply[List, roots[[All, 2]]], # ∈ Reals && # >= 0 &]] == n]


To use, just input your favorite matrix

n = 3;
m = RandomInteger[{0, 10}, {n, n}];
semiDefQ[m]

-