# parafac2: Parallel Factor Analysis-2 In multiway: Component Models for Multi-Way Data

## Description

Given a list of matrices `X[[k]] = matrix(xk,I[k],J)` for `k = seq(1,K)`, the 3-way Parafac2 model (with Mode A nested in Mode C) can be written as

 `X[[k]] = tcrossprod(A[[k]]%*%diag(C[k,]),B) + E[[k]]` subject to `crossprod(A[[k]]) = Phi`

where `A[[k]] = matrix(ak,I[k],R)` are the Mode A (first mode) weights for the `k`-th level of Mode C (third mode), `Phi` is the common crossproduct matrix shared by all `K` levels of Mode C, `B = matrix(b,J,R)` are the Mode B (second mode) weights, `C = matrix(c,K,R)` are the Mode C (third mode) weights, and `E[[k]] = matrix(ek,I[k],J)` is the residual matrix corresponding to `k`-th level of Mode C.

Given a list of arrays `X[[l]] = array(xl,dim=c(I[l],J,K))` for `l = seq(1,L)`, the 4-way Parafac2 model (with Mode A nested in Mode D) can be written as

 `X[[l]][,,k] = tcrossprod(A[[l]]%*%diag(D[l,]*C[k,]),B) + E[[k]]` subject to `crossprod(A[[l]]) = Phi`

`A[[l]] = matrix(al,I[l],R)` are the Mode A (first mode) weights for the `l`-th level of Mode D (fourth mode), `Phi` is the common crossproduct matrix shared by all `L` levels of Mode D, `D = matrix(d,L,R)` are the Mode D (fourth mode) weights, and `E[[l]] = matrix(el,I[l],J,K)` is the residual array corresponding to `l`-th level of Mode D.

Weight matrices are estimated using an alternating least squares algorithm with optional constraints.

## Usage

 ```1 2 3 4 5 6``` ```parafac2(X, nfac, nstart = 10, const = NULL, control = NULL, Gfixed = NULL, Bfixed = NULL, Cfixed = NULL, Dfixed = NULL, Gstart = NULL, Bstart = NULL, Cstart = NULL, Dstart = NULL, Gstruc = NULL, Bstruc = NULL, Cstruc = NULL, Dstruc = NULL, maxit = 500, ctol = 1e-4, parallel = FALSE, cl = NULL, output = c("best", "all"), verbose = FALSE) ```

## Arguments

 `X` For 3-way Parafac2: list of length `K` where `k`-th element is `I[k]`-by-`J` matrix or three-way data array with `dim=c(I,J,K)`. For 4-way Parafac2: list of length `L` where `l`-th element is `I[l]`-by-`J`-by-`K` array or four-way data array with `dim=c(I,J,K,L)`. Missing data are allowed (see Note). `nfac` Number of factors. `nstart` Number of random starts. `const` Constraints for each mode. Vector of length 3 or 4 with entries: 0 = unconstrained (default), 1 = orthogonal, 2 = non-negative, 3 = unimodal, 4 = monotonic, 5 = periodic, 6 = smooth. Use `control` argument to adjust options for constraints 3-6. Note: constraints 2-4 cannot be applied on Mode A. `control` List of parameters controlling options for constraints 3-6. This is passed to `const.control`, which describes the available options. `Gfixed` Fixed Mode A crossproducts (`crossprod(Gfixed)=Phi`). Only used to fit model with fixed Phi matrix. `Bfixed` Fixed Mode B weights. Only used to fit model with fixed Mode B weights. `Cfixed` Fixed Mode C weights. Only used to fit model with fixed Mode C weights. `Dfixed` Fixed Mode D weights. Only used to fit model with fixed Mode D weights. `Gstart` Starting Mode A crossproduct matrix for ALS algorithm (`crossprod(Gstart)=Phi`). Default uses random weights. `Bstart` Starting Mode B weights for ALS algorithm. Default uses random weights. `Cstart` Starting Mode C weights for ALS algorithm. Default uses random weights. `Dstart` Starting Mode D weights for ALS algorithm. Default uses random weights. `Gstruc` Structure constraints for Mode A crossproduct matrix (`crossprod(Gstruc) = Phi structure`). Default uses unstructured crossproducts. `Bstruc` Structure constraints for Mode B weights. Default uses unstructured weights. `Cstruc` Structure constraints for Mode C weights. Default uses unstructured weights. `Dstruc` Structure constraints for Mode D weights. Default uses unstructured weights. `maxit` Maximum number of iterations. `ctol` Convergence tolerance. `parallel` Logical indicating if `parLapply` should be used. See Examples. `cl` Cluster created by `makeCluster`. Only used when `parallel=TRUE`. `output` Output the best solution (default) or output all `nstart` solutions. `verbose` Logical indicating if extra information on progress should be reported. Ignored if `parallel=TRUE`.

## Value

If `output="best"`, returns an object of class `"parafac2"` with the following elements:

 `A` List of Mode A weight matrices. `B` Mode B weight matrix. `C` Mode C weight matrix. `D` Mode D weight matrix. `Phi` Mode A crossproduct matrix. `SSE` Sum of Squared Errors. `Rsq` R-squared value. `GCV` Generalized Cross-Validation. `edf` Effective degrees of freedom. `iter` Number of iterations. `cflag` Convergence flag. `const` See argument `const`. `control` See argument `control`. `fixed` Logical vector indicating whether 'fixed' weights were used for each mode. `struc` Logical vector indicating whether 'struc' constraints were used for each mode.

Otherwise returns a list of length `nstart` where each element is an object of class `"parafac2"`.

## Warnings

The ALS algorithm can perform poorly if the number of factors `nfac` is set too large.

Non-negativity constraints can be sensitive to local optima.

Non-negativity constraints can result in slower performance.

Structure constraints for override constraints in `const` input.

## Note

Default use is 10 random strarts (`nstart=10`) with 500 maximum iterations of the ALS algorithm for each start (`maxit=500`) using a convergence tolerance of 1e-4 (`ctol=1e-4`). The algorithm is determined to have converged once the change in R^2 is less than or equal to `ctol`.

Output `cflag` gives convergence information: `cflag=0` if ALS algorithm converged normally, `cflag=1` if maximum iteration limit was reached before convergence, and `cflag=2` if ALS algorithm terminated abnormally due to a problem with the constraints.

Constraints 3 (unimodality) and 4 (monotonicity) are implemented using I-splines, and constraints 5 (periodicity) and 6 (smoothness) are implemented using M-splines (see Ramsay, 1988).

Missing data should be specified as `NA` values in the input `X`. The missing data are randomly initialized and then iteratively imputed as a part of the ALS algorithm.

## Author(s)

Nathaniel E. Helwig <[email protected]>

## References

Bro, R., & De Jong, S. (1997). A fast non-negativity-constrained least squares algorithm. Journal of Chemometrics, 11, 393-401.

Harshman, R. A. (1972). PARAFAC2: Mathematical and technical notes. UCLA Working Papers in Phonetics, 22, 30-44.

Helwig, N. E. (2013). The special sign indeterminacy of the direct-fitting Parafac2 model: Some implications, cautions, and recommendations, for Simultaneous Component Analysis. Psychometrika, 78, 725-739.

Helwig, N. E. (2017). Estimating latent trends in multivariate longitudinal data via Parafac2 with functional and structural constraints. Biometrical Journal, 59(4), 783-803.

Kiers, H. A. L., ten Berge, J. M. F., & Bro, R. (1999). PARAFAC2-part I: A direct-fitting algorithm for the PARAFAC2 model. Journal of Chemometrics, 13, 275-294.

Ramsay, J. O. (1988). Monotone regression splines in action. Statistical Science, 3, 425-441.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123``` ```########## 3-way example ########## # create random data list with Parafac2 structure set.seed(3) mydim <- c(NA,10,20) nf <- 2 nk <- sample(c(50,100,200),mydim[3],replace=TRUE) Gmat <- matrix(rnorm(nf^2),nf,nf) Bmat <- matrix(runif(mydim[2]*nf),mydim[2],nf) Cmat <- matrix(runif(mydim[3]*nf),mydim[3],nf) Xmat <- Emat <- Hmat <- vector("list",mydim[3]) for(k in 1:mydim[3]){ Hmat[[k]] <- svd(matrix(rnorm(nk[k]*nf),nk[k],nf),nv=0)\$u Xmat[[k]] <- tcrossprod(Hmat[[k]]%*%Gmat%*%diag(Cmat[k,]),Bmat) Emat[[k]] <- matrix(rnorm(nk[k]*mydim[2]),nk[k],mydim[2]) } Emat <- nscale(Emat,0,sumsq(Xmat)) # SNR=1 X <- mapply("+",Xmat,Emat) # fit Parafac2 model (unconstrained) pfac <- parafac2(X,nfac=nf,nstart=1) pfac # check solution Xhat <- fitted(pfac) sse <- sumsq(mapply("-",Xmat,Xhat)) sse/(sum(nk)*mydim[2]) crossprod(pfac\$A[[1]]) crossprod(pfac\$A[[2]]) pfac\$Phi # reorder and resign factors pfac\$B[1:4,] pfac <- reorder(pfac, 2:1) pfac\$B[1:4,] pfac <- resign(pfac, mode="B") pfac\$B[1:4,] Xhat <- fitted(pfac) sse <- sumsq(mapply("-",Xmat,Xhat)) sse/(sum(nk)*mydim[2]) # rescale factors colSums(pfac\$B^2) colSums(pfac\$C^2) pfac <- rescale(pfac, mode="C", absorb="B") colSums(pfac\$B^2) colSums(pfac\$C^2) Xhat <- fitted(pfac) sse <- sumsq(mapply("-",Xmat,Xhat)) sse/(sum(nk)*mydim[2]) ########## 4-way example ########## # create random data list with Parafac2 structure set.seed(4) mydim <- c(NA,10,20,5) nf <- 3 nk <- sample(c(50,100,200),mydim[4],replace=TRUE) Gmat <- matrix(rnorm(nf^2),nf,nf) Bmat <- scale(matrix(rnorm(mydim[2]*nf),mydim[2],nf), center=FALSE) cseq <- seq(-3, 3, length=mydim[3]) Cmat <- cbind(pnorm(cseq), pgamma(cseq+3.1, shape=1, rate=1)*(3/4), pt(cseq-2, df=4)*2) Dmat <- scale(matrix(runif(mydim[4]*nf)*2,mydim[4],nf), center=FALSE) Xmat <- Emat <- Hmat <- vector("list",mydim[4]) for(k in 1:mydim[4]){ aseq <- seq(-3,3,length=nk[k]) Hmat[[k]] <- svd(cbind(sin(aseq), sin(abs(aseq)), exp(-aseq^2)),nv=0)\$u Xmat[[k]] <- array(tcrossprod(Hmat[[k]]%*%Gmat%*%diag(Dmat[k,]), krprod(Cmat,Bmat)),dim=c(nk[k],mydim[2],mydim[3])) Emat[[k]] <- array(rnorm(nk[k]*mydim[2]*mydim[3]),dim=c(nk[k],mydim[2],mydim[3])) } Emat <- nscale(Emat,0,sumsq(Xmat)) # SNR=1 X <- mapply("+",Xmat,Emat) # fit Parafac model (smooth A, unconstrained B, monotonic C, non-negative D) pfac <- parafac2(X,nfac=nf,nstart=1,const=c(6,0,4,2)) pfac # check solution Xhat <- fitted(pfac) sse <- sumsq(mapply("-",Xmat,Xhat)) sse/(sum(nk)*mydim[2]*mydim[3]) crossprod(pfac\$A[[1]]) crossprod(pfac\$A[[2]]) pfac\$Phi ## Not run: ########## parallel computation ########## # create random data list with Parafac2 structure set.seed(3) mydim <- c(NA,10,20) nf <- 2 nk <- sample(c(50,100,200),mydim[3],replace=TRUE) Gmat <- matrix(rnorm(nf^2),nf,nf) Bmat <- matrix(runif(mydim[2]*nf),mydim[2],nf) Cmat <- matrix(runif(mydim[3]*nf),mydim[3],nf) Xmat <- Emat <- Hmat <- vector("list",mydim[3]) for(k in 1:mydim[3]){ Hmat[[k]] <- svd(matrix(rnorm(nk[k]*nf),nk[k],nf),nv=0)\$u Xmat[[k]] <- tcrossprod(Hmat[[k]]%*%Gmat%*%diag(Cmat[k,]),Bmat) Emat[[k]] <- matrix(rnorm(nk[k]*mydim[2]),nk[k],mydim[2]) } Emat <- nscale(Emat,0,sumsq(Xmat)) # SNR=1 X <- mapply("+",Xmat,Emat) # fit Parafac2 model (10 random starts -- sequential computation) set.seed(1) system.time({pfac <- parafac2(X,nfac=nf)}) pfac # fit Parafac2 model (10 random starts -- parallel computation) set.seed(1) cl <- makeCluster(detectCores()) ce <- clusterEvalQ(cl,library(multiway)) system.time({pfac <- parafac2(X,nfac=nf,parallel=TRUE,cl=cl)}) pfac stopCluster(cl) ## End(Not run) ```

multiway documentation built on Nov. 17, 2017, 7:12 a.m.