1isSemidefinite <- function( m, ... )
2   UseMethod( "isSemidefinite" )
3
4isSemidefinite.default <- function( m, ... ) {
5   stop( "there is currently no default method available" )
6}
7
8## ----- test positive / negative semidefiniteness
9isSemidefinite.matrix <- function( m, positive = TRUE,
10   tol = 100 * .Machine$double.eps,
11   method = ifelse( nrow( m ) < 13, "det", "eigen" ), ... ) {
12
13   if( !is.matrix( m ) ) {
14      stop( "argument 'm' must be a matrix" )
15   } else {
16      if( nrow( m ) != ncol( m ) ) {
17         stop( "argument 'm' or each of its elements must be a _quadratic_ matrix" )
18      } else if( !isSymmetric( unname( m ), tol = 1000 * tol ) ) {
19         stop( "argument 'm' must be a symmetric matrix" )
20      }
21      # make sure that the matrix is almost exactly symmetric
22      # even if it is slightly non-symmetric
23      m <- ( m + t(m) ) / 2
24      n <- nrow( m )
25      if( !positive ) {
26         m <- -m
27      }
28      if( n >= 12 && method == "det" ) {
29         warning( "using method 'det' can take a very long time",
30            " for matrices with more than 12 rows and columns;",
31            " it is suggested to use method 'eigen' for larger matrices",
32            immediate. = TRUE )
33      }
34      if( method == "det" ) {
35         for( i in 1:n ) {
36            comb <- combn( n, i )
37            for( j in 1:ncol( comb ) ) {
38               mat <- m[ comb[ , j ], comb[ , j ], drop = FALSE ]
39               if( rcond( mat ) >= tol ) {
40                  princMin <-  det( mat )
41               } else {
42                  princMin <- 0
43               }
44               if( princMin < -tol ) {
45                  return( FALSE )
46               }
47            }
48         }
49         return( TRUE )
50      } else if( method == "eigen" ) {
51         ev <- eigen( m, only.values = TRUE )$values
52         if( is.complex( ev ) ) {
53            stop( "complex (non-real) eigenvalues,",
54               " which could be caused by a non-symmetric matrix" )
55         }
56         if( all( ev > -tol ) ) {
57            return( TRUE )
58         } else {
59            if( rcond( m ) >= tol || n == 1 ) {
60               return( FALSE )
61            } else {
62               k <- max( 1, min( sum( abs( ev ) <= tol ), n - 1 ) )
63               comb <- combn( n, n-k )
64               for( j in 1:ncol( comb ) ) {
65                  mm <- m[ comb[ , j ], comb[ , j ], drop = FALSE ]
66                  if( !semidefiniteness( mm, tol = tol, method = method  ) ) {
67                     return( FALSE )
68                  }
69               }
70               return( TRUE )
71            }
72         }
73      } else {
74         stop( "argument 'method' must be either 'det' or 'eigen'" )
75      }
76   }
77   stop( "internal error: please inform the maintainer",
78      " of the 'miscTools' package",
79      " (preferably with a reproducible example)" )
80}
81
82isSemidefinite.list <- function( m, ... ) {
83
84   if( !is.list( m ) ) {
85      stop( "argument 'm' must be a matrix or a list of matrices" )
86   } else if( !all( sapply( m, is.matrix ) ) ) {
87      stop( "all components of the list specified by argument 'm'",
88         " must be matrices" )
89   }
90
91   result <- logical( length( m ) )
92   for( t in 1:length( m ) ) {
93      result[ t ] <- isSemidefinite( m[[ t ]], ... )
94   }
95   return( result )
96}
97
98semidefiniteness <- function( m, ... ) {
99   result <- isSemidefinite( m = m, ... )
100   return( result )
101}
102