In [50]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

Inverse Wavelet Transformation

In [55]:
# constants

#filter 1
c_bits = np.array([1+np.sqrt(3), 3+np.sqrt(3), 3-np.sqrt(3), 1-np.sqrt(3)]) / (4.0*np.sqrt(2))
# #filter 2
c = np.array([-c_bits[3], -c_bits[2], c_bits[1], -c_bits[0]])


def inverse_transform(x):
    # get the filter matrix
    N = x.shape[-1]
    F = np.zeros((N,N))
    row1 = np.zeros(N)
    row1[:4] = c_bits
    
    row2 = np.zeros(N)
    row2[:4] = c

    
    F[:,::2] = linalg.circulant(row1)[:,::2]
    F[:,1::2] = linalg.circulant(row2)[:,::2]
    
    #print F.T

    # splitting matrix
    S = np.zeros((N,N))
    v = np.zeros(N/2)
    v[0] = 1
    S[::2, :N/2] = linalg.circulant(v)
    S[1::2,N/2:] = linalg.circulant(v)
    final = np.zeros(len(S[:,N-1]))
    final[1] = 1    
    S[:,N-1]= final
    
    #print S.T

    inv = np.dot(np.dot(F,S),x)
    
    
    return inv

def run_inverse():
    x = np.zeros(2**12)
    x[4] = 1
    x[29] = 1
    
    for i in range(2,13):
        inp = x.copy()[:2**i]
        x[:2**i] = inverse_transform(inp)
    
    plt.figure()
    plt.plot(x)
    plt.show()
In [56]:
run_inverse()

PCA: Numerical Verification of Eigenvalues and Covariance

In [68]:
# generate points from gaussian distribution
sigma = 1
mu = 0

N = 1000000
vals = np.random.randn(2,N)

# x1, x2, x3
x = np.vstack((vals[0], vals[1], vals[0] + vals[1]))

# get covariance matrix
C = np.sum(x*x.reshape(-1,1,N), axis=-1) / N
print "Covariance Matrix\n", C

# find the eigenvalues
w,v = np.linalg.eig(C)
print "Eigenvalues\n", w
print "Eigenvectors\n", v

print "------------------------------"
# new covariance and eigen metrics
E = v.T[:2] # don't need the 0 eigenvector
x_p = np.dot(E,x)
Cx_p = np.sum(x_p*x_p.reshape(-1,1,N), axis=-1) / N
print "New Covariance\n", Cx_p
w_xp, v_xp = np.linalg.eig(Cx_p)

print "Eigenvalues\n", w_xp
print "Eigenvectors\n", v_xp

Covariance Matrix
[[  1.00128495e+00  -5.16585612e-04   1.00076836e+00]
 [ -5.16585612e-04   1.00052408e+00   1.00000750e+00]
 [  1.00076836e+00   1.00000750e+00   2.00077586e+00]]
Eigenvalues
[  3.00116401e+00   1.00142088e+00   3.26948305e-16]
Eigenvectors
[[ -4.08481264e-01   7.06972223e-01   5.77350269e-01]
 [ -4.08015273e-01  -7.07241263e-01   5.77350269e-01]
 [ -8.16496537e-01  -2.69040196e-04  -5.77350269e-01]]
------------------------------
New Covariance
[[  3.00116401e+00   1.14683019e-15]
 [  1.14683019e-15   1.00142088e+00]]
Eigenvalues
[ 3.00116401  1.00142088]
Eigenvectors
[[  1.00000000e+00  -5.73488751e-16]
 [  5.73488751e-16   1.00000000e+00]]

ICA

In [78]:
# plot pairs of uniform random variables
N = 1000
s = np.random.rand(N,2)
print s.shape
plt.figure()
plt.scatter(s[:,0], s[:,1])
plt.show()

(1000, 2)
In [79]:
# mix them with a square matrix and plot
A = [[1,2],[3,1]]
x = np.dot(s,A)
plt.figure()
plt.scatter(x[:,0], x[:,1])
plt.show()
In [113]:
# make zero mean, diagonalize with unit variance and plot
x_mean = x - np.mean(x)
# make column vectors
x_mean = x_mean.T
c = np.cov(x_mean)
print c
w,v = np.linalg.eig(c)
x_diag = np.dot((v/np.sqrt(w)).T, x_mean)
print x_diag.shape
plt.figure()
plt.scatter(x_diag[0], x_diag[1])
plt.show()
print x_diag.shape

[[ 0.85848294  0.43511274]
 [ 0.43511274  0.42896427]]
(2, 1000)

(2, 1000)
In [116]:
# find the independant components of x with log cosh contrast function and plot
def f(x):
    return np.log(cosh(x))

def df(x):
    return np.tanh(x)

def ddf(x):
    return 1/np.cosh(x)**2
    

def do_ICA(w, x):
    i = 0
    max_count = 5
    thresh = 1e-5
    while i <= max_count:
        # get new weight
        temp = np.dot(w, x)
        w_new = np.mean(x * df(np.dot(w, x))) - w*np.mean(ddf(np.dot(w,x)))
        # normalize
        w_new = w_new / np.linalg.norm(w_new)
        # exit conditions
        if np.linalg.norm(w_new - w) < thresh:
            break
        w = w_new
        i+=1
        
    return w

w_f = do_ICA(np.random.rand(2), x_diag)
print w_f
plt.figure()
plt.scatter(s[:,0], s[:,1], c="r")
plt.scatter(x[:,0], x[:,1], c='b')
plt.scatter(x_diag[0], x_diag[1], c='g')
plt.plot([0, w_f[0]], [0, w_f[1]], marker='.', c='y')
plt.plot([0, -w_f[1]], [0, w_f[0]], marker='.', c='y')
plt.show()

[ 0.85656615 -0.51603724]
In [ ]: