Wolfram Research

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 matrix can be square or rectangular.
The columns of the argument matrix are interpreted as signals.
Here are the options taken by ResourceFunction["IndependentComponentAnalysis"]:
"InitialUnmixingMatrix" Automatic initial un-mixing matrix
"NegEntropyFactor" 1 negative entropy factor
"NonGaussianityFunction" Automatic non-Gaussianity measure function
"RowNorm" False whether to normalize the matrix rows at each step
MaxSteps 200 maximum steps of the ICA interation process
PrecisionGoal 6 precision 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

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[5]=

Here is the matrix product of the obtained factors:

In[6]:=
MatrixForm[A.S]
Out[6]=

Scope

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

In[7]:=
SeedRandom[232];
mat = RandomReal[{10, 100}, {10, 2}];
mat2 = RandomReal[{-10, 10}, {10, 7}];
mat = ArrayPad[mat, {{0, 0}, {0, 5}}] + mat2;
ListLinePlot[Transpose[mat], PlotRange -> All, ImageSize -> Medium, PlotTheme -> "Detailed"]
Out[8]=

Compute the ICA matrix factors:

In[9]:=
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[10]:=
Norm[mat - A.S]/Norm[mat]
Out[10]=

Here is the relative error for the first three columns:

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

Here are comparison plots:

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

Options

InitialUnmixingMatrix

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[13]:=
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]}]];

Here we plot the results:

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

NonGaussianityFunction

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[15]:=
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[16]:=
Grid[KeyValueMap[{#1[x], ListLinePlot[Transpose[#2[[1]]], PlotRange -> All, ImageSize -> 250, PlotTheme -> "Detailed"]} &, res], Alignment -> {Center, Left}]
Out[16]=

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

Blind signal separation

Here are signal generating functions:

In[17]:=
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[18]:=
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[19]:=
nSize = 600;
S = Table[{s1[t], s2[t], s3[t]}, {t, 0, nSize, 0.5}];
Dimensions[S]
Out[10]=
In[20]:=
Grid[{Map[ListLinePlot[#, PlotRange -> All, ImageSize -> 250] &, Transpose@S]}]
Out[20]=

Create the mixed signals matrix:

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

Find the ICA decomposition:

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

Approximation norm:

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

Visualize the ICA-identified source signals:

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

Possible Issues

Outliers presence sensitivity

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[26]:=
M1 = M + RandomChoice[{0.005, 0.995} -> {1, 0}, Dimensions[M]];
MapThread[
 ListLinePlot[{##}, PlotLegends -> {"with outliers", "original"}] &, {Transpose[M1], Transpose[M]}]
Out[16]=

Here we apply ICA to the new mixed signals:

In[27]:=
{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[28]:=
Grid[{Map[ListLinePlot[#, PlotRange -> All, ImageSize -> 250] &, Transpose[A]]}]
Out[28]=

Resource History

Source Metadata

License Information