[PYTHON] ICA mit Scipy

Ich habe zum ersten Mal seit einiger Zeit ICA geschrieben, um Numpy zu lernen. Wenn ich versuche, es zu verwenden, ist es unerwartet schwierig, mit Array und Matrix umzugehen. Gibt es einen guten Weg?

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

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

    #Initialisierung des Isolationsfilters
    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

    #Erstellen Sie ein Trennfilter basierend auf der Verstärkung des 0. Signals
    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-Gauß-Verteilung mit Laplace-Verteilung
    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()