scICA: spatial colored Independent Component Analysis
In coloredICA: Implementation of Colored Independent Component Analysis and Spatial Colored Independent Component Analysis
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
Description
This function implements the spatial colored Independent Component Analysis (scICA) algorithm, where sources are treated as spatial stochastic processes on a lattice.
Usage
1
2
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).
n1
Number of rows of the lattice.
n2
Number of columns of the lattice.
nx01
Number of rows of the lattice where the spectral density is evaluate. Default value is n1.
nx02
Number of columns of the lattice where the spectral density is evaluate. Default value is n2.
h
Kernel bandwidth used for the nonparametric estimation of the sources spectral densities.
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 scICA 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}.
Spatial colored Independent Component Analysis assumes that the independent sources are spatial stochastic processes on a lattice. 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 spatial autocorrelated sources as follows:
l(W,\boldsymbol{f}_{\boldsymbol{S}};\widetilde{X})=∑_{j=1}^p∑_{k=1}^n≤ft(\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 τ\bold{w}'_{j}C_{j}\bold{w}_{j}, where \bold{w}'_{j} is the jth column of W, C_j=∑_{k\neq j}\bold{w}_k\bold{w}'_k and τ is a tuning parameter. The matrix C_j provides an orthogonality constraint in the sense that \bold{w}'_{j}C_{j}\bold{w}_{j}=∑_{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 setting unmixing.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 setting unmixing.matrix = "newton".
Independently from the choice of the technique to minimize the objective function, the scICA 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:
nonparametric estimation of the sources spectral density through a multidimensional local linear kernel estimator \widehat{m}_{LK} (see Shen et al. (2014) for further details).
estimate the unmixing matrix W according the method selected in unmixing.estimate.
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.
Note
Note that source matrix S and spectral density matrix den are n x M matrices. Every row, that should be in a n1 x n2 grid, has been vectorized in a n vector by column, with n = n1 x n2.
Author(s)
Lee, S., Shen, H., Truong, Y. and Zanini, P.
References
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.
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.
See Also
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
## Not run:
require(fastICA)
n1=20
n2=20
M=2
# Fist source
sigma1=2
S1=matrix(0,n1,n2)
for (i in 1:n1){
S1[i,]=rnorm(n2,i*2,0.2)
}
for (j in 1:n2){
S1[,j]=S1[,j]+rnorm(n1,j*2,0.2)
}
S1=S1+matrix(rnorm(n1*n2,0,sigma1),n1,n2)
image(1:n2,1:n1,t(S1[n1:1,]),xlab="",ylab="",main="Source 1")
contour(1:n2,1:n1,t(S1[n1:1,]),add=TRUE)
# Second source
val1=1
val2=1.2
val3=1.5
val4=2
sigma2=0.1
S2=matrix(0,n1,n2)
S2[2:5,4:10]=val1
S2[3:4,5:9]=val3
S2[13:18,16:19]=val2
S2[14:17,17:18]=val4
S2=S2+matrix(rnorm(n1*n2,0,sigma2),n1,n2)
image(1:n2,1:n1,t(S2[n1:1,]),xlab="",ylab="",main="Source 2")
contour(1:n2,1:n1,t(S2[n1:1,]),add=TRUE)
# Generating data matrix X
A = rerow(matrix(runif(M^2)-0.5,M,M))
W = solve(A)
S=rbind(as.vector(S1),as.vector(S2))
X = A %*% S
# Solving Blind Source Separation problem with three different methods
cica = cICA(X,tol=0.001)
## scica = scICA(X,n1=n1,n2=n2,h=(2*pi/10),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))
image(t(S1[n1:1,]),xlab="",ylab="")
contour(t(S1[n1:1,]),add=TRUE)
image(t(S2[n1:1,]),xlab="",ylab="")
contour(t(S2[n1:1,]),add=TRUE)
image(t(matrix(Shat1[1,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat1[1,],n1,n2)[n1:1,]),add=TRUE)
image(t(matrix(Shat1[2,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat1[2,],n1,n2)[n1:1,]),add=TRUE)
## par(mfrow=c(2,2))
## image(t(S1[n1:1,]),xlab="",ylab="")
## contour(t(S1[n1:1,]),add=TRUE)
## image(t(S2[n1:1,]),xlab="",ylab="")
## contour(t(S2[n1:1,]),add=TRUE)
## image(t(matrix(Shat2[1,],n1,n2)[n1:1,]),xlab="",ylab="")
## contour(t(matrix(Shat2[1,],n1,n2)[n1:1,]),add=TRUE)
## image(t(matrix(Shat2[2,],n1,n2)[n1:1,]),xlab="",ylab="")
## contour(t(matrix(Shat2[2,],n1,n2)[n1:1,]),add=TRUE)
par(mfrow=c(2,2))
image(t(S1[n1:1,]),xlab="",ylab="")
contour(t(S1[n1:1,]),add=TRUE)
image(t(S2[n1:1,]),xlab="",ylab="")
contour(t(S2[n1:1,]),add=TRUE)
image(t(matrix(Shat3[1,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat3[1,],n1,n2)[n1:1,]),add=TRUE)
image(t(matrix(Shat3[2,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat3[2,],n1,n2)[n1:1,]),add=TRUE)
## End (Not run)
Example output
Loading required package: MASS
Loading required package: fastICA
[1] "Color ICA - Iteration 1: error is equal to 0.112973512675102"
[1] "Color ICA - Iteration 2: error is equal to 0.00647369072744824"
[1] "Color ICA - Iteration 3: error is equal to 0.000454218285733443"
[1] 0.2634561
[1] 0.2434067
coloredICA documentation built on May 1, 2019, 10:55 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note Author(s) References See Also Examples
Description
This function implements the spatial colored Independent Component Analysis (scICA) algorithm, where sources are treated as spatial stochastic processes on a lattice.
Usage
1 2 |
Arguments
Xin |
Data matrix with |
M |
Number of components to be extracted. |
Win |
Initial guess for the unmixing matrix |
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 |
unmixing.estimate |
The method used in the unmixing matrix estimation step. The two allowed choices are |
n1 |
Number of rows of the lattice. |
n2 |
Number of columns of the lattice. |
nx01 |
Number of rows of the lattice where the spectral density is evaluate. Default value is |
nx02 |
Number of columns of the lattice where the spectral density is evaluate. Default value is |
h |
Kernel bandwidth used for the nonparametric estimation of the sources spectral densities. |
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 scICA 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}.
Spatial colored Independent Component Analysis assumes that the independent sources are spatial stochastic processes on a lattice. 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 spatial autocorrelated sources as follows:
l(W,\boldsymbol{f}_{\boldsymbol{S}};\widetilde{X})=∑_{j=1}^p∑_{k=1}^n≤ft(\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 τ\bold{w}'_{j}C_{j}\bold{w}_{j}, where \bold{w}'_{j} is the jth column of W, C_j=∑_{k\neq j}\bold{w}_k\bold{w}'_k and τ is a tuning parameter. The matrix C_j provides an orthogonality constraint in the sense that \bold{w}'_{j}C_{j}\bold{w}_{j}=∑_{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 setting
unmixing.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 setting
unmixing.matrix = "newton".
Independently from the choice of the technique to minimize the objective function, the scICA 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:
nonparametric estimation of the sources spectral density through a multidimensional local linear kernel estimator \widehat{m}_{LK} (see Shen et al. (2014) for further details).
estimate the unmixing matrix W according the method selected in
unmixing.estimate.
Value
A list containing the following components:
W |
Estimate of the |
K |
pre-whitening matrix that projects data onto the first |
A |
Estimate of the |
S |
Estimate of the |
X |
Original |
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 |
Note
Note that source matrix S and spectral density matrix den are n x M matrices. Every row, that should be in a n1 x n2 grid, has been vectorized in a n vector by column, with n = n1 x n2.
Author(s)
Lee, S., Shen, H., Truong, Y. and Zanini, P.
References
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.
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.
See Also
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 | ## Not run:
require(fastICA)
n1=20
n2=20
M=2
# Fist source
sigma1=2
S1=matrix(0,n1,n2)
for (i in 1:n1){
S1[i,]=rnorm(n2,i*2,0.2)
}
for (j in 1:n2){
S1[,j]=S1[,j]+rnorm(n1,j*2,0.2)
}
S1=S1+matrix(rnorm(n1*n2,0,sigma1),n1,n2)
image(1:n2,1:n1,t(S1[n1:1,]),xlab="",ylab="",main="Source 1")
contour(1:n2,1:n1,t(S1[n1:1,]),add=TRUE)
# Second source
val1=1
val2=1.2
val3=1.5
val4=2
sigma2=0.1
S2=matrix(0,n1,n2)
S2[2:5,4:10]=val1
S2[3:4,5:9]=val3
S2[13:18,16:19]=val2
S2[14:17,17:18]=val4
S2=S2+matrix(rnorm(n1*n2,0,sigma2),n1,n2)
image(1:n2,1:n1,t(S2[n1:1,]),xlab="",ylab="",main="Source 2")
contour(1:n2,1:n1,t(S2[n1:1,]),add=TRUE)
# Generating data matrix X
A = rerow(matrix(runif(M^2)-0.5,M,M))
W = solve(A)
S=rbind(as.vector(S1),as.vector(S2))
X = A %*% S
# Solving Blind Source Separation problem with three different methods
cica = cICA(X,tol=0.001)
## scica = scICA(X,n1=n1,n2=n2,h=(2*pi/10),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))
image(t(S1[n1:1,]),xlab="",ylab="")
contour(t(S1[n1:1,]),add=TRUE)
image(t(S2[n1:1,]),xlab="",ylab="")
contour(t(S2[n1:1,]),add=TRUE)
image(t(matrix(Shat1[1,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat1[1,],n1,n2)[n1:1,]),add=TRUE)
image(t(matrix(Shat1[2,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat1[2,],n1,n2)[n1:1,]),add=TRUE)
## par(mfrow=c(2,2))
## image(t(S1[n1:1,]),xlab="",ylab="")
## contour(t(S1[n1:1,]),add=TRUE)
## image(t(S2[n1:1,]),xlab="",ylab="")
## contour(t(S2[n1:1,]),add=TRUE)
## image(t(matrix(Shat2[1,],n1,n2)[n1:1,]),xlab="",ylab="")
## contour(t(matrix(Shat2[1,],n1,n2)[n1:1,]),add=TRUE)
## image(t(matrix(Shat2[2,],n1,n2)[n1:1,]),xlab="",ylab="")
## contour(t(matrix(Shat2[2,],n1,n2)[n1:1,]),add=TRUE)
par(mfrow=c(2,2))
image(t(S1[n1:1,]),xlab="",ylab="")
contour(t(S1[n1:1,]),add=TRUE)
image(t(S2[n1:1,]),xlab="",ylab="")
contour(t(S2[n1:1,]),add=TRUE)
image(t(matrix(Shat3[1,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat3[1,],n1,n2)[n1:1,]),add=TRUE)
image(t(matrix(Shat3[2,],n1,n2)[n1:1,]),xlab="",ylab="")
contour(t(matrix(Shat3[2,],n1,n2)[n1:1,]),add=TRUE)
## End (Not run)
|
Example output
Loading required package: MASS
Loading required package: fastICA
[1] "Color ICA - Iteration 1: error is equal to 0.112973512675102"
[1] "Color ICA - Iteration 2: error is equal to 0.00647369072744824"
[1] "Color ICA - Iteration 3: error is equal to 0.000454218285733443"
[1] 0.2634561
[1] 0.2434067
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.