IsPSD

From QETLAB
Revision as of 18:25, 16 December 2014 by Nathaniel (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
IsPSD
Determines whether or not a matrix is positive semidefinite

Other toolboxes required none
Related functions IsCP
IsPPT
IsTotallyPositive
Function category Basic operation

IsPSD is a function that determines whether or not a given matrix is positive semidefinite. The input matrix can be either full or sparse and, if requested, a vector that proves that the given matrix is not positive semidefinite can be provided as output.

Syntax

  • PSD = IsPSD(X)
  • PSD = IsPSD(X,TOL)
  • [PSD,WIT] = IsPSD(X,TOL)

Argument descriptions

Input arguments

  • X: A square matrix.
  • TOL (optional, default eps^(3/4)): The numerical tolerance used when determining positive semidefiniteness. The matrix will be determined to be positive semidefinite if its minimal eigenvalue is computed to be at least -TOL.

Output arguments

  • PSD: A flag (either 1 or 0) indicating that X is or is not positive semidefinite.
  • WIT (optional): An eigenvector corresponding to the minimal eigenvalue of X. When PSD = 0, this serves as a witness that verifies that X is not positive semidefinite, since WIT'*X*WIT < 0.

Examples

Simple example with low tolerance

When X is very simple, positive semidefiniteness can be be determined exactly. The following example has the TOL = 0 (not recommended in general!) to highlight the fact that the script really is checking for positive semidefiniteness, not positive definiteness.

>> X = diag([1 0])
 
X =
 
     1     0
     0     0
 
>> IsPSD(X,0)
 
ans =
 
     1

Furthermore, if we make one of the eigenvalues even slightly negative in this case, it is detected as not positive semidefinite:

>> IsPSD(X-eps*eye(2),0)
 
ans =
 
     0

Note that in general you can not expect this kind of accuracy.

Notes

Do not request the WIT output argument unless you need it. If WIT is not requested, positive semidefiniteness is determined by attempting a Cholesky decomposition of X, which is both faster and more accurate than computing its minimum eigenvalue/eigenvector pair.

Source code

Click on "expand" to the right to view the MATLAB source code for this function.

  1. %%  ISPSD    Determines whether or not a matrix is positive semidefinite
  2. %   This function has one required argument:
  3. %     X: a square matrix
  4. %
  5. %   PSD = IsPSD(X) is either 1 or 0, indicating that X is or is not
  6. %   positive semidefinite (within reasonable numerical error).
  7. %
  8. %   This function has one optional input argument:
  9. %     TOL (default eps^(3/4))
  10. %
  11. %   [PSD,WIT] = IsPSD(X,TOL) determines whether or not X is positive
  12. %   semidefinite within the tolerance specified by TOL. WIT is the
  13. %   eigenvector corresponding to the minimal eigenvalue of X, and thus can
  14. %   act as a witness that proves X is not positive semidefinite (i.e.,
  15. %   WIT'*X*WIT < 0).
  16. %
  17. %   URL: http://www.qetlab.com/IsPSD
  18.  
  19. %   requires: opt_args.m
  20. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  21. %   package: QETLAB
  22. %   last updated: December 24, 2014
  23.  
  24. function [psd,wit] = IsPSD(X,varargin)
  25.  
  26. if(size(X,1) ~= size(X,2))
  27.     psd = 0;
  28.     wit = 0;
  29.     return
  30. end
  31.  
  32. % set optional argument defaults: tol=eps^(3/4)
  33. [tol] = opt_args({ eps^(3/4) },varargin{:});
  34.  
  35. % Allow this function to be called within CVX optimization problems.
  36. if(isa(X,'cvx'))
  37.     cvx_begin sdp quiet
  38.     subject to
  39.     	X >= 0;
  40.     cvx_end
  41.     psd = 1-min(cvx_optval,1); % CVX-safe way to map (0,Inf) to (1,0)
  42.  
  43. % If the function is just being called on a non-CVX variable, just check
  44. % positive semidefiniteness normally (which is much faster).
  45. else
  46.     % only check the Hermitian part of X
  47.     X = (X+X')/2;
  48.  
  49.     % if the user requested the smallest eigenvector, compute it
  50.     if(nargout > 1)
  51.         if(isreal(X))
  52.             eigs_id = 'SA';
  53.         else
  54.             eigs_id = 'SR';
  55.         end
  56.         [wit,eigval] = eigs(X,1,eigs_id);
  57.         psd = (eigval >= -tol);
  58.  
  59.     % otherwise, use chol to determine positive semidefiniteness, which is both faster and more accurate
  60.     else
  61.         [~,p] = chol(X+(tol+eps)*speye(length(X)));
  62.         psd = (p == 0);
  63.     end
  64. end