Nothing
#' The functional CP-TPA algorithm
#'
#' This function implements the functional CP-TPA (FCP-TPA) algorithm, that
#' calculates a smooth PCA for 3D tensor data (i.e. \code{N} observations of 2D
#' images with dimension \code{S1 x S2}). The results are given in a
#' CANDECOMP/PARAFRAC (CP) model format \deqn{X = \sum_{k = 1}^K d_k \cdot u_k
#' \circ v_k \circ w_k}{X = \sum d_k u_k \%o\% v_k \%o\% w_k} where
#' \eqn{\circ}{\%o\%} stands for the outer product, \eqn{d_k} is a scalar and
#' \eqn{u_k, v_k, w_k} are eigenvectors for each direction of the tensor. In
#' this representation, the outer product \eqn{v_k \circ w_k}{v_k \%o\% w_k} can
#' be regarded as the \eqn{k}-th eigenimage, while \eqn{d_k \cdot u_k}{d_k u_k}
#' represents the vector of individual scores for this eigenimage and each
#' observation.
#'
#' The smoothness of the eigenvectors \eqn{v_k, w_k} is induced by penalty
#' matrices for both image directions, that are weighted by smoothing parameters
#' \eqn{\alpha_{vk}, \alpha_{wk}}. The eigenvectors \eqn{u_k} are not smoothed,
#' hence the algorithm does not induce smoothness along observations.
#'
#' Optimal smoothing parameters are found via a nested generalized cross
#' validation. In each iteration of the TPA (tensor power algorithm), the GCV
#' criterion is optimized via \code{\link[stats]{optimize}} on the interval
#' specified via \code{alphaRange$v} (or \code{alphaRange$w}, respectively).
#'
#' The FCP_TPA algorithm is an iterative algorithm. Convergence is assumed if
#' the relative difference between the actual and the previous values are all
#' below the tolerance level \code{tol}. The tolerance level is increased
#' automatically, if the algorithm has not converged after \code{maxIter} steps
#' and if \code{adaptTol = TRUE}. If the algorithm did not converge after
#' \code{maxIter} steps (or \code{2 * maxIter}) steps, the function throws a
#' warning.
#'
#' @param X The data tensor of dimensions \code{N x S1 x S2}.
#' @param K The number of eigentensors to be calculated.
#' @param penMat A list with entries \code{v} and \code{w}, containing a
#' roughness penalty matrix for each direction of the image. The algorithm
#' does not induce smoothness along observations (see Details).
#' @param alphaRange A list of length 2 with entries \code{v} and \code{w} ,
#' containing the range of smoothness parameters to test for each direction.
#' @param verbose Logical. If \code{TRUE}, computational details are given on
#' the standard output during calculation of the FCP_TPA.
#' @param tol A numeric value, giving the tolerance for relative error values in
#' the algorithm. Defaults to \code{1e-4}. It is automatically multiplied by
#' 10 after \code{maxIter} steps, if \code{adaptTol = TRUE}.
#' @param maxIter A numeric value, the maximal iteration steps. Can be doubled,
#' if \code{adaptTol = TRUE}.
#' @param adaptTol Logical. If \code{TRUE}, the tolerance is adapted (multiplied
#' by 10), if the algorithm has not converged after \code{maxIter} steps and
#' another \code{maxIter} steps are allowed with the increased tolerance, see
#' Details. Use with caution. Defaults to \code{TRUE}.
#'
#' @return \item{d}{A vector of length \code{K}, containing the numeric weights
#' \eqn{d_k} in the CP model.} \item{U}{A matrix of dimensions \code{N x K},
#' containing the eigenvectors \eqn{u_k} in the first dimension.} \item{V}{A
#' matrix of dimensions \code{S1 x K}, containing the eigenvectors \eqn{v_k}
#' in the second dimension.} \item{W}{A matrix of dimensions \code{S2 x K},
#' containing the eigenvectors \eqn{w_k} in the third dimension.}
#'
#' @references G. I. Allen, "Multi-way Functional Principal Components
#' Analysis", IEEE International Workshop on Computational Advances in
#' Multi-Sensor Adaptive Processing, 2013.
#'
#' @seealso \code{\link{fcptpaBasis}}
#'
#' @importFrom stats runif
#'
#' @export FCP_TPA
#'
#' @examples
#' # set.seed(1234)
#'
#' N <- 100
#' S1 <- 75
#' S2 <- 75
#'
#' # define "true" components
#' v <- sin(seq(-pi, pi, length.out = S1))
#' w <- exp(seq(-0.5, 1, length.out = S2))
#'
#' # simulate tensor data with dimensions N x S1 x S2
#' X <- rnorm(N, sd = 0.5) %o% v %o% w
#'
#' # create penalty matrices (penalize first differences for each dimension)
#' Pv <- crossprod(diff(diag(S1)))
#' Pw <- crossprod(diff(diag(S2)))
#'
#' # estimate one eigentensor
#' res <- FCP_TPA(X, K = 1, penMat = list(v = Pv, w = Pw),
#' alphaRange = list(v = c(1e-4, 1e4), w = c(1e-4, 1e4)),
#' verbose = TRUE)
#'
#' # plot the results and compare to true values
#' plot(res$V)
#' points(v/sqrt(sum(v^2)), pch = 20)
#' legend("topleft", legend = c("True", "Estimated"), pch = c(20, 1))
#'
#' plot(res$W)
#' points(w/sqrt(sum(w^2)), pch = 20)
#' legend("topleft", legend = c("True", "Estimated"), pch = c(20, 1))
FCP_TPA <- function(X, K, penMat, alphaRange, verbose = FALSE, tol = 1e-4, maxIter = 15, adaptTol = TRUE)
{
### check input parameters
check_FCP_TPA_input(X, K, penMat, alphaRange, verbose, tol, maxIter, adaptTol)
dimX <- dim(X)
### initialize all relevant values
# initialize the norm-one vectors u, v, w randomly
u <- stats::runif(dimX[1], min = -1, max = 1); u <- u/normVec(u)
v <- stats::runif(dimX[2], min = -1, max = 1); v <- v/normVec(v)
w <- stats::runif(dimX[3], min = -1, max = 1); w <- w/normVec(w)
# initialize smoothing parameters: start with minimal values (almost no smoothing)
alphaV <- min(alphaRange$v)
alphaW <- min(alphaRange$w)
# eigendecomposition of penalty matrices
eV <- eigen(penMat$v)
eW <- eigen(penMat$w)
# digonal matrices
Iv <- diag(dimX[2])
Iw <- diag(dimX[3])
# initialize results
d <- rep(NA, K)
U <- matrix(NA, ncol = K, nrow = dimX[1])
V <- matrix(NA, ncol = K, nrow = dimX[2])
W <- matrix(NA, ncol = K, nrow = dimX[3])
for(k in seq_len(K))
{
if(verbose)
cat("\nk = ", k, "\n")
# initialize "old" versions
uOld <- 0*u
vOld <- 0*v
wOld <- 0*w
tolOld <- tol
# number of iterations
iter <- 0
# repeat until convergence
while(any(c(normVec(u - uOld)/normVec(u), normVec(v - vOld)/normVec(v), normVec(w - wOld)/normVec(w)) > tol))
{
# update u
uOld <- u
u <- ttv(X, list(v, w), dim = 2:3) /
as.numeric(crossprod(v, v + alphaV*penMat$v %*%v )*crossprod(w, w + alphaW * penMat$w %*% w))
# update v
vOld <- v
v <- solve(Iv + alphaV * penMat$v, ttv(X, list(u, w), dim = c(1,3)))/
as.numeric(crossprod(u, u)*crossprod(w, w + alphaW * penMat$w %*% w))
# update alphaV
alphaV <- findAlphaVopt(alphaRange = alphaRange$v, data = X, u = u, w = w,
alphaW = alphaW, OmegaW = penMat$w,
GammaV = eV$vectors, lambdaV = eV$values)
# update w
wOld <- w
w <- solve(Iw + alphaW * penMat$w, ttv(X, list(u, v), dim = 1:2)) /
as.numeric(crossprod(u, u)* crossprod(v, v + alphaV*penMat$v %*% v ))
# update alphaW
alphaW <- findAlphaWopt(alphaRange = alphaRange$w, data = X, u = u, v = v,
alphaV = alphaV, OmegaV = penMat$v,
GammaW = eW$vectors, lambdaW = eW$values)
iter <- iter + 1
if(iter > maxIter)
{
if((iter < 2*maxIter) & adaptTol) # try another maxIter with increased tolerance
{
tol <- 10*tol
}
else
{
uOld <- u
vOld <- v
wOld <- w
warning("FCP-TPA algorithm did not converge; iteration ", k , " stopped.")
}
}
}
if(verbose)
cat("Absolute error:\n
u: ", normVec(u - uOld), ", v: ", normVec(v - vOld), ", w: ", normVec(w - wOld), ", alphaV: ", alphaV, ", alphaW: ", alphaW, "\n")
# reset tolerance if necessary
if(adaptTol & (iter >= maxIter) )
tol <- tolOld
# scale vectors to have norm one
u <- u/normVec(u)
v <- v/normVec(v)
w <- w/normVec(w)
# calculcate results
d[k] <- ttv(X, list(u,v,w), dim = 1:3)
U[,k] <- u
V[,k] <- v
W[,k] <- w
# update X
X <- X - d[k] * u %o% v %o% w #u[,1] %o% v[,1] %o% w[,1] # outer product
}
return(list(d = d, U = U, V = V, W = W))
}
#' Check input of FCP TPA function
#'
#' @return No return value, called for side effects
#'
#' @keywords internal
check_FCP_TPA_input <- function(X, K, penMat, alphaRange, verbose, tol, maxIter, adaptTol)
{
# X
if(!is.array(X))
stop("Parameter 'X' must be passed as an array.")
else
{
if(length(dim(X)) != 3)
stop("Parameter 'X' must have three dimensions.")
}
dimX <- dim(X)
# K
if(! (is.numeric(K) & length(K) == 1))
stop("Parameter 'K' must be passed as a number.")
# penMat
if(!is.list(penMat))
stop("Parameter 'penMat' must be passed as a list.")
else
{
if(!identical(names(penMat), c("v", "w")))
stop("Function FCP_TPA: penMat must be a list of matrices with entries 'v' and 'w'.")
# penMat$v: correct dimensions and symmetric
if(any(dim(penMat$v) != rep(dimX[2],2)))
stop("Function FCP_TPA: the penalization matrix for dimension v must be ", dimX[2], " x ", dimX[2], ".")
if(!isSymmetric(penMat$v))
stop("Function FCP_TPA: the penalization matrix for dimension v must be symmetric.")
# penMat$w: correct dimensions and symmetric
if(any(dim(penMat$w) != rep(dimX[3],2)))
stop("Function FCP_TPA: the penalization matrix for dimension w must be ", dimX[3], " x ", dimX[3], ".")
if(!isSymmetric(penMat$w))
stop("Function FCP_TPA: the penalization matrix for dimension w must be symmetric.")
}
# alphaRange
if(!is.list(alphaRange))
stop("Parameter 'alphaRange' must be passed as a list.")
else
{
if(!identical(names(alphaRange), c("v", "w")))
stop("Function FCP_TPA: alphaRange must be a list of vectors with entries 'v' and 'w'.")
if(length(alphaRange$v) != 2)
stop("Function FCP_TPA: alphaRange$v must be a vector of length 2.")
if(any(alphaRange$v < 0))
stop("Function FCP_TPA: Values for alphaV must not be negative.")
if(length(alphaRange$w) != 2)
stop("Function FCP_TPA: alphaRange$w must be a vector of length 2.")
if(any(alphaRange$w < 0))
stop("Function FCP_TPA: Values for alphaW must not be negative.")
}
# verbose
if(!is.logical(verbose))
stop("Parameter 'verbose' must be passed as a logical.")
# tol
if(! (is.numeric(tol) & length(tol) == 1))
stop("Parameter 'tol' must be passed as a number.")
# maxIter
if(! (is.numeric(maxIter) & length(maxIter) == 1))
stop("Parameter 'maxIter' must be passed as a number.")
# adaptTol
if(!is.logical(adaptTol))
stop("Parameter 'adaptTol' must be passed as a logical.")
return()
}
#' Calculate the euclidean norm of a vector
#'
#' @param x The vector for which the norm is to be calculated
#'
#' @return The euclidean norm of \code{x}.
#'
#' @keywords internal
normVec <- function(x)
{
sqrt(sum(x^2))
}
#' Find the optimal smoothing parameters in FCP_TPA using GCV
#'
#' These functions find the optimal smoothing parameters \eqn{\alpha_v,
#' \alpha_w} for the two image directions (v and w) in the FCP_TPA algorithm
#' based on generalized cross-validation, which is nested in the tensor power
#' algorithm. Given a range of possible values of \eqn{\alpha_v} (or
#' \eqn{\alpha_w}, respectively), the optimum is found by optimizing the GCV
#' criterion using the function \code{\link[stats]{optimize}}.
#'
#' @param alphaRange A numeric vector with two elements, containing the minimal
#' and maximal value for the smoothing parameter that is to be optimized.
#' @param data The tensor containing the data, an array of dimensions \code{N x
#' S1 x S2}.
#' @param u,v,w The current value of the eigenvectors \eqn{u_k, v_k, w_k} (not
#' normalized) of dimensions \code{N}, \code{S1} and \code{S2}.
#' @param alphaV,alphaW The current value of the smoothing parameter for the
#' other image direction (\eqn{\alpha_w} for \code{findAlphaVopt} and
#' \eqn{\alpha_v} for \code{findAlphaWopt}), which is kept as fixed.
#' @param OmegaV, OmegaW A matrix of dimension \code{S1 x S1} (\code{OmegaV} in \code{findAlphaWopt}) or \code{S2 x S2} (\code{OmegaW} in \code{findAlphaVopt}), the penalty matrix for
#' other image direction.
#' @param GammaV,GammaW A matrix of dimension \code{S1 x S1} (\code{GammaV} in \code{findAlphaVopt}) or \code{S2 x S2} (\code{GammaW} in \code{findAlphaWopt}), containing the
#' eigenvectors of the penalty matrix for the image direction for which the optimal smoothing parameter is to be found.
#' @param lambdaV, lambdaW A numeric vector of length \code{S1}(\code{lambdaV} in \code{findAlphaVopt}) or \code{S2} (\code{lambdaW} in \code{findAlphaWopt}), containing the
#' eigenvalues of the penalty matrix for the image direction for which the optimal smoothing parameter is to be found.
#'
#' @return The optimal \eqn{\alpha_v} (or \eqn{\alpha_w}, respectively), found by optimizing the GCV criterion
#' within the given range of possible values.
#'
#' @references G. I. Allen (2013), "Multi-way Functional Principal Components
#' Analysis", IEEE International Workshop on Computational Advances in
#' Multi-Sensor Adaptive Processing.
#' @references J. Z. Huang, H. Shen and A. Buja (2009), "The Analysis of Two-Way
#' Functional Data Using Two-Way Regularized Singular Value Decomposition".
#' Journal of the American Statistical Association, Vol. 104, No. 488, 1609 --
#' 1620.
#'
#' @keywords internal
#'
#' @seealso \code{\link{FCP_TPA}}, \code{\link{gcv}}
#'
#' @importFrom stats optimize
findAlphaVopt <- function(alphaRange, data, u, w, alphaW, OmegaW, GammaV, lambdaV)
{
z <- crossprod(GammaV, ttv(data, list(u,w), c(1,3))) / (normVec(u) * normVec(w))
wOw <- as.numeric(crossprod(w, OmegaW %*% w)) # without as.numeric, this is a 1x1 matrix
eta <- 1/(1 + alphaW * wOw/normVec(w))
res <- stats::optimize(f=gcv, interval=c(min(alphaRange), max(alphaRange)),
n = length(lambdaV), z = z, eta = eta, lambda = lambdaV)$minimum
return(res)
}
#' @describeIn findAlphaVopt
#'
#' @importFrom stats optimize
findAlphaWopt <- function(alphaRange, data, u, v, alphaV, OmegaV, GammaW, lambdaW)
{
z <- crossprod(GammaW, ttv(data, list(u,v), c(1,2))) / (normVec(u) * normVec(v))
vOv <- as.numeric(crossprod(v, OmegaV %*% v))
eta <- 1/(1 + alphaV * vOv/normVec(v)) # without as.numeric, this is a 1x1 matrix
res <- stats::optimize(f=gcv, interval=c(min(alphaRange), max(alphaRange)),
n = length(lambdaW), z = z, eta = eta, lambda = lambdaW)$minimum
return(res)
}
#' Generalized cross-validation for the FCP-TPA algorithm
#'
#' These function calculates the generalized cross-validation criterion for the
#' smoothing parameters \eqn{\alpha_v} or \eqn{\alpha_w} that are used in the
#' \code{\link{FCP_TPA}} algorithm. As the criterion is symmetric in \eqn{v} and
#' \eqn{w}, this function implements a generic criterion, which is called by
#' \code{\link{findAlphaVopt}}, \code{\link{findAlphaWopt}} with the correct values.
#'
#' The criterion can be evaluated in a numerically efficient way, adopting the ideas in Huang, Shen and Buja (2008) to three-ways tensors. TODO!
#'
#' @param alpha The current value of the smoothing parameter.
#' @param n The length of the dimension, for which the smoothing parameter is to
#' be optimized.
#' @param z A vector of length \code{n}. See Details.
#' @param eta A vector of length \code{n}. See Details.
#' @param lambda A vector of length \code{n}, containing the eigenvalues of the
#' penalty matrix corresponding the the current image direction.
#'
#' @return The value of the GCV criterion.
#'
#' @references G. I. Allen, "Multi-way Functional Principal Components
#' Analysis", IEEE International Workshop on Computational Advances in
#' Multi-Sensor Adaptive Processing, 2013.
#'
#' @keywords internal
#'
#' @seealso \code{\link{FCP_TPA}}
gcv <- function(alpha, n, z, eta, lambda)
{
eS <- 1/(1 + alpha*lambda) # eigenvalues of S
res <- sum((z* (1 - eta* eS))^2)/n
res <- res/(1 - eta/n * sum(eS))^2
return(res)
}
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.