ICA via distance covariance, tests of mutual independence, and other ICA functions

Description

Functions related to multivariate measures of independence and ICA:
-estimate independent components by minimizing distance covariance;
-conduct a test of mutual independence based on distance covariance;
-estimate independent components via infomax (a popular method but generally performs poorer than steadyICA or ProDenICA but is useful for comparisons);
-order independent components by skewness;
-match independent components from multiple estimates;
-other functions useful in ICA.

Details

Author(s)

Benjamin B. Risk and Nicholas A. James and David S. Matteson.
Maintainer: Benjamin Risk <bbr28@cornell.edu>

References

Bernaards, C. & Jennrich, R. (2005) Gradient projection algorithms and software for arbitrary rotation criteria in factor analysis. Educational and Psychological Measurement 65, 676-696

Matteson, D. S. & Tsay, R. Independent component analysis via U-Statistics. <http://www.stat.cornell.edu/~matteson/#ICA>

Szekely, G., Rizzo, M. & Bakirov, N. Measuring and testing dependence by correlation of distances. (2007) The Annals of Statistics, 35, 2769-2794.

Tichavsky, P. & Koldovsky, Z. Optimal pairing of signal components separated by blind techniques. (2004) Signal Processing Letters 11, 119-122.

See Also

fastICA ProDenICA::ProDenICA

Examples

#see steadyICA

Convert an orthogonal matrix to its angular parameterization.

Description

Convert a d x d orthogonal matrix to a sequence of d*(d-1)/2 Givens rotations.

Usage

W2theta(W)

Arguments

Details

A d x d orthogonal matrix can be decomposed into a series of d*(d-1)/2 Givens rotation matrices, where each matrix is parameterized by a single angle.

Value

A vector of length d*(d-1)/2 comprised of the angles.

Author(s)

David S. Matteson

References

Golub, G. & Van Loan, C. 1996. Matrix computations. Johns Hopkins University Press.

See Also

theta2W

Examples

theta = c(pi/6,pi/4,pi/2)

(W = theta2W(theta))

#Recover theta:
W2theta(W)

Complete Measure of Mutual Multivariate Independence

Description

Calculates a complete empirical measure of mutual multivariate independence. Makes use of the utils::combn function.

Usage

compInd(S,group=1:ncol(S),alpha=1)

Arguments

Value

Returns a scalar equal to the empirical multivariate distance between the observed samples, and their grouped counterpart.

Note

Suppose that the each component belongs to exactly one of C groups. This method makes use of the utils::combn and combinat::permn functions. As a result it will be both computationally and memory intensive, even for small to moderate n and small C.

Author(s)

Nicholas James

References

Chasalow, Scott (2012) combinat: Combinatorics Utilities <http://CRAN.R-project.org/package=combinat

See Also

dcovustat, energy::dcov

Examples

library(steadyICA)
library(combinat)
set.seed(100)
S = matrix(rnorm(40),ncol=4)
group = c(1,2,3,3)
compInd(S,group,1)

ICA via distance covariance for 2 components

Description

This algorithm finds the rotation which minimizes the distance covariance between two orthogonal components via the angular parameterization of a 2x2 orthogonal matrix with the function stats::optimize. The results will be (approximately) equivalent to steadyICA but this function is much faster (but does not extend to higher dimensions).

Usage

dcovICA(Z, theta.0 = 0)

Arguments

Value

Author(s)

David Matteson and Benjamin Risk

References

Matteson, D. S. & Tsay, R. Independent component analysis via U-Statistics. <http://www.stat.cornell.edu/~matteson/#ICA>

See Also

steadyICA, optimize

Examples

library(JADE)
library(ProDenICA)
set.seed(123)
simS = cbind(rjordan(letter='j',n=1024),rjordan(letter='m',n=1024))
simM = mixmat(p=2)
xData = simS%*%simM
xWhitened = whitener(xData)

#Define true unmixing matrix as true M multiplied by the estimated whitener:
#Call this the target matrix:
W.true <- solve(simM%*%xWhitened$whitener) 


a=Sys.time()
est.dCovICA = dcovICA(Z = xWhitened$Z,theta.0=0)
Sys.time()-a

#See the example with steadyICA for an explanation
#of the parameterization used in amari.error:
amari.error(t(est.dCovICA$W),W.true)

##NOTE: also try theta.0 = pi/4 since there may be local minima
  ## Not run: est.dcovICA = dcovICA(Z = xWhitened$Z,theta.0=pi/4)
  amari.error(t(est.dcovICA$W),W.true)
## End(Not run)

a=Sys.time()
est.steadyICA = steadyICA(X=xWhitened$Z,verbose=TRUE)
Sys.time()-a
amari.error(t(est.steadyICA$W),W.true)
##theta parameterization with optimize is much faster

Calculate distance covariance via U-statistics

Description

Calculates the square of the U-statistic formulation of distance covariance. This is faster than the function 'dcov' in the R package 'energy' and requires less memory. Note that negative values are possible in this version.

Usage

dcovustat(x,y,alpha=1)

Arguments

Value

Returns the distance covariance U-statistic.

Note

The value returned by dcovustat is equal to the square of the value returned by energy::dcov in the limit.

In dcovustat, a vector of length n is stored; in energy::dcov, an n x n matrix is stored. Thus, dcovustat requires far less memory and works for very large datasets.

Even though dcovustat converges to the square of the distance covariance of the random variables x and y, it can be negative.

Author(s)

David Matteson

References

Matteson, D. S. & Tsay, R. Independent component analysis via U-Statistics. <http://www.stat.cornell.edu/~matteson/#ICA>

Szekely, G., Rizzo, M. & Bakirov, N. Measuring and testing dependence by correlation of distances. (2007) The Annals of Statistics, 35, 2769-2794.

See Also

multidcov, energy::dcov

Examples

x = rnorm(5000)
y = rbinom(5000,1,0.5)
y = y - 1*(y==0)
z = y*exp(-x) #some non-linear dependence

dcovustat(x[1:1000],y[1:1000]) #close to zero

a = Sys.time()
dcovustat(x[1:1000],z[1:1000]) #greater than zero
a = Sys.time() - a

#measures of linear dependence close to zero:
cov(x,z)
cor(rank(x),rank(z))


## Not run: 
#dcovustat differs from energy::dcov but are equal in the limit
library(energy)
b = Sys.time()
(dcov(x[1:1000],z[1:1000]))^2
b = Sys.time() - b
as.double(b)/as.double(a) #dcovustat is much faster

## energy::dcov and dcovustat become approximately equal as n increases:
c = Sys.time()
dcovustat(x,z)
c = difftime(Sys.time(), c, sec)
d = Sys.time()
(dcov(x,z)^2)
d = difftime(Sys.time(), d, sec)
as.double(d)/as.double(c) 

## End(Not run)

match mixing matrices or ICs and calculate their Frobenius distance

Description

The ICA model is only identifiable up to signed permutations of the ICs. This function provides a similarity measure between two mixing matrices for the model X = S M + E, where X is n x p, S is n x d, and M is d x p. The input is either two mixing matrices M1 and M2 or two matrices of independent components S1 and S2. For M1 and M2, frobICA() finds the signed row permutation of M2 that minimizes the Frobenius norm between M1 and M2 using the Hungarian method. For S1 and S2, frobICA() finds the signed column permutation of S2 that minimizes the Frobenius norm between S1 and S2. This function allows the mixing matrices (or independent components) to have differing numbers of rows (respectively, columns) such that the similarity measure is defined by the matching rows (resp., columns), and the non-matching rows (resp., columns) are discarded.

Usage

frobICA(M1 = NULL, M2 = NULL, S1 = NULL, S2 = NULL, standardize = FALSE)

Arguments

Details

frobICA(M1,M2) = 0 if there exists a signed permutation of the rows of M2 such that M1 = P%*%M2, where P is a d x q signed permutation matrix, i.e., composed of 0, 1, and -1, with d <= q; the function also allows d > q, in which case frobICA(M1,M2) = 0 if there exists a P such that P%*% M1 = M2. Unlike other ICA performance measures, this function can accomodate non-square mixing matrices.

Value

returns the Frobenius norm divided by p*min(d,q) (or n*min(d,q)) of the matched mixing matrices (resp., matched independent components).

Note

If standardize=TRUE, then scales the rows of M1 and M2 to have unit norm or the columns of S1 and S2 to have zero mean and sample variance equal to one. The user can supply either M1 and M2 or S1 and S2 but not both.

Author(s)

Benjamin Risk

References

Kuhn, H. The Hungarian Method for the assignment problem Naval Research Logistics Quarterly, 1955, 2, 83 - 97

Risk, B.B., D.S. Matteson, D. Ruppert, A. Eloyan, B.S. Caffo. In review, 2013. Evaluating ICA methods with an application to resting state fMRI.

See Also

JADE::MD clue::solve_LSAP matchICA

Examples

mat1 <- matrix(rnorm(4*6),nrow=4)
perm <- matrix(c(-1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1),4,4)
mat2 <- perm%*%mat1
sqrt(sum((mat1-mat2)^2))
frobICA(M1=mat1,M2=mat2)

#Another example showing invariance to permutations:
covMat <- t(mat1)%*%mat1
mvsample <- matrix(rnorm(400),100,4)%*%mat1
frobICA(M1=cov(mvsample),M2=covMat)
frobICA(M1=cov(mvsample),M2=covMat[sample(1:6),])

#Example using independent components:
nObs=300
simS<-cbind(rgamma(nObs, shape = 1, scale = 2),
            rgamma(nObs, shape = 3, scale = 2),
            rgamma(nObs, shape = 3, scale = 2),
            rgamma(nObs, shape = 9, scale = 0.5))
              
#not necessary in this example, but this should be done when used with ICA:            
simS <- apply(simS,2,scale) 
frobICA(S1=simS,S2=simS%*%perm)
## Not run: 
#returns an error if S1 and S2 are not explicitly defined:
frobICA(simS,simS%*%perm)            

## End(Not run)            

Symmetric multivariate distance covariance for grouped components

Description

Calculate either the symmetric or asymmetric multivariate distance covariance statistic for a given grouping of the components.

Usage

gmultidcov(S,group=1:ncol(S),alpha=1,symmetric=TRUE)

Arguments

Details

Suppose that the groups are numbered 1,2,...,C and that group is a vector indicating group membership for each component. If symmetric==TRUE, calculates: sum_i=1^C dcovustat(S[,group==i],S[,group!=i]) If symmetric==FALSE, calculates: sum_i=1^C-1 dcovustat(S[,group==i],S[,group>i])

Value

Returns a scalar equal to the multivariate distance covariance statistic for grouped components of S.

Author(s)

Nicholas James

See Also

dcovustat, energy::dcov

Examples

library(steadyICA)
S = matrix(rnorm(300),ncol=3)
group = c(1,2,2)
gmultidcov(S,group,TRUE) # close to zero
gmultidcov(S,group,FALSE) # sill close to zero

Sigma = matrix(c(1,0.7,0,0.7,1,-0.2,0,-0.2,1),ncol=3)
X = MASS::mvrnorm(100,rep(0,3),Sigma)
gmultidcov(X,group,TRUE) # further from zero
gmultidcov(X,group,FALSE) # further from zero

Estimates independent components via infomax

Description

Estimate independent components using the infomax criteria, which is equivalent to maximum likelihood using the logistic density, exp(-S)/(1+exp(-S))^2.

Usage

infomaxICA(X, n.comp, W.list = NULL, whiten = FALSE, maxit = 500, eps = 1e-08, 
alpha.eps = 1e-08, verbose = FALSE, restarts=0) 

Arguments

Details

This is an R version of ICA using the infomax criteria that provides an alternative to Matlab code (ftp://ftp.cnl.salk.edu/pub/tony/sep96.public), but with a few modifications. First, we use the full data (the so-called offline algorithm) in each iteration rather than an online algorithm with batches. Secondly, we use an adaptive method to choose the step size (based upon Bernaards and Jennrich 2005), which speeds up convergence. We also omitted the bias term (intercept) included in the original formulation because we centered our data.

Value

Note

In contrast to most other ICA methods, W is not contrained to be orthogonal.

Author(s)

Benjamin Risk

References

Bell, A. & Sejnowski, T. An information-maximization approach to blind separation and blind deconvolution Neural computation, Neural computation, 1995, 7, 1129-1159.

Bernaards, C. A. and Jennrich, R. I. (2005) Gradient Projection Algorithms and Software for Arbitrary Rotation Criteria in Factor Analysis, Educational and Psychological Measurement 65, 676-696. <http://www.stat.ucla.edu/research/gpa>

Examples

## Example when p > d. The MD function and amari measures 
# are not defined for M. We can compare the 
# "true W inverse", which is the mixing matrix multiplied 
# by the whitening matrix; alternatively, we can use
# multidcov::frobICA. These two approaches are
# demonstrated below:

set.seed(999)
nObs <- 1024
nComp <- 3

# simulate from gamma distributions with
# varying amounts of skewness:
simS<-cbind(rgamma(nObs, shape = 1, scale = 2),
            rgamma(nObs, shape = 3, scale = 2),
            rgamma(nObs, shape = 9, scale = 0.5))

#standardize by expected value and variance:
simS[,1] = (simS[,1] - 1*2)/sqrt(1*2^2)
simS[,2] = (simS[,2] - 3*2)/sqrt(3*2^2)
simS[,3] = (simS[,3] - 9*0.5)/sqrt(9*0.5^2)

# slightly revised 'mixmat' function (from ProDenICA)
# for p>=d: uses fastICA and ProDenICA parameterization:
myMixmat <- function (p = 2, d = NULL) {
  if(is.null(d)) d = p
  a <- matrix(rnorm(d * p), d, p)
  sa <- La.svd(a)
  dL <- sort(runif(d) + 1)
  mat <- sa$u%*%(sa$vt*dL)
  attr(mat, "condition") <- dL[d]/dL[1]
  mat
}

simM <- myMixmat(p = 6, d = nComp)
xData <- simS%*%simM
xWhitened <- whitener(xData, n.comp = nComp)
#Define a 'true' W (uses the estimated whitening matrix):
W.true <- solve(simM%*%xWhitened$whitener)

estInfomax <- infomaxICA(X = xData, n.comp = nComp, whiten = TRUE, verbose = TRUE)

frobICA(estInfomax$M,simM)
library(JADE)
MD(t(estInfomax$W),t(solve(W.true)))
amari.error(t(estInfomax$W),t(solve(W.true)))

match independent components using the Hungarian method

Description

The ICA model is only identifiable up to signed permutations of the ICs. This function finds the signed permutation of a matrix S such that ||S%*%P - template|| is minimized. Optionally also matches the mixing matrix M.

Usage

matchICA(S, template, M = NULL)

Arguments

Value

Returns the signed permutation of S that is matched to the template. If the optional argument M is provided, returns a list with the permuted S and M matrices.

Author(s)

Benjamin Risk

References

Kuhn, H. The Hungarian Method for the assignment problem Naval Research Logistics Quarterly, 1955, 2, 83 - 97

Risk, B.B., D.S. Matteson, D. Ruppert, A. Eloyan, B.S. Caffo. In review, 2013. Evaluating ICA methods with an application to resting state fMRI.

See Also

frobICA clue::solve_LSAP

Examples

set.seed(999)
nObs <- 1024
nComp <- 3
# simulate from some gamma distributions:
simS<-cbind(rgamma(nObs, shape = 1, scale = 2),
            rgamma(nObs, shape = 3, scale = 2),
            rgamma(nObs, shape = 9, scale = 0.5))

simM <- matrix(rnorm(9),3)
pMat <- matrix(c(0,-1,0,1,0,0,0,0,-1),3)
permS <- simS%*%pMat
permM <- t(pMat)%*%simM

matchedS <- matchICA(S = permS, template = simS, M = permM)
sum(abs(matchedS$S - simS))
sum(abs(simM - matchedS$M))

Symmetric multivariate distance covariance

Description

Calculate either the symmetric or asymmetric multivariate distance covariance statistic.

Usage

multidcov(S,symmetric=TRUE,alpha=1)

Arguments

Details

If symmetric==TRUE, calculates: sum_i=1^d dcovustat(S[,i],S[,-i]) If symmetric==FALSE, calculates: sum_i=1^d-1 dcovustat(S[,i],S[,(i+1):d])

Value

returns a scalar equal to the multivariate distance covariance statistic for the columns of S

Author(s)

David Matteson

See Also

dcovustat, energy::dcov

Examples

nObs <- 1024
nComp <- 3

simM <- matrix(rnorm(nComp*nComp),nComp)

# simulate some data:
simS<-cbind(rgamma(nObs, shape = 1, scale = 2),
            rgamma(nObs, shape = 3, scale = 2),
            rgamma(nObs, shape = 9, scale = 0.5))


simS <- scale(simS) #Standardize variance for identifiability

#mix the sources:
xData <- simS %*% simM

multidcov(simS) #close to zero
multidcov(whitener(xData)$Z) #should be larger than simS
multidcov(xData) #greater than zero

Permutation test for mutual independence.

Description

Calculates an approximate p-values based upon a permutation test for mutual independence.

Usage

permTest(S, group=1:ncol(S), R=199, FUN=c('gmultidcov','compInd'), ...)

Arguments

Details

Suppose that the groups are numbered 1,2,...,C and that group is a vector indicating group membership for each component. If symmetric==TRUE, calculates: sum_i=1^C dcovustat(S[,group==i],S[,group!=i]) If symmetric==FALSE, calculates: sum_i=1^C-1 dcovustat(S[,group==i],S[,group>i])

If no additional arguments are supplied for FUN then the default values are used. In the case of gmultidcov, values for alpha and symmetric can be supplied. While for compInd only the value of alpha is needed.

Value

Returns an approximate p-values based upon a permutation test.

Author(s)

Nicholas James

See Also

dcovustat, energy::dcov

force ICs to have positive skewness and order by skewness

Description

The ICA model is only identifiable up to signed permutations. This function provides a canonical ordering for ICA that is useful for fMRI or studies where signals are skewed. Multiplies columns of S that are left-skewed by -1 to force right skewness. Optionally orders the columns by descending skewness.

Usage

rightskew(S, M = NULL, order.skew = TRUE)

Arguments

Value

Returns the matrix S such that all columns have positive skewness. If optional argument M is supplied, returns a list with the new S and corresponding M.

Author(s)

Benjamin Risk

References

Eloyan, A. & Ghosh, S. A Semiparametric Approach to Source Separation using Independent Component Analysis Computational Statistics and Data Analysis, 2013, 58, 383 - 396.

Examples

nObs = 1024
simS<-cbind(rgamma(nObs, shape = 1, scale = 2),
            rgamma(nObs, shape = 9, scale = 0.5),
            -1*rgamma(nObs, shape = 3, scale = 2))

apply(simS,2,function(x){ 
  (sum((x - mean(x))^3)/length(x))/(sum((x - mean(x))^2)/length(x))^(3/2)})

canonicalS <- rightskew(simS)

apply(canonicalS,2,function(x){
 (sum((x - mean(x))^3)/length(x))/(sum((x - mean(x))^2)/length(x))^(3/2)})

Estimate independent components by minimizing distance covariance

Description

The model is: X = S M + E, where X is n x p and has mean zero, S is n x d, M is d x p, and E is measurement error. For whitened data, we have Z = S t(W), where W is orthogonal. We find the matrix M such that S minimizes the distance covariance dependency measure.

Usage

steadyICA(X, n.comp = ncol(X), w.init = NULL, PIT = FALSE, bw = 'SJ', adjust = 1,
whiten = FALSE, irlba = FALSE, symmetric = FALSE, eps = 1e-08, alpha.eps = 1e-08,
maxit = 100, method = c('Cpp','R'), verbose = FALSE)

Arguments

Value

Author(s)

Benjamin Risk

References

Matteson, D. S. & Tsay, R. Independent component analysis via U-Statistics. <http://www.stat.cornell.edu/~matteson/#ICA>

See Also

multidcov

Examples

set.seed(999)
nObs <- 1024
nComp <- 3
# simulate from some gamma distributions:
simS<-cbind(rgamma(nObs, shape = 1, scale = 2),
            rgamma(nObs, shape = 3, scale = 2),
            rgamma(nObs, shape = 9, scale = 0.5))

#standardize by expected value and variance:
simS[,1] = (simS[,1] - 1*2)/sqrt(1*2^2)
simS[,2] = (simS[,2] - 3*2)/sqrt(3*2^2)
simS[,3] = (simS[,3] - 9*0.5)/sqrt(9*0.5^2)

# slightly revised 'mixmat' function (from ProDenICA)
# for p>=d: uses fastICA and ProDenICA parameterization:
myMixmat <- function (p = 2, d = NULL) {
  if(is.null(d)) d = p
  a <- matrix(rnorm(d * p), d, p)
  sa <- La.svd(a)
  dL <- sort(runif(d) + 1)
  mat <- sa$u%*%(sa$vt*dL)
  attr(mat, "condition") <- dL[d]/dL[1]
  mat
}

simM <- myMixmat(p = 6, d = nComp)
xData <- simS%*%simM
xWhitened <- whitener(xData, n.comp = nComp)

#estimate mixing matrix:
est.steadyICA.v1 = steadyICA(X = xData,whiten=TRUE,n.comp=nComp,verbose = TRUE)

#Define the 'true' W:
W.true <- solve(simM%*%xWhitened$whitener)

frobICA(M1=est.steadyICA.v1$M,M2=simM)
frobICA(S1=est.steadyICA.v1$S,S2=simS)

## Not run: 
#now initiate from target:
est.steadyICA.v2 = steadyICA(X = xData, w.init= W.true, n.comp = nComp, whiten=TRUE, verbose=TRUE)

#estimate using PIT steadyICA such that dimension reduction is via ICA:
est.steadyICA.v3 = steadyICA(X = xData, w.init=ginv(est.steadyICA.v2$M),
PIT=TRUE, n.comp = nComp, whiten=FALSE, verbose=TRUE)

frobICA(M1=est.steadyICA.v2$M,M2=simM)
frobICA(M1=est.steadyICA.v3$M,M2=simM) 
frobICA(S1=est.steadyICA.v2$S,S2=simS)

#tends to be lower than PCA-based (i.e., whitening) methods:
frobICA(S1=est.steadyICA.v3$S,S2=simS) 

# JADE uses a different parameterization and different notation.
# Using our parameterization and notation, the arguments for 
# JADE::amari.error correspond to:
amari.error(t(W.hat), W.true)

library(JADE)

amari.error(t(est.steadyICA.v1$W), W.true) 
amari.error(t(est.steadyICA.v2$W), W.true)
##note that a square W is not estimated if PIT=TRUE and whiten=FALSE

#Compare performance to fastICA:
library(fastICA)
est.fastICA = fastICA(X = xData, n.comp = 3, tol=1e-07)
amari.error(t(est.fastICA$W), W.true)
##steadyICA usually outperforms fastICA

##Compare performance to ProDenICA:
library(ProDenICA)
est.ProDenICA = ProDenICA(x = xWhitened$Z, k = 3, maxit=40,trace=TRUE)
amari.error(t(est.ProDenICA$W), W.true)
##ProDenICA and steadyICA tend to be similar when sources
##are continuously differentiable

## End(Not run)

Convert angles to an orthogonal matrix.

Description

Convert d*(d-1)/2 angles from a sequence of Givens rotations to a d x d orthogonal matrix.

Usage

theta2W(theta)

Arguments

Value

A d x d orthogonal matrix resulting from the sequence of d*(d-1)/2 Givens rotation matrices.

Author(s)

David S. Matteson

References

Golub, G. & Van Loan, C. 1996. Matrix computations. Johns Hopkins University Press.

See Also

W2theta

Examples

#Generate orthogonal matrix:
mat <- matrix(rnorm(9),3,3)
W = svd(mat)$u

theta <- W2theta(W)

#Recovers W:
theta2W(theta)

Whitening function

Description

Subtract column means and transform columns such that the empirical covariance is equal to the identity matrix. Uses the SVD.

Usage

whitener(X, n.comp = ncol(X), center.row = FALSE, irlba = FALSE)

Arguments

Value

Note

The use of the option 'irlba = TRUE' requires the package irlba and is very useful for large p. The function irlba only calculates the first n.comp eigenvectors and is much faster than svd for p >> n.comp, for e.g., in groupICA of fMRI data.

Author(s)

Benjamin Risk

See Also

svd, irlba::irlba

Examples

simData <- cbind(rnorm(1000,1,2),rnorm(1000,-1,3),rnorm(1000,4,1))
simMVN <- simData%*%matrix(rnorm(12),3,4)
simWhiten <- whitener(simMVN,n.comp = 3)
colMeans(simWhiten$Z)
cov(simWhiten$Z)