FCP_TPA {MFPCA}

The functional CP-TPA algorithm

Description

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

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

where \circ 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 \circ w_k can be regarded as the k-th eigenimage, while d_k \cdot u_k represents the vector of individual scores for this eigenimage and each observation.

Usage

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.

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.

Value

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.

References

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

See Also

fcptpaBasis

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