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 [ ]: