Nothing
# compute A^n where A is a square matrix, allowing non-integer and negative powers
# For integer values, see the technique in...
#https://en.wikipedia.org/wiki/Exponentiation_by_squaring
#' General Matrix Power
#'
#' Calculates the \code{n}-th power of a square matrix, where \code{n} can be a
#' positive or negative integer or a fractional power.
#'
#' @details
#' If \code{n<0}, the method is applied to \eqn{A^{-1}}.
#' When \code{n} is an
#' integer, the function uses the Russian peasant method, or repeated squaring
#' for efficiency.
#' Otherwise, it uses the spectral decomposition of \code{A},
#' \eqn{\mathbf{A}^n = \mathbf{V} \mathbf{D}^n \mathbf{V}^{T}}
#' requiring a symmetric matrix.
#'
#' @aliases mpower %^%
#' @param A A square matrix. Must also be symmetric for non-integer powers.
#' @param n matrix power
#' @return Returns the matrix \eqn{A^n}
#'
#' @author Michael Friendly
#' @references \url{https://en.wikipedia.org/wiki/Exponentiation_by_squaring}
#' @seealso Packages \pkg{corpcor} and \pkg{expm} define similar functions.
#' @keywords array
#' @examples
#'
#' M <- matrix(sample(1:9), 3,3)
#' mpower(M,2)
#' mpower(M,4)
#'
#' # make a symmetric matrix
#' MM <- crossprod(M)
#' mpower(MM, -1)
#' Mhalf <- mpower(MM, 1/2)
#' all.equal(MM, Mhalf %*% Mhalf)
#'
#'
#' @export
mpower <- function(A,n){
is.wholenumber <-
function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
if(!ncol(A)==nrow(A)) stop("Input must be a square matrix")
p <- nrow(A)
if (n==0) return (diag(p))
if (n==1) return (A)
if (n < 0 ) {
A <- solve(A)
n <- abs(n)
}
if (is.wholenumber(n)) {
result <- diag(p)
while (n > 0) {
if (n %% 2 != 0) {
result <- result %*% A
n <- n - 1
}
A <- A %*% A
n <- n / 2
}
}
else {
if( any( abs(A-t(A)) > 100*.Machine$double.eps ) )
stop("Input must be a symmetric matrix for non-integer powers")
result <- with(eigen(A), vectors %*% (values^n * t(vectors)))
}
dimnames(result) <- dimnames(A)
return(result)
}
#' Shorthand operator for mpower
#'
#' @rdname mpower
#' @export
"%^%" <- function(A,n) mpower(A,n)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.