###Functions that operate on Matrices and Arrays
#'List Hamadard Product
#'
#'Returns the Hamadard (element-wise) product from a list of matrices or vectors. Commonly used for n-mode products and various Tensor decompositions.
#'@name hamadard_list
#'@rdname hamadard_list
#'@aliases hamadard_list
#'@export
#'@param L list of matrices or vectors
#'@return matrix that is the Hamadard product
#'@seealso \code{\link{kronecker_list}}, \code{\link{khatri_rao_list}}
#'@note The modes/dimensions of each element in the list must match.
#'@examples
#'lizt <- list('mat1' = matrix(runif(40),ncol=4),
#' 'mat2' = matrix(runif(40),ncol=4),
#' 'mat3' = matrix(runif(40),ncol=4))
#'dim(hamadard_list(lizt))
hamadard_list <- function(L){
isvecORmat <- function(x){is.matrix(x) || is.vector(x)}
stopifnot(all(unlist(lapply(L,isvecORmat))))
retmat <- L[[1]]
for (i in 2:length(L)){
retmat <- retmat*L[[i]]
}
retmat
}
#'List Kronecker Product
#'
#'Returns the Kronecker product from a list of matrices or vectors. Commonly used for n-mode products and various Tensor decompositions.
#'@name kronecker_list
#'@rdname kronecker_list
#'@aliases kronecker_list
#'@export
#'@param L list of matrices or vectors
#'@return matrix that is the Kronecker product
#'@seealso \code{\link{hamadard_list}}, \code{\link{khatri_rao_list}}, \code{\link{kronecker}}
#'@examples
#'smalllizt <- list('mat1' = matrix(runif(12),ncol=4),
#' 'mat2' = matrix(runif(12),ncol=4),
#' 'mat3' = matrix(runif(12),ncol=4))
#'dim(kronecker_list(smalllizt))
kronecker_list <- function(L){
isvecORmat <- function(x){is.matrix(x) || is.vector(x)}
stopifnot(all(unlist(lapply(L,isvecORmat))))
retmat <- L[[1]]
for(i in 2:length(L)){
retmat <- kronecker(retmat,L[[i]])
}
retmat
}
#'Khatri-Rao Product
#'
#'Returns the Khatri-Rao (column-wise Kronecker) product of two matrices. If the inputs are vectors then this is the same as the Kronecker product.
#'@name khatri_rao
#'@rdname khatri_rao
#'@aliases khatri_rao
#'@export
#'@param x first matrix
#'@param y second matrix
#'@return matrix that is the Khatri-Rao product
#'@seealso \code{\link{kronecker}}, \code{\link{khatri_rao_list}}
#'@note The number of columns must match in the two inputs.
#'@examples
#'dim(khatri_rao(matrix(runif(12),ncol=4),matrix(runif(12),ncol=4)))
khatri_rao <- function(x,y){
if (!(is.matrix(x)&&is.matrix(y))) stop("Arguments must be matrices.")
if (dim(x)[2]!=dim(y)[2]) stop("Arguments must have same number of columns.")
retmat <- matrix(0,nrow=dim(x)[1]*dim(y)[1],ncol=dim(x)[2])
for (j in 1:ncol(retmat)) retmat[,j] <- kronecker(x[,j],y[,j])
retmat
}
#'List Khatri-Rao Product
#'
#'Returns the Khatri-Rao product from a list of matrices or vectors. Commonly used for n-mode products and various Tensor decompositions.
#'@name khatri_rao_list
#'@rdname khatri_rao_list
#'@aliases khatri_rao_list
#'@export
#'@param L list of matrices or vectors
#'@param reverse whether or not to reverse the order
#'@return matrix that is the Khatri-Rao product
#'@seealso \code{\link{khatri_rao}}
#'@note The number of columns must match in every element of the input list.
#'@examples
#'smalllizt <- list('mat1' = matrix(runif(12),ncol=4),
#' 'mat2' = matrix(runif(12),ncol=4),
#' 'mat3' = matrix(runif(12),ncol=4))
#'dim(khatri_rao_list(smalllizt))
khatri_rao_list <- function(L,reverse=FALSE){
stopifnot(all(unlist(lapply(L,is.matrix))))
ncols <- unlist(lapply(L,ncol))
stopifnot(length(unique(ncols))==1)
ncols <- ncols[1]
nrows <- unlist(lapply(L,nrow))
retmat <- matrix(0,nrow=prod(nrows),ncol=ncols)
if (reverse) L <- rev(L)
for(j in 1:ncols){
Lj <- lapply(L,function(x) x[,j])
retmat[,j] <- kronecker_list(Lj)
}
retmat
}
#'Tensor Times Matrix (m-Mode Product)
#'
#'Contracted (m-Mode) product between a Tensor of arbitrary number of modes and a matrix. The result is folded back into Tensor.
#'@name ttm
#'@rdname ttm
#'@aliases ttm
#'@details By definition, \code{rs_unfold(ttm(tnsr,mat),m) = mat\%*\%rs_unfold(tnsr,m)}, so the number of columns in \code{mat} must match the \code{m}th mode of \code{tnsr}. For the math on the m-Mode Product, see Kolda and Bader (2009).
#'@export
#'@param tnsr Tensor object with K modes
#'@param mat input matrix with same number columns as the \code{m}th mode of \code{tnsr}
#'@param m the mode to contract on
#'@param transpose if mat should be transposed before multiplication
#'@return a Tensor object with K modes
#'@seealso \code{\link{ttl}}, \code{\link{rs_unfold-methods}}
#'@note The \code{m}th mode of \code{tnsr} must match the number of columns in \code{mat}. By default, the returned Tensor does not drop any modes equal to 1.
#'@references T. Kolda, B. Bader, "Tensor decomposition and applications". SIAM Applied Mathematics and Applications 2009.
#'@examples
#'tnsr <- new("Tensor",3L,c(3L,4L,5L),data=runif(60))
#'mat <- matrix(runif(50),ncol=5)
#'ttm(tnsr,mat,m=3)
ttm<-function(tnsr,mat,m=NULL,transpose=FALSE){
stopifnot(is.matrix(mat))
if(is.null(m)) stop("m must be specified")
mat_dims <- dim(mat)
if ( transpose ) mat_dims <- rev(mat_dims)
modes_in <- tnsr@modes
stopifnot(modes_in[m]==mat_dims[2])
modes_out <- modes_in
modes_out[m] <- mat_dims[1]
tnsr_m <- rs_unfold(tnsr,m=m)@data
retarr_m <- if (transpose) crossprod(mat,tnsr_m) else mat%*%tnsr_m
rs_fold(retarr_m,m=m,modes=modes_out)
}
#'Tensor Times List
#'
#'Contracted (m-Mode) product between a Tensor of arbitrary number of modes and a list of matrices. The result is folded back into Tensor.
#'@name ttl
#'@rdname ttl
#'@aliases ttl
#'@details Performs \code{ttm} repeated for a single Tensor and a list of matrices on multiple modes. For instance, suppose we want to do multiply a Tensor object \code{tnsr} with three matrices \code{mat1}, \code{mat2}, \code{mat3} on modes 1, 2, and 3. We could do \code{ttm(ttm(ttm(tnsr,mat1,1),mat2,2),3)}, or we could do \code{ttl(tnsr,list(mat1,mat2,mat3),c(1,2,3))}. The order of the matrices in the list should obviously match the order of the modes. This is a common operation for various Tensor decompositions such as CP and Tucker. For the math on the m-Mode Product, see Kolda and Bader (2009).
#'@export
#'@param tnsr Tensor object with K modes
#'@param list_mat a list of matrices
#'@param ms a vector of modes to contract on (order should match the order of \code{list_mat})
#'@param transpose if matrices should be transposed before multiplication
#'@return Tensor object with K modes
#'@seealso \code{\link{ttm}}
#'@note The returned Tensor does not drop any modes equal to 1.
#'@references T. Kolda, B. Bader, "Tensor decomposition and applications". SIAM Applied Mathematics and Applications 2009.
#'@examples
#'tnsr <- new("Tensor",3L,c(3L,4L,5L),data=runif(60))
#'lizt <- list('mat1' = matrix(runif(30),ncol=3),
#' 'mat2' = matrix(runif(40),ncol=4),
#' 'mat3' = matrix(runif(50),ncol=5))
#'ttl(tnsr,lizt,ms=c(1,2,3))
ttl<-function(tnsr,list_mat,ms=NULL,transpose=FALSE){
if(is.null(ms)||!is.vector(ms)) stop ("m modes must be specified as a vector")
if(length(ms)!=length(list_mat)) stop("m modes length does not match list_mat length")
num_mats <- length(list_mat)
if(length(unique(ms))!=num_mats) warning("consider pre-multiplying matrices for the same m for speed")
mat_nrows <- vector("list", num_mats)
mat_ncols <- vector("list", num_mats)
for(i in 1:num_mats){
mat <- list_mat[[i]]
m <- ms[i]
mat_dims <- dim(mat)
if (transpose) mat_dims <- rev(mat_dims)
modes_in <- tnsr@modes
stopifnot(modes_in[m]==mat_dims[2])
modes_out <- modes_in
modes_out[m] <- mat_dims[1]
tnsr_m <- rs_unfold(tnsr,m=m)@data
retarr_m <- if (transpose) crossprod(mat, tnsr_m) else mat%*%tnsr_m
tnsr <- rs_fold(retarr_m,m=m,modes=modes_out)
}
tnsr
}
#'Tensor Multiplication (T-MULT)
#'
#'Implements T-MULT based on block circulant matrices (Kilmer et al. 2013) for 3-tensors.
#'
#'@details Uses the Fast Fourier Transform (FFT) speed up suggested by Kilmer et al. 2013 instead of explicitly constructing the block circulant matrix. For the mathematical details of T-MULT, see Kilmer et al. (2013).
#'@export
#'@name t_mult
#'@rdname t_mult
#'@aliases t_mult
#'@param x a 3-tensor
#'@param y another 3-tensor
#'@return tensor product between \code{x} and \code{y}
#'@note This only works (so far) between 3-Tensors.
#'@references M. Kilmer, K. Braman, N. Hao, and R. Hoover, "Third-order tensors as operators on matrices: a theoretical and computational framework with applications in imaging". SIAM Journal on Matrix Analysis and Applications 2013.
#'@examples
#'tnsr <- new("Tensor",3L,c(3L,4L,5L),data=runif(60))
#'tnsr2 <- new("Tensor",3L,c(4L,3L,5L),data=runif(60))
#'t_mult(tnsr, tnsr2)
t_mult <- function(x,y){
if((x@num_modes>3)||(y@num_modes>3)) stop("Tensor Multiplication currently only implemented for 3-Tensors")
modes_x <- x@modes
modes_y <- y@modes
if(modes_x[2]!=modes_y[1]) stop("Mode 2 of x and Mode 1 of y must match")
n3 <- modes_x[3]
if(n3!=modes_y[3]) stop("Modes 3 of x and y must match")
#fft's for x and y
fft_x <- aperm(apply(x@data,MARGIN=1:2,fft),c(2,3,1))
fft_y <- aperm(apply(y@data,MARGIN=1:2,fft),c(2,3,1))
#multiply the faces (this is terribad! TO-DO: think of better way!)
fft_ret <- array(0,dim=c(modes_x[1],modes_y[2],n3))
for(i in 1:n3){
first <- fft_x[,,i,drop=FALSE]
second <- fft_y[,,i,drop=FALSE]
fft_ret[,,i]<-matrix(first,nrow=dim(first)[1])%*%matrix(second,,nrow=dim(second)[1])
}
#ifft and return as Tensor
ifft <- function(x){suppressWarnings(as.numeric(fft(x,inverse=TRUE))/length(x))}
as.tensor(aperm(apply(fft_ret,MARGIN=1:2,ifft),c(2,3,1)),drop=FALSE)
}
#####Special Tensors
#'Tensor with Random Entries
#'
#'Generate a Tensor with specified modes with iid normal(0,1) entries.
#'@export
#'@name rand_tensor
#'@rdname rand_tensor
#'@aliases rand_tensor
#'@param modes the modes of the output Tensor
#'@param drop whether or not modes equal to 1 should be dropped
#'@return a Tensor object with modes given by \code{modes}
#'@note Default \code{rand_tensor()} generates a 3-Tensor with modes \code{c(3,4,5)}.
#'@examples
#'rand_tensor()
#'rand_tensor(c(4,4,4))
#'rand_tensor(c(10,2,1),TRUE)
rand_tensor <- function(modes=c(3,4,5),drop=FALSE){
as.tensor(array(rnorm(prod(modes)), dim=modes),drop=drop)
}
###Matrix Foldings
#'General Folding of Matrix
#'
#'General folding of a matrix into a Tensor. This is designed to be the inverse function to \code{\link{unfold-methods}}, with the same ordering of the indices. This amounts to following: if we were to unfold a Tensor using a set of \code{row_idx} and \code{col_idx}, then we can fold the resulting matrix back into the original Tensor using the same \code{row_idx} and \code{col_idx}.
#'@export
#'@details This function uses \code{aperm} as the primary workhorse.
#'@name fold
#'@rdname fold
#'@aliases fold
#'@param mat matrix to be folded into a Tensor
#'@param row_idx the indices of the modes that are mapped onto the row space
#'@param col_idx the indices of the modes that are mapped onto the column space
#'@param modes the modes of the output Tensor
#'@return Tensor object with modes given by \code{modes}
#'@seealso \code{\link{unfold-methods}}, \code{\link{k_fold}}, \code{\link{unmatvec}}
#'@references T. Kolda, B. Bader, "Tensor decomposition and applications". SIAM Applied Mathematics and Applications 2009.
#'@examples
#'tnsr <- new("Tensor",3L,c(3L,4L,5L),data=runif(60))
#'matT3<-unfold(tnsr,row_idx=2,col_idx=c(3,1))
#'identical(fold(matT3,row_idx=2,col_idx=c(3,1),modes=c(3,4,5)),tnsr)
fold <- function(mat, row_idx = NULL, col_idx = NULL, modes=NULL){
#checks
rs <- row_idx
cs <- col_idx
if(is.null(rs)||is.null(cs)) stop("row space and col space indices must be specified")
if(is.null(modes)) stop("Tensor modes must be specified")
if(!is(mat,"Tensor")){
if(!is.matrix(mat)) stop("mat must be of class 'matrix'")
}else{
stopifnot(mat@num_modes==2)
mat <- mat@data
}
num_modes <- length(modes)
stopifnot(num_modes==length(rs)+length(cs))
mat_modes <- dim(mat)
if((mat_modes[1]!=prod(modes[rs])) || (mat_modes[2]!=prod(modes[cs]))) stop("matrix nrow/ncol does not match Tensor modes")
#rearranges into array
iperm <- match(1:num_modes,c(rs,cs))
as.tensor(aperm(array(mat,dim=c(modes[rs],modes[cs])),iperm))
}
#'k-mode Folding of Matrix
#'
#'k-mode folding of a matrix into a Tensor. This is the inverse funtion to \code{k_unfold} in the m mode. In particular, \code{k_fold(k_unfold(tnsr, m),m,getModes(tnsr))} will result in the original Tensor.
#'@export
#'@details This is a wrapper function to \code{\link{fold}}.
#'@name k_fold
#'@rdname k_fold
#'@aliases k_fold
#'@param mat matrix to be folded into a Tensor
#'@param m the index of the mode that is mapped onto the row indices
#'@param modes the modes of the output Tensor
#'@return Tensor object with modes given by \code{modes}
#'@references T. Kolda, B. Bader, "Tensor decomposition and applications". SIAM Applied Mathematics and Applications 2009.
#'@seealso \code{\link{k_unfold-methods}}, \code{\link{fold}}, \code{\link{unmatvec}}
#'@examples
#'tnsr <- new("Tensor",3L,c(3L,4L,5L),data=runif(60))
#'matT2<-k_unfold(tnsr,m=2)
#'identical(k_fold(matT2,m=2,modes=c(3,4,5)),tnsr)
k_fold <- function(mat,m=NULL,modes=NULL){
if(is.null(m)) stop("mode m must be specified")
if(is.null(modes)) stop("Tensor modes must be specified")
num_modes <- length(modes)
rs <- m
cs <- (1:num_modes)[-m]
fold(mat,row_idx=rs,col_idx=cs,modes=modes)
}
#'Unmatvec Folding of Matrix
#'
#'The inverse operation to \code{\link{matvec-methods}}, turning a matrix into a Tensor. For a full account of matrix folding/unfolding operations, consult Kolda and Bader (2009).
#'@export
#'@name unmatvec
#'@rdname unmatvec
#'@aliases unmatvec
#'@param mat matrix to be folded into a Tensor
#'@param modes the modes of the output Tensor
#'@return Tensor object with modes given by \code{modes}
#'@references T. Kolda, B. Bader, "Tensor decomposition and applications". SIAM Applied Mathematics and Applications 2009.
#'@seealso \code{\link{matvec-methods}}, \code{\link{fold}}, \code{\link{k_fold}}
#'@examples
#'tnsr <- new("Tensor",3L,c(3L,4L,5L),data=runif(60))
#'matT1<-matvec(tnsr)
#'identical(unmatvec(matT1,modes=c(3,4,5)),tnsr)
unmatvec <- function(mat,modes=NULL){
if(is.null(modes)) stop("Tensor modes must be specified")
num_modes <- length(modes)
cs <- 2
rs <- (1:num_modes)[-2]
fold(mat,row_idx=rs,col_idx=cs,modes=modes)
}
#'Row Space Folding of Matrix
#'
#'DEPRECATED. Please see \code{\link{k_fold}}.
#'@export
#'@param mat matrix to be folded
#'@param m the mode corresponding to rs_unfold
#'@param modes the original modes of the tensor
#'@name rs_fold
#'@rdname rs_fold
#'@aliases rs_fold
rs_fold <- function(mat,m=NULL,modes=NULL){
if(is.null(m)) stop("mode m must be specified")
if(is.null(modes)) stop("Tensor modes must be specified")
num_modes <- length(modes)
rs <- m
cs <- (1:num_modes)[-m]
fold(mat,row_idx=rs,col_idx=cs,modes=modes)
}
#'Column Space Folding of Matrix
#'
#'DEPRECATED. Please see \code{\link{unmatvec}}
#'@export
#'@param mat matrix to be folded
#'@param m the mode corresponding to cs_unfold
#'@param modes the original modes of the tensor
#'@name cs_fold
#'@rdname cs_fold
#'@aliases cs_fold
cs_fold <- function(mat,m=NULL,modes=NULL){
if(is.null(m)) stop("mode m must be specified")
if(is.null(modes)) stop("Tensor modes must be specified")
num_modes <- length(modes)
cs <- m
rs <- (1:num_modes)[-m]
fold(mat,row_idx=rs,col_idx=cs,modes=modes)
}
###Invisible Functions (undocumented)
#Creates a superdiagonal tensor
.superdiagonal_tensor <- function(num_modes,len,elements=1L){
modes <- rep(len,num_modes)
arr <- array(0, dim = modes)
if(length(elements)==1) elements <- rep(elements,len)
for (i in 1:len){
txt <- paste("arr[",paste(rep("i", num_modes),collapse=","),"] <- ", elements[i],sep="")
eval(parse(text=txt))
}
as.tensor(arr)
}
#3-Tensor Kilmer et. al (2013) identity
.identity_tensor3d <- function(modes){
if(length(modes)!=3L) stop("identity tensor only implemented for 3d so far")
n <- modes[1]
stopifnot(n==modes[2])
arr <- array(0,dim=modes)
arr[,,1] <- diag(1,n,n)
as.tensor(arr)
}
#Simple timing functions
.tic <- function (gcFirst = TRUE,overwrite=TRUE) {
if(gcFirst) gc(FALSE)
tic <- proc.time()
ticExists <- ".tic"%in%ls(all.names=TRUE,envir=baseenv())
if(overwrite||!ticExists){
assign(".tic", tic, envir=baseenv())
}
else{
stop("Another timing function running")
}
invisible(tic)
}
.toc <- function (pr=FALSE) {
toc <- proc.time()
tic <- get(".tic", envir=baseenv())
if(pr) print(toc - tic)
invisible(toc - tic)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.