Nothing
#' @title Matrix Logarithm
#' @description \code{matlog} computes the (principal) matrix logarithm of a square matrix.
#' @details This function adapts function \code{logm} from package \pkg{expm}
#' to be able to handle complex matrices, by decomposing the original matrix
#' into real and purely imaginary matrices and creating a real
#' block matrix that function \code{logm} can successfully process.
#' If the original matrix is real, \code{matlog} calls \code{logm} directly for maximum efficiency.
#' Hence, for real matrices, \code{matlog} can compute the matrix logarithm in the same instances as \code{logm};
#' for complex matrices, \code{matlog} can compute the matrix logarithm as long as all real
#' eigenvalues are positive: zero eigenvalues imply singularity (and therefore the log
#' does not exist) and negative eigenvalues can be problematic as it may be hard and
#' numerically unstable to calculate Jordan blocks. See references below.
#' @references For more on the matrix logarithm, visit
#' \url{https://en.wikipedia.org/wiki/Logarithm_of_a_matrix}
#' @param A a square matrix, real or complex.
#' @param ... arguments passed to or from other methods.
#' @return The matrix logarithm of A. Method used may be chosen from the options available in \code{logm}.
#' @author Albert Dorador
#' @export
#' @import expm
#' @seealso \code{\link[expm]{logm}}
#' @examples
#' A <- matrix(c(1, 2, 2+3i, 5), ncol = 2) # complex matrix
#' B <- matrix(c(2, 0, 3, 4), ncol = 2) # real matrix with positive eigenvalues
#'
#' matlog(A)
#' matlog(A, "Eigen") # uses Eigen method in function logm
#' matlog(B)
#' matlog(B, "Eigen") # uses Eigen method in function logm
#'
matlog <- function(A, ...)
{
d <- dim(A)
n <- d[1]
p <- d[2]
if (n != p)
stop("Supplied matrix is not square")
if (is.complex(A)){
e.val <- eigen(A, only.values = TRUE)$values
safe <- TRUE
i <- 1
while (safe && i <= n){
if(class(Imzap(e.val[i])) != "complex"){ #i.e. if ith eigenvalue is real
if(Imzap(e.val[i]) <= 0){ #0 would make the matrix singular, and negative is likely to cause trouble
safe <- FALSE
}
}
i <- i + 1
}
if (safe){ #i.e. only if taking the log is safe
# Extract real and imaginary parts
Ar <- Re(A)
Ai <- Im(A)
# Construct real but extended matrix
L <- rbind( cbind(Ar,-Ai), cbind(Ai,Ar))
# Compute exponential of that
lL <- logm(L, ...)
# Construct logm(A) by extracting real and imaginary parts from blocks
lA <- lL[1:n,1:n] + 1i*lL[(1:n)+n,1:n]
lA <- matrix(Imzap(lA), ncol = n) # to guarantee best-looking solutions
return(lA)
}else{
stop("The matrix logarithm may not exist for this matrix.")
}
}else{
logm(A, ...)
}
}
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.