Step by step...: Independent Component Analysis. ICA Neural Network. Java implementation

Independent component analysis (ICA) aims to solve problem of signals separation from their linear mixture. ICA is a special case of blind source separation, when separation performed without the aid of information (or with very little information) about the source signals or the process of signal mixing. Although blind source separation problem in general is underdetermined, the useful solution can be obtained under a certain assumptions.

ICA model assumes that there are

independent signals

and some mixing matrix

We can observe only signals

that are represented by linear superposition of

The main difficulty lies in the fact that both

and

are unknown. To solve the problem ICA uses assumption that signals in mixture are statistically independent and has non-Gaussian probability distributions.

In general, two random variables

and

are statistically independent if information about one variable says nothing about other variable. From mathematical point of view statistical independence means that 2-D probability density function

is the product of 1-D probability density functions:

For statistically independent signals the covariance matrix of odd functions

and

is diagonal: all mutual covariances are zeros:

and all self covariances are non-zeros:

There are a several methods exists to solve ICA problem. One of them is uses special Artificial Neural Network. The structure of such Neural Network for signal separation is shown on Figure 1.

Figure1 - Neural network structure

The mixed signals

are only sources of information for network. The output vector

is obtained by linear transformation of input signals

by system of weights

The matrix

is defined by learning rule. There are a lot of learning rules performed in real-time (see sources). One of well-performed rules is natural gradient-based rule:

where

is learning coefficient. Functions

are odd. In practice there are various choices of these functions. Often

and

I present here Java implementation of ICA Neural Network.

Class NeuralICA performs independent component analysis using neural network. It contains following methods:

Public methods:

1. Method double[] mult (double[][]a, double []b) receives matrix a and vector b, performs matrix-vector multiplication and returns result vector.

2. Method double[][] mult( double[][] a, double[][] b ) receives matrixes a and b, performs matrix-matrix multiplication and returns result matrix.

3. Method double[][] sum (double[][]a, double [][]b) receives matrices a and b,

performs matrix element-wise summation and returns result matrix.

4. Method void update (double[][]a, double [][]b) receives matrices a and b,

performs element-wise summation and rewrites matrix a with result.

5. Method double[] ica(double [] x) receives vector x of observations for a certain time moment and returns update of output vector of separated signals for this time. It performs two operations:

Static methods:

1. Method static double[] getCol(double [][] a, int col) receives matrix a and column index col and returns the column from matrix a with index col.

2. Method static void setCol(double [][] a, double[] b, int col) receives matrix a, vector b and column index col and sets the column in matrix a with values from vector b.

3. Method static void main(String args[]) runs Neural network ICA algorithm on test signals mixture.

Private methods:

1. Method double[] F(double[] y) computes function

2. Method double[] G(double[] y) computes function

3.Method double[][] dW(double[]f, double[]g) receives functions

and

and performs

and returns

Class NeuralICA:

package com.blogspot.shulgadim.ica;

import java.util.Random;

public class NeuralICA {
    private double [][] W;  
    private double [] y;        
    private double[][]I;            
    private double eta = 0.01;
        
    public NeuralICA(int N){        
        I = new double[N][N];
        for(int i=0; i<I.length; i++){              
            for(int j=0; j<I[0].length; j++){
                if(i==j){
                    I[i][j] = 1.0;
                }
            }
        }                   
        W = new double[N][N];           
        for(int i=0; i<N; i++){             
            for(int j=0; j<N; j++){     
                if(i==j){
                    W[i][j] = 1.0;
                }
            }
        }                       
        y = new double[N];
    }
    
    
    public double[] ica(double [] x){               
        y = mult(W,x);              
        update(W,  mult(dW(F(y),G(y)), W ) );               
        return y;
    }
    
    public  double[] mult (double[][]a, double []b){
        int M = a.length;
        int N = a[0].length;
        double[] y = new double[M];
        for(int i=0; i<M; i++){                 
            for(int j=0; j<N; j++){
                y[i] += a[i][j]*b[j];
            }           
        }
        return y;
    }
 
    
    public  double[][] sum (double[][]a, double [][]b){
        int M = a.length;
        int N = a[0].length;        
        double[][] c = new double[M][N];
        for(int i=0; i<M; i++){                 
            for(int j=0; j<N; j++){
                c[i][j] = a[i][j] + b[i][j];
            }           
        }   
        return c;
    }

    public double[][] mult( double[][] a, double[][] b ){
        int m = a.length;
        int n = a[0].length;
        int o = b[0].length;
        double[][] matres = new double[m][];
        for ( int i = 0; i < m; ++i ){
            matres[i] = new double[o];
            for ( int j = 0; j < o; ++j ){
                matres[i][j] = 0.0f;
                for ( int k = 0; k < n; ++k )
                    matres[i][j] += a[i][k] * b[k][j];
            }
        }
        return(matres);
    }
    
    public  void update (double[][]a, double [][]b){
        int M = a.length;
        int N = a[0].length;        
        for(int i=0; i<M; i++){                 
            for(int j=0; j<N; j++){
                a[i][j] += b[i][j];
            }           
        }   
    }
    
    private double[] F(double[] y){     
        double[] f = new double[y.length];
        for(int i=0; i<y.length; i++){          
            f[i] = y[i];
        }       
        return f;       
    }
    
    private double[] G(double[] y){     
        double[] g = new double[y.length];
        for(int i=0; i<y.length; i++){          
            g[i] = Math.tanh(10*y[i]);          
        }       
        return g;
    }
    
    private double[][] dW(double[]f, double[]g){                
        double[][] A = new double[f.length][f.length];
        for(int i=0; i<f.length; i++){              
            for(int j=0; j<f.length; j++){              
                A[i][j] = eta*(I[i][j] - f[i]*g[j]);                
            }
        }                       
        return A;
    }
    
    public static void main(String args[]){
        int K = 5000;       
        double [][] s = new double[3][K];   
        double dt = 0.001;      
        Random rand = new Random();
        for(int k=0; k< K; k++){
            s[0][k] = rand.nextGaussian();        
            s[1][k] = 
            Math.pow(10, 0)*Math.sin(500*dt*k)*Math.sin(30*dt*k);
            s[2][k] = Math.pow(10, 0)*Math.sin(234*dt*k);
            
        }
        new MultiPlot("s",s);
        double[][] A = {{0.20, -0.61, 0.62},
                          {0.91, -0.89, -0.33},
                          {0.59,  0.76, 0.60}};             
        double [][] x = new double[3][K];       
        for(int k=0; k<K; k++){
            for(int i=0; i<3; i++){             
                for(int j=0; j<3; j++){                                     
                    x[i][k] = x[i][k] + A[i][j]*s[j][k];                                                                        
                }
            }
            
        }       
        new MultiPlot("x",x);                           
        NeuralICA ica = new NeuralICA(3);       
        double [][] y = new double[3][K];       
        for(int k=0; k<K; k++){
            double[] res = ica.ica(getCol(x,k));
            setCol(y,res,k);
            
        }       
        new MultiPlot("y",y);       
    }
    
    public static double[] getCol(double [][] a, int col){
        int M = a.length;
        double [] y = new double[M];
        for(int i =0; i<M; i++){
            y[i] = a[i][col];
        }
        return y;
    }
    
    public static void setCol(double [][] a, double[] b, int col){
        int M = b.length;       
        for(int i =0; i<M; i++){
            a[i][col] = b[i];
        }       
    }
}

Class MultiPlot plots signals in one window for convenient viewing. It uses JFreeChart library.

package com.blogspot.shulgadim.ica;

import java.awt.*;
import javax.swing.*;
import org.jfree.chart.*;
import org.jfree.chart.plot.*;
import org.jfree.data.xy.*;

public class MultiPlot {    
    private String name;           
    private JFrame frame;   
    public MultiPlot(String name, double y[][]) {
        this.name = name;     
        frame = new JFrame("MuliPlot");                                                        
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);        
        JPanel panel = new JPanel();
        panel.setLayout(new BoxLayout(panel, BoxLayout.Y_AXIS));        
        for(int row=0; row< y.length; row++){           
            XYDataset dataset = createDataset(y,row);
            JFreeChart chart = createChart(dataset, row);
            ChartPanel chartPanel = new ChartPanel(chart);        
            chartPanel.setPreferredSize(new Dimension(1400,100));
            panel.add(chartPanel);                
        }        
        frame.getContentPane().add(new JScrollPane(panel));        
        frame.setSize(1200, 600);        
        frame.setVisible(true);        
    }
            
    public MultiPlot(String name, double y[][][]) {
        this.name = name;     
        frame = new JFrame("MuliPlot");                                                        
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);        
        JPanel panel = new JPanel();
        panel.setLayout(new BoxLayout(panel, BoxLayout.Y_AXIS));        
        int M = y.length;
        int N = y[0].length;                                   
        int row = 0;
        for(int i=0; i<M; i++){
            for(int j=0; j<N; j++){
                XYDataset dataset = createDataset(y,i,j);
                JFreeChart chart = createChart(dataset, row);
                ChartPanel chartPanel = new ChartPanel(chart);        
                chartPanel.setPreferredSize(new Dimension(1400,100));
                panel.add(chartPanel);
                row++;
            }
        }             
        frame.getContentPane().add(new JScrollPane(panel));        
        frame.setSize(1200, 600);        
        frame.setVisible(true);        
    }
    
    
    private XYDataset createDataset(double[][] y, int row) {
        XYSeries series = new XYSeries("");
        for(int i=0; i<y[0].length;i++){
            series.add(i,y[row][i]);   
        }        
        XYSeriesCollection dataset = new XYSeriesCollection();
        dataset.addSeries(series);                                    
        return dataset;        
    }
        
    private XYDataset createDataset(double[][][] y, int i, int j) {
        XYSeries series = new XYSeries("");
        int K = y[0][0].length;
        for(int k=0; k<K;k++){
            series.add(k,y[i][j][k]);   
        }        
        XYSeriesCollection dataset = new XYSeriesCollection();
        dataset.addSeries(series);                                    
        return dataset;      
    }
    
    private JFreeChart createChart(final XYDataset dataset, int row) {
               
        final JFreeChart chart = ChartFactory.createXYLineChart(
            name + String.valueOf(row),                         //chart title
            "",                        // x axis label
            "",        // y axis label
            dataset,                    // data
            PlotOrientation.VERTICAL,
            false,                      // include legend
            false,                      // tooltips
            false                       // urls
        );       
        chart.setBackgroundPaint(Color.LIGHT_GRAY);        
        XYPlot plot = chart.getXYPlot();        
        plot.setBackgroundPaint(Color.WHITE);             
        plot.setDomainGridlinePaint(Color.LIGHT_GRAY);
        plot.setRangeGridlinePaint(Color.LIGHT_GRAY);                         
        return chart;        
    }
}

Figure 2 shows three input signals: one noise signal and two deterministic signals. Figure 3 shows the signals mixture and Figure 4 shows signals separation performing in time. Fromg Figure 4 it is clear, that at the beging of processing the neural network is learned, then is starts separating signals quite distinctly.