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