Commit 1cf822ab authored by Corson N. Areshenkoff's avatar Corson N. Areshenkoff

Initial commit

parent 1c22872d
Package: spdm
Type: Package
Title: Functions for working with symmetric positive-definite matrices
Version: 0.1.0
Author: Corson N. Areshenkoff
Maintainer: Corson N. Areshenkoff <[email protected]>
Description: Means, differencing, and centering for symmetric positive-definite matrices.
Depends: R (>= 3.0.2),
Rcpp (>= 0.12.0)
Imports: expm, PDSCE, huge
License: GPL-3
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.0.1
Suggests: knitr,
rmarkdown
VignetteBuilder: knitr
exportPattern("^[[:alpha:]]+")
importFrom(expm,expm)
importFrom(expm,logm)
importFrom(expm,sqrtm)
importFrom(huge,huge)
importFrom(huge,huge.select)
importFrom(PDSCE,pdsoft.cv)
#' Check if a matrix is symmetric and positive definite
#' @param ... One or more numeric matrices
#' @param tol Numerical tolerance for checking symmetry and positive definiteness
#' @details Function checks whether the matrix \code{x} is symmetric and positive
#' definite. Symmetry is evaluated up to an entrywise tolerance of \code{tol}, so
#' that differences smaller than \code{tol} are ignored. Positive-definiteness is
#' checked by computing the eigenvalues of \code{x} using \code{eigen}, setting
#' eigenvalues smaller than \code{tol} in absolute value to zero, and then checking
#' whether any are less than or equal to zero.
#'
#' @return A Boolean value indicating whether the matrix is symmetric and positive
#' definite
is.spd <- function(..., tol = 1e-8){
data <- list(...)
check <- sapply(data, function(i) {
# Test if matrix
if (!is.matrix(i)){
return(FALSE)
}
# Test if symmetric
if (!isSymmetric(i, tol = tol)){
return(FALSE)
}
# Test if positive definite
eig <- eigen(i, symmetric = T, only.values = T)$values
eig[abs(eig) <= tol] <- 0
if (any(eig <= 0)){
return(FALSE)
}
return(TRUE)
})
}
is.valid.input <- function(x, type){
if (type == 'spd.mat'){
if (!is.numeric(x) | !is.matrix(x)) return(FALSE)
if (any(Im(x) != 0)) return(FALSE)
if (!is.spd(x)) return(FALSE)
} else if (type == 's.mat') {
if (!is.numeric(x) | !is.matrix(x)) return(FALSE)
if (any(Im(x) != 0)) return(FALSE)
if (dim(x)[1] != dim(x)[2]) return(FALSE)
if (!isSymmetric(x)) return(FALSE)
} else if (type == 'spd.list'){
if (!is.list(x)) return(FALSE)
if (!length(x) > 1) return(FALSE)
if (any(sapply(x, is.null))) return(FALSE)
chk <- sapply(x, function(i) {
if (!is.numeric(i) | !is.matrix(i)) return(FALSE)
if (any(Im(i) != 0)) return(FALSE)
if (dim(i)[1] != dim(i)[2]) return(FALSE)
if (!is.spd(i)) return(FALSE)
return(TRUE)
})
if (any(!chk)) return(FALSE)
if (length(unique(sapply(x, function(i) dim(i)[1]))) > 1) return(FALSE)
} else if (type == 's.list') {
if (!is.list(x)) return(FALSE)
if (!length(x) > 1) return(FALSE)
if (any(sapply(x, is.null))) return(FALSE)
chk <- sapply(x, function(i) {
if (!is.numeric(i) | !is.matrix(i)) return(FALSE)
if (any(Im(i) != 0)) return(FALSE)
if (dim(i)[1] != dim(i)[2]) return(FALSE)
if (!isSymmetric(i)) return(FALSE)
return(TRUE)
})
if (any(!chk)) return(FALSE)
if (length(unique(sapply(x, function(i) dim(i)[1]))) > 1) return(FALSE)
}
return(TRUE)
}
input.type <- function(x){
if (!is.matrix(x) & !is.list(x)){
stop('Input must be either a matrix or a list of matrices')
}
# Check input type
type <- c('spd.mat', 's.mat', 'spd.list', 's.list')
chk <- sapply(type, function(i) {
is.valid.input(x, i)
})
if (!any(chk)){
stop('Invalid input type. Probably, input matrix fails to be symmetric')
} else {
return(type[which(chk)])
}
}
#' Center a set of spd matrices
#'
#' Functions centers a set of symmetric, positive-definite matrices so that they
#' have mean equal to the identity matrix.
#'
#' @param x Either a list of symmetric, positive-definite matrices, an individual amtrix
#' @param y Value(s) for centering. Either a list of spd matrices of the same
#' length as x, a single spd matrix, or NULL. If NULL, uses the mean of \code{x}
#' @details Functions transports the matrices in \code{x} so that \code{y} lies at
#' the identity. If \code{x} is a list and \code{y = NULL} (default), function sets
#' \code{y} equal to the mean of \code{x}. Otherwise, \code{y} can be set to any
#' spd matrix. When \code{x} is a matrix, \code{y} must be provided. If \code{y}
#' is a list of the same length as {x}, the function centers each element of
#' \code{x} by the corresponding element of \code{y}.
#'
#' Centering is accomplished by the parallel transport specified by XXX (2XXX), by
#' first computing the inverse square root of \code{y}, and then centering \code{x}
#' by \code{y.inv.sqrt x y.inv.sqrt}.
#'
#' @return An object of the same type as \code{x}, containing the centered matrices..
spd.center <- function(x, y = NULL){
if (is.list(x)){
# If y is a list, do elementwise
if (is.list(y)){
if (length(x) != length(y)){
stop('If a list, y must have the same length as x')
}
# inv-sqrm
y.inv.sqrt <- lapply(y, function(i) sqrtm(solve(i)))
# Center
x.centered <- lapply(1:length(x), function(i){
y.inv.sqrt[i] %*% x[i] %*% y.inv.sqrt[i]
})
}
# If no center specified, use mean
if (is.null(y)){
y <- spd.mean(x)
}
# inv-sqrm
y.inv.sqrt <- sqrtm(solve(y))
# Center
x.centered <- lapply(x, function(i)
y.inv.sqrt %*% i %*% y.inv.sqrt)
} else if (is.matrix(x)) {
# If no center specified, abort
if (is.null(y)){
stop('Must provide y if x is a single matrix')
}
# inv-sqrm
y.inv.sqrt <- sqrtm(solve(y))
# Center
x.centered <- y.inv.sqrt %*% x %*% y.inv.sqrt
}
return(x.centered)
}
#' Distance between two symmetric, positive-definite matrices
#'
#' Function implements several distance measures between symmetric,
#' positive-definite matrices, not all of which are proper metrics.
#'
#' @param x,y Symmetric, positive-definite matrices
#' @param method The distance measure. See details.
#' @details Allowable distance measures are
#' \itemize{
#' \item{"frobenius": }{The Frobenius norm of the difference \code{x-y}. Not affinely invariant.}
#' \item{"cholesky": }{The Frobenius norm of the difference between the cholesky factors
#' of \code{x} and \code{y}. Not affinely invariant.}
#' \item{"affine": }{An affinely invariant distance measure, computed by taking
#' the Frobenius norm of the projection of \code{y} into the tangent space at \code{x}}
#' \item{"riemannian": }{The Riemmanian distance proposed by Barachant, et al. (2013)}
#' }
spd.dist <- function(x, y, method = 'riemannian'){
if (!'spd.mat' %in% input.type(x) | !'spd.mat' %in% input.type(y)){
return('Both inputs must be positive definite matrices')
}
if (method == 'frobenius'){
return(norm(x - y, type = 'F'))
}
if (method == 'cholesky'){
return(norm(chol(x) - chol(y), type = 'F'))
}
if (method == 'affine'){
isr <- solve(sqrtm(x))
return(norm(logm(isr %*% y %*% isr), type = 'F'))
}
if (method == 'riemannian'){
return(norm(logm(solve(x) %*% y), type = 'F'))
}
}
#' Regularized covariance estimation
#'
#' Function performs regularized covariance estimation using either soft-thresholding
#' or the graphical lasso, involving l1 penalized estimation of the inverse
#' covariance (precision) matrix. In either case, the optimal tuning parameter
#' is selected by cross-validation. Currently, the function's internals are set
#' to my own default settings, for simplicity.
#'
#' @param y A data matrix, where rows are observations and columns are variables
#' @param method Either soft thresholding ('threshold') or graphical lasso ('sparse')
#' @return A covariance matrix
spd.estimate <- function(y, method = 'threshold'){
if (method == 'threshold'){
fit <- pdsoft.cv(y)$sigma
} else if (method == 'sparse') {
n <- dim(y)[1]
est <- huge(y, method = 'glasso', scr = T, scr.num = n/log(n),
cov.output = T)
sel <- huge.select(est, criterion = 'stars')
fit <- as.matrix(sel$opt.cov)
}
return(fit)
}
#' Exponential map onto the space of SPD matrices
#'
#' Function projects one or more vectors onto
#' the space of SPD matrices space at a point \code{p}
#'
#' @param x A symmetric matrix
#' @param p An spd matrix from whose tangent space to perform the projection. If
#' unspecified, the identity matrix is chosen.
#' @details Function uses the exponential map to project a symmetric
#' matrix onto the space of SPD matrices at a matrix \code{p}.
spd.expmap <- function(x, p = NULL){
# Check x input
if (!'s.mat' %in% input.type(x)){
stop('x must be a symmetric matrix')
}
# Set and check p
if (is.null(p)){
p <- diag(rep(1, dim(x)[1]))
} else if (!'spd.mat' %in% input.type(p)){
stop('p must be a symmetric positive definite matrix')
}
p.sqrt <- sqrtm(p)
p.inv.sqrt <- solve(p.sqrt)
return(p.sqrt %*% expm(p.inv.sqrt %*% x %*% p.inv.sqrt) %*% p.sqrt)
}
#' Fractional anisotropy of a symmetrix, positive-definite matrix
#'
#' @param x A symmetric, positive-definite matrix
spd.fa <- function(x){
# Check x input
if (!'s.mat' %in% input.type(x)){
stop('x must be a symmetric matrix')
}
k <- dim(x)[1]
l <- eigen(x, symmetric = T, only.values = T)$values
l.bar <- mean(l)
fa <- sqrt( (k)/(k-1) * sum( (l - l.bar)^2) / sum(l^2) )
return(fa)
}
#' Log-map onto the tangent space
#'
#' Function projects one or more symmetric, positive-definite matrices onto
#' the tangent space at a point \code{p}
#'
#' @param x An spd matrix
#' @param p An spd matrix on whose tangent space to perform the projection. If
#' unspecified, the identity matrix is chosen.
#' @details Function uses the log-map to project a symmetric, positive-definite
#' matrix onto the tangent space of the positive cone at a matrix \code{p}.
spd.logmap <- function(x, p = NULL){
# Check x input
if (!'spd.mat' %in% input.type(x)){
stop('x must be a symmetric positive definite matrix')
}
# Set and check p
if (is.null(p)){
p <- diag(rep(1, dim(x)[1]))
} else if (!'s.mat' %in% input.type(p)){
stop('p must be a symmetric positive definite matrix')
}
p.sqrt <- sqrtm(p)
p.inv.sqrt <- solve(p.sqrt)
return( p.sqrt %*% logm(p.inv.sqrt %*% x %*% p.inv.sqrt) %*% p.sqrt )
}
#' Compute the mean of a set of spd matrices
#'
#' Function computes the mean of a set of symmetric positive-definite matrices.
#' Several methods are implemented, as described in \code{Details}.
#'
#' @param x A list of symmetric, positive-definite matrices
#' @param method The type of mean to compute. See details
#' @param ... Further arguments. See details
#' @details Function computes the mean of a set of symmetrix, positive-definite
#' matrices. Several methods are implemented:
#' \itemize{
#' \item{"euclidean": }{The ordinary arithmetic mean -- the sum of the matrices in x,
#' divided by the number of matrices. This is guaranteed to be an spd matrix, but
#' does not necessarily preserve the spectral characteristics of the individual matrices.}
#' \item{"logeuclidean": }{Computed by taking the arithmetic mean of the logarithms
#' of the matrices in \code{x}, and then projecting back onto the space of spd matrices.
#' In general, better behaved than the arithmetic mean.}
#' \item{"riemannian": }{The geometric mean proposed by Barachant et al. (2013).
#' Is approximated iteratively, and accepts two additional argument: \code{max.iter},
#' giving the maximum number of iterations (defaults to 20), and \code{tol}, giving
#' the maximum Frobenius norm of the approximation error (defaults to 0.1)}
#' }
#' @return The mean of the matrices in \code{x}.
spd.mean <- function(x, method = 'euclidean', tol = .1, max.iter = 20){
if (!'spd.list' %in% input.type(x)){
return('Invalid input type')
}
if (method == 'euclidean'){
return(Reduce(`+`, x)/length(x))
}
if (method == 'logeuclidean'){
logmean <- Reduce(`+`, lapply(x, spd.logmap))/length(x)
return(expm(logmean))
}
if (method == 'riemannian'){
# Initialize to arithmetic mean
init <- Reduce(`+`, x) / length(x)
# Iterate until convergence
err <- 2*tol
iter <- 1
while (err > tol & iter <= max.iter){
# Project onto tangent space at current mean estimate
proj <- lapply(x, spd.logmap, p = init)
tan.mean <- Reduce(`+`, proj) / length(proj)
tan.mean <- (tan.mean + t(tan.mean))/2
# Update mean estimate
init <- spd.expmap(tan.mean, p = init)
# Calculate error
err <- norm(tan.mean, type = 'F')
iter <- iter + 1
}
if (iter >= max.iter & err >= tol){
warning('Mean may not have converged')
}
return(init)
}
}
#' Vectorize and unvectorize a matrix
#'
#' Function implements several distance measures between symmetric,
#' positive-definite matrices, not all of which are proper metrics.
#'
#' @param x Either a square matrix, or a vector contains the lower triangular
#' elements of x (including the diagonal elements)
#' @param scaling If TRUE, scales the off diagonal elements by sqrt(2) to preserve
#' the norm of the resulting vector
#' @details If input is a square matrix, converts the lower triangular elements
#' (including the diagonal) into a numeric vector. If input is a numeric vector,
#' converts the input to a symmetric matrix. Note that, if the input is a vector,
#' its length must be a triangular number.
spd.vectorize <- function(x, scaling = F){
# If input is a matrix, vectorize
if (is.matrix(x)){
if (dim(x)[1] != dim(x)[2]){
stop('Inout must be a square matrix')
}
n <- dim(x)[1]
if (scaling){
x[lower.tri(x)] <- sqrt(2) * x[lower.tri(x)]
}
return(x[lower.tri(x, diag = T)])
# If input is a vector, check triangular, then vectorize
} else if (is.vector(x)) {
n <- length(x)
if (sqrt(8*n + 1) %% 1 != 0){
stop('Length of x must be a triangular number')
}
k <- (sqrt(8*n + 1) - 1)/2
y <- matrix(0, k, k)
y[lower.tri(y, diag = T)] <- x
if (scaling){
y[lower.tri(y)] <- y[lower.tri(y)]
}
return(y + t(y) - diag(diag(y)))
}
}
#' Whiten a symmetric positive definite matrix
#'
#' Approximately whitens an spd matrix x by a matrix p. Specifically,
#' conjugates x by inverse.sqrt(p), which, assuming that p is relatively close
#' to x, is close to the identity matrix.
#'
#' @param x,p Symmetric, positive-definite matrices
#' @param dewhiten If TRUE, reverse the whitening transform
spd.whiten <- function(x, p, dewhiten = F) {
# Check input
if (!is.spd(x) | !is.spd(p)){
stop('x must be a symmetric matrix')
}
# Whiten
if (!dewhiten){
p.inv.sqrt <- solve(sqrt(p))
return( p.inv.sqrt %*% x %*% p.inv.sqrt )
} else {
p.sqrt <- sqrt(p)
return( p.sqrt %*% x %*% p.sqrt )
}
}
# spdm
Routines for the analysis of symmetric, positive-definite matrices
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/is.spd.R
\name{is.spd}
\alias{is.spd}
\title{Check if a matrix is symmetric and positive definite}
\usage{
is.spd(..., tol = 1e-08)
}
\arguments{
\item{...}{One or more numeric matrices}
\item{tol}{Numerical tolerance for checking symmetry and positive definiteness}
}
\value{
A Boolean value indicating whether the matrix is symmetric and positive
definite
}
\description{
Check if a matrix is symmetric and positive definite
}
\details{
Function checks whether the matrix \code{x} is symmetric and positive
definite. Symmetry is evaluated up to an entrywise tolerance of \code{tol}, so
that differences smaller than \code{tol} are ignored. Positive-definiteness is
checked by computing the eigenvalues of \code{x} using \code{eigen}, setting
eigenvalues smaller than \code{tol} in absolute value to zero, and then checking
whether any are less than or equal to zero.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/spd.center.R
\name{spd.center}
\alias{spd.center}
\title{Center a set of spd matrices}
\usage{
spd.center(x, y = NULL)
}
\arguments{
\item{x}{Either a list of symmetric, positive-definite matrices, an individual amtrix}
\item{y}{Value(s) for centering. Either a list of spd matrices of the same
length as x, a single spd matrix, or NULL. If NULL, uses the mean of \code{x}}
}
\value{
An object of the same type as \code{x}, containing the centered matrices..
}
\description{
Functions centers a set of symmetric, positive-definite matrices so that they
have mean equal to the identity matrix.
}
\details{
Functions transports the matrices in \code{x} so that \code{y} lies at
the identity. If \code{x} is a list and \code{y = NULL} (default), function sets
\code{y} equal to the mean of \code{x}. Otherwise, \code{y} can be set to any
spd matrix. When \code{x} is a matrix, \code{y} must be provided. If \code{y}
is a list of the same length as {x}, the function centers each element of
\code{x} by the corresponding element of \code{y}.
Centering is accomplished by the parallel transport specified by XXX (2XXX), by
first computing the inverse square root of \code{y}, and then centering \code{x}
by \code{y.inv.sqrt x y.inv.sqrt}.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/spd.dist.R
\name{spd.dist}
\alias{spd.dist}
\title{Distance between two symmetric, positive-definite matrices}
\usage{
spd.dist(x, y, method = "riemannian")
}
\arguments{
\item{x, y}{Symmetric, positive-definite matrices}
\item{method}{The distance measure. See details.}
}
\description{
Function implements several distance measures between symmetric,
positive-definite matrices, not all of which are proper metrics.
}
\details{
Allowable distance measures are
\itemize{
\item{"frobenius": }{The Frobenius norm of the difference \code{x-y}. Not affinely invariant.}
\item{"cholesky": }{The Frobenius norm of the difference between the cholesky factors
of \code{x} and \code{y}. Not affinely invariant.}
\item{"affine": }{An affinely invariant distance measure, computed by taking
the Frobenius norm of the projection of \code{y} into the tangent space at \code{x}}
\item{"riemannian": }{The Riemmanian distance proposed by Barachant, et al. (2013)}
}
}