FCP_TPA() R function from [MFPCA]

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. N observations of 2D images with dimension S1 x S2). The results are given in a CANDECOMP/PARAFRAC (CP) model format [REMOVE_ME] $X = \sum_{k = 1}^K d_k \cdot u_k\circ v_k \circ w_kX = \sum d_k u_k %o% v_k %o% w_k [REMOVE_ME_2]$ where $%o%$ stands for the outer product, $d_k$ is a scalar and $u_k, v_k, w_k$ are eigenvectors for each direction of the tensor. In this representation, the outer product $v_k %o% w_k$ can be regarded as the $k$ -th eigenimage, while $d_k u_k$

represents the vector of individual scores for this eigenimage and each observation.


FCP_TPA(
  X,
  K,
  penMat,
  alphaRange,
  verbose = FALSE,
  tol = 1e-04,
  maxIter = 15,
  adaptTol = TRUE
)

Arguments

X: The data tensor of dimensions N x S1 x S2.
K: The number of eigentensors to be calculated.
penMat: A list with entries v and w, containing a roughness penalty matrix for each direction of the image. The algorithm does not induce smoothness along observations (see Details).
alphaRange: A list of length 2 with entries v and w , containing the range of smoothness parameters to test for each direction.
verbose: Logical. If TRUE, computational details are given on the standard output during calculation of the FCP_TPA.
tol: A numeric value, giving the tolerance for relative error values in the algorithm. Defaults to 1e-4. It is automatically multiplied by 10 after maxIter steps, if adaptTol = TRUE.
maxIter: A numeric value, the maximal iteration steps. Can be doubled, if adaptTol = TRUE.
adaptTol: Logical. If TRUE, the tolerance is adapted (multiplied by 10), if the algorithm has not converged after maxIter steps and another maxIter steps are allowed with the increased tolerance, see Details. Use with caution. Defaults to TRUE.

Returns

d: A vector of length K, containing the numeric weights $d_k$ in the CP model. - U: A matrix of dimensions N x K, containing the eigenvectors $u_k$ in the first dimension. - V: A matrix of dimensions S1 x K, containing the eigenvectors $v_k$

in the second dimension. - W: A matrix of dimensions S2 x K, containing the eigenvectors $w_k$ in the third dimension.

Description

X = \sum_{k = 1}^K d_k \cdot u_k\circ v_k \circ w_kX = \sum d_k u_k %o% v_k %o% w_k

where $%o%$ stands for the outer product, $d_k$ is a scalar and $u_k, v_k, w_k$ are eigenvectors for each direction of the tensor. In this representation, the outer product $v_k %o% w_k$ can be regarded as the $k$ -th eigenimage, while $d_k u_k$

represents the vector of individual scores for this eigenimage and each observation.

Details

The smoothness of the eigenvectors $v_k, w_k$ is induced by penalty matrices for both image directions, that are weighted by smoothing parameters $\alpha_{vk}, \alpha_{wk}$ . The eigenvectors $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 optimize on the interval specified via alphaRange$v (or 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 tol. The tolerance level is increased automatically, if the algorithm has not converged after maxIter steps and if adaptTol = TRUE. If the algorithm did not converge after maxIter steps (or 2 * maxIter) steps, the function throws a warning.

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))

References

G. I. Allen, "Multi-way Functional Principal Components Analysis", IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2013.

FCP_TPA function