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