-
Notifications
You must be signed in to change notification settings - Fork 16
/
mpower.R
32 lines (30 loc) · 1.27 KB
/
mpower.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#' Matrix Power
#'
#' A simple function to demonstrate calculating the power of a square symmetric matrix in terms of its eigenvalues and eigenvectors.
#'
#' @details The matrix power \code{p} can be a fraction or other non-integer. For example, \code{p=1/2} and
#' \code{p=1/3} give a square-root and cube-root of the matrix.
#'
#' Negative powers are also allowed. For example, \code{p=-1} gives the inverse and \code{p=-1/2}
#' gives the inverse square-root.
#'
#' @param A a square symmetric matrix
#' @param p matrix power, not necessarily a positive integer
#' @param tol tolerance for determining if the matrix is symmetric
#' @return \code{A} raised to the power \code{p}: \code{A^p}
#' @seealso The \code{{\%^\%}} operator in the \pkg{expm} package is far more efficient
#' @export
#' @examples
#' C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
#' C
#' mpower(C, 2)
#' zapsmall(mpower(C, -1))
#' solve(C) # check
mpower <- function(A, p, tol=sqrt(.Machine$double.eps)) {
if (!is.numeric(A) || !is.matrix(A) || nrow(A) != ncol(A) || any(abs(A - t(A)) > tol))
stop("A must be a numeric, square, symmetric matrix")
e <- eigen(A)
M <- e$vectors # matrix for changing basis
d <- e$values # eigen values
return(M %*% diag(d^p) %*% t(M))
}