icaimax {ica}

ICA via Infomax Algorithm

Description

Computes ICA decomposition using Bell and Sejnowski's (1995) Information-Maximization (Infomax) approach with various options.

Usage

icaimax(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), 
        alg = "newton", fun = "tanh", signs = rep(1, nc), signswitch = TRUE, 
        rate = 1, rateanneal = NULL)

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="newton" for Newton iteration, and alg="gradient" for gradient descent.

fun

Nonlinear (squashing) function to use for algorithm: fun="tanh" for hyperbolic tangent, fun="log" for logistic, and fun="ext" for extended Infomax.

signs

Vector of length nc such that signs[j] = 1 if j-th component is super-Gaussian and signs[j] = -1 if j-th component is sub-Gaussian. Only used if fun="ext". Ignored if signswitch=TRUE.

signswitch

If TRUE, the signs vector is automatically determined from the data; otherwise a confirmatory ICA decomposition is calculated using input signs vector. Only used if fun="ext".

rate

Learing rate for gradient descent algorithm. Ignored if alg="newton".

rateanneal

Annealing angle and proportion for gradient descent learing rate (see Details). Ignored if alg="newton".

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

Infomax The Infomax approach finds the orthogonal rotation matrix R that (approximately) maximizes the joint entropy of a nonlinear function of the estimated source signals. See Bell and Sejnowski (1995) and Helwig and Hong (2013) for specifics of algorithms.

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

signs

Component signs (same as input).

rate

Learning rate (same as input).

converged

Logical indicating if algorithm converged.

Author(s)

Nathaniel E. Helwig <[e-mail]>

References

Bell, A.J. & Sejnowski, T.J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6), 1129-1159. doi:10.1162/neco.1995.7.6.1129

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. doi:10.1016/j.jneumeth.2012.12.009

See Also

icafast for FastICA

icajade for ICA via JADE

Examples

##########   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 Infomax with 2 components
imod <- icaimax(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 Infomax with 2 components
imod <- icaimax(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 Infomax with 2 components
imod <- icaimax(Xmat, nc = 2)
cor(Amat, imod$S)