IsTotallyPositive

From QETLAB
Jump to: navigation, search
IsTotallyPositive
Determines whether or not a matrix is totally positive

Other toolboxes required none
Related functions IsTotallyNonsingular
Function category Miscellaneous

IsTotallyPositive is a function that determines whether or not a given matrix is totally positive (i.e., all of its square submatrices have positive determinant). The input matrix can be either full or sparse.

Syntax

  • ITP = IsTotallyPositive(X)
  • ITP = IsTotallyPositive(X,SUB_SIZES)
  • ITP = IsTotallyPositive(X,SUB_SIZES,TOL)
  • [ITP,WIT] = IsTotallyPositive(X,SUB_SIZES,TOL)

Argument descriptions

Input arguments

  • X: A matrix.
  • SUB_SIZES (optional, default 1:min(size(X))): A vector specifying the sizes of submatrices to be checked to have positive determinant.
  • TOL (optional, default length(X)*eps(norm(X,'fro'))): The numerical tolerance used when determining positivity.

Output arguments

  • ITP: A flag (either 1 or 0) indicating that X is or is not totally positive.
  • WIT (optional): If ITP = 0 then WIT specifies a submatrix of X that has either negative or non-real determinant. More specifically, WIT is a matrix with 2 rows such that det(X(WIT(1,:),WIT(2,:))) is negative or non-real.

Examples

Vandermonde matrices

It is known that Vandermonde matrices are totally positive, under certain restrictions on the nodes:

>> IsTotallyPositive(vander(5:-1:1))
 
ans =
 
     1

Notes

In practice, this function is practical for matrices of size up to about 15-by-15.

Source code

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

  1. %%  ISTOTALLYPOSITIVE    Determines whether or not a matrix is totally positive
  2. %   This function has one required argument:
  3. %     X: a matrix
  4. %
  5. %   ITP = IsTotallyPositive(X) is either 1 or 0, indicating that X is or
  6. %   is not totally positive (i.e., all of its square submatrices have
  7. %   positive real determinant).
  8. %
  9. %   This function has two optional input arguments:
  10. %     SUB_SIZES (default 1:min(size(X)))
  11. %     TOL (default max(size(X))*eps(norm(X,'fro')))
  12. %
  13. %   [ITP,WIT] = IsTotallyPositive(X,SUB_SIZES,TOL) determines whether or
  14. %   not every r-by-r submatrix of X has positive determinant, where r
  15. %   ranges over all values in the vector SUB_SIZES, and positivity is
  16. %   determined within a tolerance of TOL. If ITP = 0 then WIT is a matrix
  17. %   with two rows and r columns that specifies an r-by-r submatrix of X
  18. %   that does not have positive determinant.
  19. %
  20. %   URL: http://www.qetlab.com/IsTotallyPositive
  21.  
  22. %   requires: opt_args.m
  23. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  24. %   package: QETLAB
  25. %   last updated: December 13, 2012
  26.  
  27. function [itp,wit] = IsTotallyPositive(X,varargin)
  28.  
  29. sX = size(X);
  30. wit = 0;
  31.  
  32. % set optional argument defaults: sub_sizes=1:min(size(X)), tol=max(size(X))*eps(norm(X,'fro')) (same as rank)
  33. [sub_sizes,tol] = opt_args({ 1:min(sX), max(sX)*eps(norm(X,'fro')) },varargin{:});
  34.  
  35. for j = 1:length(sub_sizes)
  36.     % 1x1 determinants can be computed quickly, so do them separately
  37.     if(sub_sizes(j) == 1)
  38.         [r,c] = find(min(real(X),-abs(imag(X))) < -tol,1);
  39.         if(min(size(r)) > 0)
  40.             itp = 0;
  41.             wit = [r;c];
  42.             return;
  43.         end
  44.  
  45.     % larger determinants are slower; just loop on through!
  46.     else
  47.         sub_ind_r = nchoosek(1:sX(1),sub_sizes(j));
  48.         if(sX(1) == sX(2)) % nchoosek is slightly slow, so only call it once if we can get away with it
  49.             sub_ind_c = sub_ind_r;
  50.         else
  51.             sub_ind_c = nchoosek(1:sX(2),sub_sizes(j));
  52.         end
  53.         sub_ind_len_r = size(sub_ind_r,1);
  54.         sub_ind_len_c = size(sub_ind_c,1);
  55.         for kr = 1:sub_ind_len_r
  56.             for kc = 1:sub_ind_len_c
  57.                 d = det(X(sub_ind_r(kr,:),sub_ind_c(kc,:)));
  58.                 if(d < -tol || abs(imag(d)) > tol)
  59.                     itp = 0;
  60.                     wit = [sub_ind_r(kr,:);sub_ind_c(kc,:)];
  61.                     return;
  62.                 end
  63.             end
  64.         end
  65.     end
  66. end
  67.  
  68. itp = 1;