In [23]:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import time
sns.set(style='whitegrid', rc={'figure.figsize':(14,8)})
In [92]:
PERIOD = 0.25
N = 2500
x = np.linspace(0, PERIOD, N)
freq1 = 697
wave1 = np.sin(2 * np.pi * freq1 * x)
freq2 = 1209
wave2 = np.sin(2 * np.pi * freq2 * x)
t = wave1 + wave2
plt.plot(t)
plt.show()
(b) Perform Direct Cosine Transform¶
In [31]:
def DCT(t):
N = len(t)
D = np.zeros((N, N))
D[0, :] = np.sqrt(1 / N)
for i in range(1, N):
for j in range(N):
D[i, j] = np.sqrt(2 / N) * np.cos((np.pi * (2*j + 1) * i) / (2 * N))
f = D.dot(t)
return f, D
In [32]:
start = time.time()
f, D = DCT(t)
print("Naive DCT took %.1f seconds" % (time.time() - start))
In [33]:
plt.plot(f)
plt.show()
(c) Inverse DCT¶
In [35]:
D_inv = D.T
t_hat = D_inv.dot(f)
plt.plot(t_hat)
plt.title("Reconstructed signal after DCT")
plt.show()
diff = np.abs(t - t_hat)
print("Max difference from original signal: ", np.max(diff))
plt.plot(diff)
plt.title("Difference")
plt.show()
(d) Random sample %5 of the signal¶
In [96]:
M = int(N * 0.05)
indexes = sorted(np.random.choice(N, M, replace=False))
x_sample = [x[i] for i in indexes]
t_sample = [t[i] for i in indexes]
D_sample = D[:, indexes]
print(D_sample.shape)
plt.plot(x_sample, t_sample)
plt.plot(x_sample, t_sample, 'o')
plt.title("5% Sample of Signal")
plt.show()
(e) DCT using gradient descent¶
$$ E = \sum_{k = 0}^{M - 1} \Big( t'_k - \sum_{i = 0}^{N - 1} D_{ik} f'_i \Big)^2 $$
In [99]:
def sample_error(t_sample, f_sample):
s = 0
for k in range(M):
s += (t_sample[k] - np.dot(D_sample[:, k], f_sample)) ** 2
return s
In [100]:
f_sample = np.random.normal(0, 1, N)
print("Sample error for random guess: ", sample_error(t_sample, f_sample))
$$ \frac{\partial E}{\partial f_i} = -2 \sum_{k = 0}^{M - 1} \Big( t'_k - \sum_{i = 0}^{N - 1} D_{ik} f'_i \Big) D_{ik} $$
In [101]:
def sample_error_grad(t_sample, f_sample):
g = np.zeros(N)
for i in range(N):
for k in range(M):
g[i] += (t_sample[k] - np.dot(D_sample[:, k], f_sample)) * D_sample[i, k]
g *= -2
return g
In [114]:
f_sample = np.random.normal(0, 1, N)
last_error = sample_error(t_sample, f_sample)
learning_rate = 0.1
iterations = 0
while learning_rate > 1e-6:
iterations += 1
g = sample_error_grad(t_sample, f_sample)
new_f = f_sample - learning_rate * g
new_error = sample_error(t_sample, new_f)
print(iterations, "\t\t", new_error)
if last_error - new_error < 1e-4:
learning_rate *= 0.1
print("New learning rate: ", learning_rate)
if new_error < last_error:
f_sample = new_f
last_error = new_error
In [115]:
print("Sample error after", iterations, "iterations of GD:", sample_error(t_sample, f_sample))
plt.plot(f_sample)
plt.title("Coefficients")
plt.show()
t_hat = D_inv.dot(f_sample)
plt.plot(t_hat)
plt.title("Reconstructed signal after subsampled DCT")
plt.show()
(f) Same, with L2 norm regularizarion¶
In [104]:
def sample_error_l2(t_sample, f_sample):
s = np.sum(f_sample ** 2)
for k in range(M):
s += (t_sample[k] - np.dot(D_sample[:, k], f_sample)) ** 2
return s
In [105]:
f_sample = np.random.normal(0, 1, N)
print("Sample error for random guess: ", sample_error_l2(t_sample, f_sample))
$$ \frac{\partial E_{L2}}{\partial f_i} = -2 \sum_{k = 0}^{M - 1} \Big( t'_k - \sum_{i = 0}^{N - 1} D_{ik} f'_i \Big) D_{ik} + 2f'_i $$
In [106]:
def sample_error_l2_grad(t_sample, f_sample):
g = np.zeros(N)
for i in range(N):
g[i] = 2 * f_sample[i]
for k in range(M):
g[i] -= 2 * (t_sample[k] - np.dot(D_sample[:, k], f_sample)) * D_sample[i, k]
return g
In [112]:
f_sample = np.random.normal(0, 1, N)
last_error = sample_error_l2(t_sample, f_sample)
learning_rate = 0.1
iterations = 0
while learning_rate > 1e-6:
iterations += 1
g = sample_error_l2_grad(t_sample, f_sample)
new_f = f_sample - learning_rate * g
new_error = sample_error_l2(t_sample, new_f)
print(iterations, "\t\t", new_error)
if last_error - new_error < 1e-4:
learning_rate *= 0.1
print("New learning rate: ", learning_rate)
if new_error < last_error:
f_sample = new_f
last_error = new_error
In [113]:
print("Sample error after", iterations, "iterations of GD:", sample_error_l2(t_sample, f_sample))
plt.plot(f_sample)
plt.title("Coefficients")
plt.show()
t_hat = D_inv.dot(f_sample)
plt.plot(t_hat)
plt.title("Reconstructed signal after subsampled DCT (L2 regularized)")
plt.show()
In [109]:
f_l2 = f_sample.copy()
Same, with L1 reguarization¶
In [82]:
def sample_error_l1(t_sample, f_sample):
s = np.sum(np.abs(f_sample))
for k in range(M):
s += (t_sample[k] - np.dot(D_sample[:, k], f_sample)) ** 2
return s
In [83]:
f_sample = np.random.normal(0, 1, N)
print("Sample error for random guess: ", sample_error_l1(t_sample, f_sample))
$$ \frac{\partial E_{L1}}{\partial f_i} = -2 \sum_{k = 0}^{M - 1} \Big( t'_k - \sum_{i = 0}^{N - 1} D_{ik} f'_i \Big) D_{ik} + sign(f'_i) $$
Where I define $sign$ as:
sign($x$) = -1 if $x < tolerance$, 1 if $x > tolerance$, else 0
In [127]:
def sample_error_l1_grad(t_sample, f_sample, tol=1e-6):
f_copy = f_sample.copy()
f_copy[np.abs(f_copy) < tol] = 0
g = np.sign(f_copy)
for i in range(N):
for k in range(M):
g[i] -= 2 * (t_sample[k] - np.dot(D_sample[:, k], f_sample)) * D_sample[i, k]
return g
In [128]:
f_sample = np.random.normal(0, 1, N)
last_error = sample_error_l1(t_sample, f_sample)
learning_rate = 0.1
iterations = 0
while learning_rate > 1e-6:
iterations += 1
g = sample_error_l1_grad(t_sample, f_sample)
new_f = f_sample - learning_rate * g
new_error = sample_error_l1(t_sample, new_f)
print(iterations, "\t\t", new_error)
if last_error - new_error < 1e-4:
learning_rate *= 0.1
print("New learning rate: ", learning_rate)
if new_error < last_error:
f_sample = new_f
last_error = new_error
In [129]:
print("Sample error after", iterations, "iterations of GD:", sample_error_l1(t_sample, f_sample))
plt.plot(f_sample)
plt.title("Coefficients")
plt.show()
t_hat = D_inv.dot(f_sample)
plt.plot(t_hat)
plt.title("Reconstructed signal after subsampled DCT (L1 regularized)")
plt.show()
Inverse Wavelet Transform¶
In [159]:
c0 = (1 + np.sqrt(3)) / (4 * np.sqrt(2))
c1 = (3 + np.sqrt(3)) / (4 * np.sqrt(2))
c2 = (3 - np.sqrt(3)) / (4 * np.sqrt(2))
c3 = (1 - np.sqrt(3)) / (4 * np.sqrt(2))
def invert_wavelet_transform(v):
N = len(v)
split = np.zeros((N, N))
for i in range(N // 2):
split[i, i * 2] = 1
split[i + (N // 2), i * 2 + 1] = 1
v2 = split.T.dot(v)
transform = np.zeros((N, N))
for i in range(0, N, 2):
transform[i, i] = c0
transform[i, i + 1] = c1
transform[i, (i + 2) % N] = c2
transform[i, (i + 3) % N] = c3
transform[i + 1, i] = c3
transform[i + 1, i + 1] = -c2
transform[i + 1, (i + 2) % N] = c1
transform[i + 1, (i + 3) % N] = -c0
v3 = transform.T.dot(v2)
return v3
In [160]:
x = np.zeros(1)
w = np.zeros(2 ** 12)
w[4] = 1
w[29] = 1
while len(x) < 2 ** 12:
tmp_w = w[:len(x)]
v = np.concatenate([x, tmp_w])
x = invert_wavelet_transform(v)
In [161]:
plt.plot(x)
plt.show()
(13.3)¶
(a) Analytically compute covariance:¶
$$ X_1 \sim N(0, 1) \\ X_2 \sim N(0, 1) \\ X_3 = X1 + X2 $$
The diagonal: $$ C_{11} = Var(X_1) = 1 \\ C_{22} = Var(X_2) = 1 \\ C_{33} = Var(X_3) = Var(X1) + Var(X1) = 2 $$
The rest: $$ C_{12} = \langle(X_1 - \langle X_1\rangle)(X_2 - \langle X_2\rangle\rangle = \langle X_1 X_2\rangle = 0 \\ C_{13} = \langle(X_1 - \langle X_1\rangle)(X_3 - \langle X_3\rangle\rangle = \langle X_1^2 + X_1 X_2\rangle = \langle X_1^2\rangle + \langle X_1 X_2\rangle = 1 \\ C_{23} = \langle(X_2 - \langle X_2\rangle)(X_3 - \langle X_3\rangle\rangle = \langle X_2^2 + X_1 X_2\rangle = \langle X_2^2\rangle + \langle X_1 X_2\rangle = 1 $$
$$ C = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 2 \end{pmatrix} $$
(b) Eigenvalues¶
$$ det \begin{pmatrix} 1 - \lambda & 0 & 1 \\ 0 & 1 - \lambda & 1 \\ 1 & 1 & 2 - \lambda \end{pmatrix} = 0 $$
$$ (1 - \lambda)((1 - \lambda)(2 - \lambda) - 1) - (1 - \lambda) = 0 \implies \\ (1 - \lambda)(2 -3 \lambda + \lambda^2 - 2) = 0 \implies \\ \lambda(1 - \lambda)(\lambda - 3) = 0 $$
Eigenvalues: $$ \lambda_1 = 0 \\ \lambda_2 = 1 \\ \lambda_3 = 3 $$
(c) Numerically validate:¶
In [202]:
np.set_printoptions(precision=2)
np.set_printoptions(suppress=True)
N = 10000
x1 = np.random.normal(0, 1, N)
x2 = np.random.normal(0, 1, N)
x3 = x1 + x2
data = np.array([x1, x2, x3])
cov = np.cov(data)
print("Covariance matrix:")
print(cov)
eigVals, eigVecs = np.linalg.eig(cov)
print("\nEigenvalues:")
print(eigVals)
(d) Transform¶
In [203]:
trans_data = eigVecs[:, :2].T.dot(data)
trans_cov = np.cov(trans_data)
print("Transformed Covariance matrix:")
print(trans_cov)
trans_eigVals, _ = np.linalg.eig(trans_cov)
print("\nEigenvalues:")
print(trans_eigVals)
In [297]:
N = 10000
s1 = np.random.uniform(0, 1, N)
s2 = np.random.uniform(0, 1, N)
plt.plot(s1, s2, 'o')
plt.show()
(b)¶
In [298]:
s = np.array([s1, s2])
A = np.array([[1, 2], [3, 1]])
x = A.dot(s)
plt.plot(x[0], x[1], 'o')
plt.show()
(c)¶
In [299]:
# zero mean
x = x - np.mean(x)
# diagnolize
cov = np.cov(x)
eigVals, eigVecs = np.linalg.eig(cov)
x = eigVecs.T.dot(x)
# scale to 1
x = (x.T / np.sqrt(eigVals)).T
# verify:
print("Cov matrix:")
print(np.cov(x))
plt.plot(x[0], x[1], 'o')
plt.show()
(d) ICA¶
In [333]:
def f(x):
return np.log(np.cosh(x))
# according to Wolfram Alpha..
def f_tag(x):
return np.tanh(x)
# again, according to Wolfram Alpha
def f_tag_tag(x):
return 1 / (np.cosh(x) ** 2)
w = np.random.normal(0, 1, 2)
iterations = 0
while iterations < 100:
iterations += 1
w_new = np.mean(x * f_tag(np.dot(w, x))) - w * np.mean(f_tag_tag(np.dot(w, x)))
w_new = w_new / np.linalg.norm(w_new)
if np.linalg.norm(w + w_new) < 1e-4:
break
w = w_new.copy()
print("Finished after", iterations, "iterations")
print(w)
ortho_w = np.array([1, -w[0] / w[1]])
ortho_w = ortho_w / np.linalg.norm(ortho_w)
print(ortho_w)
In [334]:
plt.plot(x[0], x[1], 'o')
plt.plot([0, w[0]], [0, w[1]], color='red')
plt.plot([0, ortho_w[0]], [0, ortho_w[1]], color='red')
plt.show()
