Function Repository Resource:

IndependentComponentAnalysis

Source Notebook

Decompose a matrix into Independent Component Analysis matrix factors

Contributed by: Anton Antonov

ResourceFunction["IndependentComponentAnalysis"][mat,k]

decomposes the matrix mat into k components.

ResourceFunction["IndependentComponentAnalysis"][mat,k,opts]

decomposes mat using the options opts.

Details and Options

The Independent Component Analysis (ICA) is a matrix factorization method that utilizes Gram-Schmidt orthogonalization by considering two vectors (two signals) to be orthogonal if their difference is Gaussian white noise.
The argument mat can be a square or rectangular matrix.
The columns of mat are interpreted as signals.
ResourceFunction["IndependentComponentAnalysis"] accepts the following options:
"InitialUnmixingMatrix"Automaticinitial un-mixing matrix
"NegEntropyFactor"1negative entropy factor
"NonGaussianityFunction"Automaticnon-Gaussianity measure function
"RowNorm"Falsewhether to normalize the matrix rows at each step
MaxSteps200maximum steps of the ICA interation process
PrecisionGoal6precision goal
For some data, at least several executions are required in order to get components that provide satisfactory interpretation.
The negative entropy factor is specified with the option "NegEntropyFactor". The negative entropy factor controls how aggressively ICA learns from the "landscape" produced by the non-Gaussianity measure function (specified with "NonGaussianityFunction").
The function ResourceFunction["IndependentComponentAnalysis"] uses SingularValueDecomposition internally.

Examples

Basic Examples (3) 

Here is a random integer matrix:

In[1]:=
SeedRandom[7]
mat = RandomInteger[10, {4, 3}];
MatrixForm[mat]
Out[3]=

Here are the Independent Component Analysis matrix factors:

In[4]:=
{A, S} = ResourceFunction["IndependentComponentAnalysis"][mat, 3];
Row[{MatrixForm[A], MatrixForm[S]}]
Out[4]=

Here is the matrix product of the obtained factors:

In[5]:=
MatrixForm[A . S]
Out[8]=

Scope (6) 

Here is a random matrix with its first two columns having much larger magnitudes than the rest:

In[9]:=
SeedRandom[232];
mat = RandomReal[{10, 100}, {10, 2}];
mat2 = RandomReal[{-10, 10}, {10, 7}];
mat = ArrayPad[mat, {{0, 0}, {0, 5}}] + mat2;

Plot the values:

In[10]:=
ListLinePlot[Transpose[mat], PlotRange -> All, ImageSize -> Medium, PlotTheme -> "Detailed"]
Out[10]=

Compute the ICA matrix factors:

In[11]:=
SeedRandom[22];
{A, S} = ResourceFunction["IndependentComponentAnalysis"][mat, 3, MaxSteps -> 100, PrecisionGoal -> 8];

Here is the relative error of the approximation by the obtained matrix factorization:

In[12]:=
Norm[mat - A . S]/Norm[mat]
Out[12]=

Here is the relative error for the first three columns:

In[13]:=
Norm[mat[[All, 1 ;; 3]] - (A . S)[[All, 1 ;; 3]]]/
 Norm[mat[[All, 1 ;; 3]]]
Out[13]=

Here are comparison plots:

In[14]:=
Block[{opts = {PlotRange -> All, ImageSize -> 250, PlotTheme -> "Detailed"}},
 Row[{ListLinePlot[Transpose[mat], opts, PlotLabel -> "Orginal"], ListLinePlot[Transpose[A . S], opts, PlotLabel -> "Reconstructed"]}]
 ]
Out[14]=

Options (4) 

InitialUnmixingMatrix (2) 

Using the option "InitialUnmixingMatrix", we can influence the results by providing initial un-mixing matrix for ICA's iteration process. Here we compute ICA with different initial un-mixing matrices:

In[15]:=
k = 3;
SeedRandom[8];
res = Association[
   MapThread[
    Row[{#2, ":", Spacer[2], If[MatrixQ[#], MatrixForm[#1], ToString[#1]]}] -> ResourceFunction["IndependentComponentAnalysis"][mat, k, "InitialUnmixingMatrix" -> #1] &, {{Automatic, Automatic, RandomVariate[NormalDistribution[1, 10], {k, k}], SingularValueDecomposition[mat, k][[2]]}, Range[4]}]];

Plot the results:

In[16]:=
Grid[KeyValueMap[{#1, ListLinePlot[Transpose[#2[[1]]], PlotRange -> All, ImageSize -> 250, PlotTheme -> "Detailed"]} &, res], Alignment -> {Center, Left}]
Out[16]=

NonGaussianityFunction (2) 

ICA is essentially a Gram-Schmidt process that considers two signals to be orthogonal if their difference is Gaussian noise.

Here we compute ICA with different non-Gaussianity measure functions:

In[17]:=
res = Association[
   # -> BlockRandom[
       SeedRandom[22];
       ResourceFunction["IndependentComponentAnalysis"][mat, 3, "NonGaussianityFunction" -> #, MaxSteps -> 2]
       ] & /@ {Log[
       Cosh[#]] &, -Exp[-#^2/2] &, (Log[Cosh[1.2 #]]/1.2 &)}];

Plot the results:

In[18]:=
Grid[KeyValueMap[{#1[x], ListLinePlot[Transpose[#2[[1]]], PlotRange -> All, ImageSize -> 250, PlotTheme -> "Detailed"]} &, res], Alignment -> {Center, Left}]
Out[18]=

Note that we used a small number of steps in order to obtain similar but different results. Using a sufficiently large number of steps (say, 20) we get almost the same results (with the matrix mat).

Applications (8) 

Blind signal separation (8) 

Here are signal generating functions:

In[19]:=
Clear[s1, s2, s3]
s1[t_] := Sin[600 \[Pi] t/10000 + 6*Cos[120 \[Pi] t/10000]] + 1.2
s2[t_] := Sin[\[Pi] t/10] + 1.2
s3[t_?NumericQ] := (((QuotientRemainder[t, 23][[2]] - 11)/9)^5 + 2.8)/
   2 + 0.2

Here is a mixing matrix:

In[20]:=
A = {{0.44, 0.2, 0.31}, {0.45, 0.8, 0.23}, {0.12, 0.32, 0.71}};

Form a matrix of the signals:

In[21]:=
nSize = 600;
S = Table[{s1[t], s2[t], s3[t]}, {t, 0, nSize, 0.5}];
Dimensions[S]
Out[21]=

Plot the values:

In[22]:=
Grid[{Map[ListLinePlot[#, PlotRange -> All, ImageSize -> 250] &, Transpose@S]}]
Out[22]=

Create and plot the mixed signals matrix:

In[23]:=
M = S . A;
Grid[{Map[ListLinePlot[#, PlotRange -> All, ImageSize -> 250] &, Transpose@M]}]
Out[23]=

Find the ICA decomposition:

In[24]:=
{A, S} = ResourceFunction["IndependentComponentAnalysis"][M, 3];

Approximation norm:

In[25]:=
Norm[M - A . S]
Out[25]=

Visualize the ICA-identified source signals:

In[26]:=
Grid[{Map[ListLinePlot[#, PlotRange -> All, ImageSize -> 250] &, Transpose[A]]}]
Out[26]=

Possible Issues (3) 

Outliers presence sensitivity (3) 

The ICA computations can be sensitive to the presence of outliers. Here we add a small number of outliers to the original mixed signals:

In[27]:=
M1 = M + RandomChoice[{0.005, 0.995} -> {1, 0}, Dimensions[M]];
MapThread[
 ListLinePlot[{##}, PlotLegends -> {"with outliers", "original"}] &, {Transpose[M1], Transpose[M]}]
Out[27]=

Here we apply ICA to the new mixed signals:

In[28]:=
{A, S} = ResourceFunction["IndependentComponentAnalysis"][M1, 3];

We see that the found components are not that good compared to the components found for the original mixed signals data:

In[29]:=
Grid[{Map[ListLinePlot[#, PlotRange -> All, ImageSize -> 250] &, Transpose[A]]}]
Out[29]=

Publisher

Anton Antonov

Version History

  • 1.0.0 – 18 November 2019

Source Metadata

License Information