Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
Is your matrix symbolic or numeric? –  rm -rf 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
1  
Given its all integers, you can check exact roots of the characteristic polynomial. –  Daniel Lichtblau Nov 19 '13 at 20:43
add comment

4 Answers

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;;]]]
])
share|improve this answer
    
Is diag_/;diag== 0. the same as 0.? –  ybeltukov Jan 15 at 9:49
    
Hi Hsien-Ching, good to see you are here too! –  Sjoerd C. de Vries Jan 15 at 11:30
1  
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 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 at 5:11
add comment

One can add a small addition:

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.

share|improve this answer
    
This should work as long as the entries in the matrix are all small. Otherwise... maybe not. –  Oleksandr R. Jan 15 at 18:10
    
@OleksandrR. I think it should works for any matrices. See my update, do you agree with it? –  ybeltukov Jan 16 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 at 4:06
add comment

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[
                And @@ Thread[eigs >= 0],
                False,
                GreaterEqual::nord
            ],
            GreaterEqual::nord
        ]
    ]
share|improve this answer
1  
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}} –  rm -rf Nov 19 '13 at 20:07
add comment

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]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.