icafast: ICA via FastICA Algorithm

Description

Computes ICA decomposition using Hyvarinen's (1999) FastICA algorithm with various options.

Usage

icafast(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), 
        alg = "par", fun = "logcosh", alpha = 1)

Value

S: Matrix of source signal estimates (S = Y %*% R).
M: Estimated mixing matrix.
W: Estimated unmixing matrix (W = crossprod(R, Q)).
Y: Whitened data matrix.
Q: Whitening matrix.
R: Orthogonal rotation matrix.
vafs: Variance-accounted-for by each component.
iter: Number of algorithm iterations.
alg: Algorithm used (same as input).
fun: Contrast function (same as input).
alpha: Tuning parameter (same as input).
converged: Logical indicating if algorithm converged.

Arguments

X: Data matrix with n rows (samples) and p columns (variables).
nc: Number of components to extract.
center: If TRUE, columns of X are mean-centered before ICA decomposition.
maxit: Maximum number of algorithm iterations to allow.
tol: Convergence tolerance.
Rmat: Initial estimate of the nc-by-nc orthogonal rotation matrix.
alg: Algorithm to use: alg="par" to estimate all nc components in parallel (default) or alg="def" for deflation estimation (i.e., projection pursuit).
fun: Contrast function to use for negentropy approximation: fun="logcosh" for log of hyperbolic cosine, fun="exp" for Gaussian kernel, and fun="kur" for kurtosis.
alpha: Tuning parameter for "logcosh" contrast function (1 <= alpha <= 2).

Author

Nathaniel E. Helwig <helwig@umn.edu>

Details

ICA Model The ICA model can be written as X = tcrossprod(S, M) + E, where S contains the source signals, M is the mixing matrix, and E contains the noise signals. Columns of X are assumed to have zero mean. The goal is to find the unmixing matrix W such that columns of S = tcrossprod(X, W) are independent as possible.

Whitening Without loss of generality, we can write M = P %*% R where P is a tall matrix and R is an orthogonal rotation matrix. Letting Q denote the pseudoinverse of P, we can whiten the data using Y = tcrossprod(X, Q). The goal is to find the orthongal rotation matrix R such that the source signal estimates S = Y %*% R are as independent as possible. Note that W = crossprod(R, Q).

FastICA The FastICA algorithm finds the orthogonal rotation matrix R that (approximately) maximizes the negentropy of the estimated source signals. Negentropy is approximated using $J (s) = [E (G (s)) - E (G (z))]^{2}$ where E denotes the expectation, G is the contrast function, and z is a standard normal variable. See Hyvarinen (1999) or Helwig and Hong (2013) for specifics of fixed-point algorithm.

References

Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. tools:::Rd_expr_doi("https://doi.org/10.1016/j.jneumeth.2012.12.009")

Hyvarinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626-634. tools:::Rd_expr_doi("10.1109/72.761722")

Examples

Run this code

##########   EXAMPLE 1   ##########

# generate noiseless data (p == r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via FastICA with 2 components
imod <- icafast(Xmat, nc = 2)
acy(Bmat, imod$M)
cor(Amat, imod$S)



##########   EXAMPLE 2   ##########

# generate noiseless data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2)
Xmat <- tcrossprod(Amat, Bmat)

# ICA via FastICA with 2 components
imod <- icafast(Xmat, nc = 2)
cor(Amat, imod$S)



##########   EXAMPLE 3   ##########

# generate noisy data (p != r)
set.seed(123)
nobs <- 1000
Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs))
Bmat <- matrix(2 * runif(200), 100, 2)
Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100)
Xmat <- tcrossprod(Amat,Bmat) + Emat

# ICA via FastICA with 2 components
imod <- icafast(Xmat, nc = 2)
cor(Amat, imod$S)

Run the code above in your browser using DataLab

Get 50% off unlimited learning