coloredICA (version 1.0.0)

cICA: colored Independent Component Analysis

Description

This function implements the colored Independent Component Analysis (cICA) algorithm, where sources are treated as temporal stochastic processes.

Usage

cICA(Xin, M = dim(Xin)[1], Win = diag(M), tol = 1e-04, maxit = 20, nmaxit = 1, 
unmixing.estimate = "eigenvector", maxnmodels = 100)

Arguments

Xin

Data matrix with p rows (representing variables) and n columns (representing observations).

M

Number of components to be extracted.

Win

Initial guess for the unmixing matrix W. Dimensions need to be M x M.

tol

Tolerance used to establish the convergence of the algorithm.

maxit

Maximum number of iterations.

nmaxit

If the algorithm does not converge, it is run again with a new initial guess for the unmixing matrix W. This operation is done nmaxit times.

unmixing.estimate

The method used in the unmixing matrix estimation step. The two allowed choices are "eigenvector" and "newton" (see Details).

maxnmodels

Maximum number of models tested in the spectral density estimation step of the algorithm (see Details).

Value

• A list containing the following components:
• W
  Estimate of the M x M unmixing matrix in the whitened space.
• K
  pre-whitening matrix that projects data onto the first M principal components. Dimensions are M x p.
• A
  Estimate of the p x M mixing matrix.
• S
  Estimate of the M x n source matrix.
• X
  Original p x n data matrix.
• iter
  number of iterations.
• NInv
  number of times the algorithm is rerun after it does not achieve convergence.
• den
  Estimate of the spectral density of the sources. Dimensions are M x n.

Details

In the Independent Component Analysis approach, the data matrix $X$ is considered to be a linear combination of independent components, i.e. $X = AS$, where rows of $S$ contain the unobserved realizations of the independent components and $A$ is a linear mixing matrix. According to classical ICA procedures data matrix $X$ is centered and, then, whitened by projecting the data onto its principal component directions, i.e. $X \rightarrow KX = \widetilde{X}$ where $K$ is a M x p pre-whitening matrix. The cICA algorithm then estimates the unmixing matrix $W$, with $W\widetilde{X} = S$, according to the procedure described below. Then, defining $\widetilde{W}=WK$, the mixing matrix $A$ is recovered through $A=\widetilde{W}^T(\widetilde{W}\widetilde{W}^T)^{-1}$. Colored Independent Component Analysis assumes that the independent sources are temporal stochastic processes. To perform ICA, the Whittle log-likelihood is exploited. In particular the log-likelihood is written in function of the unmixing matrix $W$ and the spectral densities $f_{S_j}$ of the autocorrelated sources as follows: $$l(W,\boldsymbol{f}_{\boldsymbol{S}};\widetilde{X})=\sum_{j=1}^p\sum_{k=1}^n\left(\frac{\boldsymbol{e}_j^T W \widetilde{\boldsymbol{f}}(r_k,\widetilde{X})W^T \boldsymbol{e}_j}{f_{S_j}(r_k)}+\ln f_{S_j}(r_k)\right) + n\ln|\det(W)|.$$ Due to whitening, $W$ is orthogonal and the last term of the objective function can be dropped. The orthogonality of the unmixing matrix $W$ can be imposed in two different ways, setting the argument unmixing.estimate. In this way the estimate of the unmixing matrix $W$ can be found according two different procedures:

• as described in Shen et al. (2014). A penalty term is added to the objective function. In particular$\tau\bold{w}'_{j}C_{j}\bold{w}_{j}$, where$\bold{w}'_{j}$is the$j$th column of$W$,$C_j=\sum_{k\neq j}\bold{w}_k\bold{w}'_k$and$\tau$is a tuning parameter. The matrix$C_j$provides an orthogonality constraint in the sense that$\bold{w}'_{j}C_{j}\bold{w}_{j}=\sum_{k\neq j}$. In this way the objective function assumes a symmetric and positive-definite form and the argmin correspond to the lower eigenvalue. This choice is obtained settingunmixing.matrix = "eigenvector".
• as described in Lee et al. (2011). The orthogonality constraint is considered performing the minimization of the objective function according a Newton-Raphson method with Lagrange multiplier. This choice is obtained settingunmixing.matrix = "newton".

Independently from the choice of the technique to minimize the objective function, the cICA algorithm is based on an iterative procedure. While the Amari error is greater than tol and the number of iteration is less or equal than maxit, the two following steps are repeated:

• parametric estimation of the sources spectral densities using the Yule-Walker method, evaluatingmaxmodelsmodels.
• estimate the unmixing matrix$W$according the method selected inunmixing.estimate.

References

Lee, S., Shen, H., Truong, Y., Lewis, M., Huang, X. (2011). Independent Component Analysis Involving Autocorrelated Sources With an Application to Funcional Magnetic Resonance Imaging. Journal of the American Statistical Association, 106, 1009--1024. Shen, H., Truong, Y., Zanini, P. (2014). Independent Component Analysis for Spatial Processes on a Lattice. MOX report 38/2014, Department of Mathematics, Politecnico di Milano.

See Also

scICA

Examples

## Not run:

require(fastICA)

T=256
n1=16
n2=16
M=2

S1 = arima.sim(list(order=c(0,0,2),ma=c(1,0.25)),T)
S2 = arima.sim(list(order=c(1,0,0), ar=-0.5),T,rand.gen = function(n, ...) (runif(n)-0.5)*sqrt(3))

A = rerow(matrix(runif(M^2)-0.5,M,M))
W = solve(A)
S=rbind(S1,S2)
X = A %*% S

cica = cICA(X,tol=0.001)
## scica = scICA(X,n1=n1,n2=n2,h=0.8,tol=0.001)
fica = fastICA(t(X),2)

amari_distance(t(A),t(cica$A))
## amari_distance(t(A),t(scica$A))
amari_distance(t(A),fica$A)

Shat1=cica$S
## Shat2=scica$S
Shat3=t(fica$S)

par(mfrow=c(2,2))
plot(S[1,],type="l",lwd=2)
plot(S[2,],type="l",lwd=2)
plot(Shat1[1,],type="l",lwd=2,col="red")
plot(Shat1[2,],type="l",lwd=2,col="red")

## par(mfrow=c(2,2))
## plot(S[1,],type="l",lwd=2)
## plot(S[2,],type="l",lwd=2)
## plot(Shat2[1,],type="l",lwd=2,col="green")
## plot(Shat2[2,],type="l",lwd=2,col="green")

par(mfrow=c(2,2))
plot(S[1,],type="l",lwd=2)
plot(S[2,],type="l",lwd=2)
plot(Shat3[1,],type="l",lwd=2,col="blue")
plot(Shat3[2,],type="l",lwd=2,col="blue")

## End (Not run)