[PYTHON] ICA with Scipy

I wrote ICA for the first time in a while while studying numpy. When I try to use it, it is unexpectedly troublesome to handle array and matrix. Is there any good way?

fast_ica.py


#!/usr/bin/env python
# -*- coding:utf-8 -*-
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt

def g( val, alpha = 0.01 ):
    return np.tanh( alpha * val )

def gd( val, alpha = 0.01 ):
    return alpha * ( 1.0 - np.asarray( np.tanh( alpha * val ) ) )

def sym_orth( sig ):
    a,b,c = np.linalg.svd( sig * sig.T )
    return np.dot( np.dot( a, np.dot( np.diag( np.sqrt( 1./ b ) ), c ) ), sig )

def fast_ica( sig ):
    m,n = sig.shape

    #Whitening
    u,s,v = np.linalg.svd( sig, full_matrices=False)
    v = np.sqrt( n ) * v

    #Isolation filter initialization
    w = np.random.rand( m, m )
    w = w / np.matrix( [ [ np.linalg.norm( w[ :, x ] ) ] * m for x in xrange( m ) ] ).T
    w = sym_orth( w )

    #ICA
    err = 1.0
    it = 0
    max_iter = 3000
    while 1e-12 < min( err , it < max_iter ):
        last_w = w
        w1 =  np.dot( v , g( np.dot( w.T , v ) ).T )
        w2 = np.multiply( np.dot( gd( np.dot( w.T, v ) ), np.ones( (n, 1 ) ) ) , w )
        w = ( w1 - w2 ) / n
        w = sym_orth( w )
        err = np.sum( np.abs( np.abs(w) - np.abs(last_w) ) )
        it += 1

    #Create a separation filter based on the gain of the 0th signal
    W = np.dot( np.diag( np.ravel( u[0,] ) ), np.dot( np.diag( s ) , w.T ) )
    return np.dot( W, v )  / np.sqrt( n )

def main():
    #Super Gaussian distribution with Laplace distribution
    s = np.random.laplace( 0., 1.0, ( 2, 10000 ) )
    a = np.asmatrix( np.random.rand( 2,2 ) )
    x = a * s
    y = np.asarray( fast_ica( x ) )
    plt.scatter( y[0,], y[1,] )
    plt.show()
    return

if __name__ == '__main__':
    main()

Recommended Posts

ICA with Scipy
LPC with Scipy
CORDIC with Scipy
Create filter with scipy
Normarize data with Scipy
Use OpenBLAS with numpy, scipy
Extract peak values with scipy
Harmonic mean with Python Harmonic mean (using SciPy)
Calculate sample distribution with Scipy (discrete distribution)
Generate a normal distribution with SciPy
scipy Voronoi
Install Scipy
scipy stumbles with pip install on python 2.7.8
Use multithreaded BLAS / LAPACK with numpy / scipy
Using Intel MKL with NumPy / SciPy (November 2019 version)
One liner to make Lena images with scipy
Create 3D scatter plot with SciPy + matplotlib (Python)