Chapter 13 - Transformations

(13.1) Prove that the DFT is unitary.

In the matrix form, each element of the Discrete Fourier Transform is $$M_{fn} = \frac{e^{2\pi ifn/N}}{\sqrt{N}} $$

We need to show that $M^\dagger M=I$ or $(M^\dagger M)_{fn} = \delta_{fn}$.

We have $$ \begin{align} (M^\dagger M)_{fn} &= \sum_{k=0}^{N-1}M^{\dagger}_{fk}M_{kn} \\ &=\frac{1}{N}\sum_{k=0}^{N-1}e^{2 \pi ik(f-n)/N} \end{align}$$

For $f=n$, we have $(M^\dagger M)_{ff} =\frac{1}{N}\sum_{k=0}^{N-1}1 = 1 $.

For $f\neq n$, we have $$\begin{align} (M^\dagger M)_{fn} &= \frac{1}{N}\sum_{k=0}^{N-1}(e^{2 \pi i(f-n)})^{k/N} \\ &= \frac{1}{N}\sum_{k=0}^{N-1}[\cos(2 \pi (f-n))+i \sin(2 \pi (f-n))]^k \\ &= 0 \end{align}$$ , since $f-n$ always integer.

Combining the results, we come to the conclusion that $(M^\dagger M)_{fn} = \delta_{fn}$ and thus $M$ unitary, i.e. the DFT is a unitary transformation.

(12.2) Calculate the inverse wavelet transform, using Daubechies fourth-order coefficients, of a vector of length $2^{12}$, with a $1$ in the 5th and 30th places and zeros elsewhere.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="ticks")
import numpy as np

L=2**12
# Fourth-order coefficients
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))

v = np.zeros((L))
v[4] = 1
v[29] = 1
v_next = np.zeros_like(v)

# I calculate the inverse wavelet transform by the matrix in page 155, and applying the cascade algorithm on the 
# reverse direction, starting from two x elements (2^1), all the way up to the full (2^12) x elements 
for j in range(1, 12):
    N = 2**(j+1)  # Current number of elements in correct positions
    # Assuming that we start as the last column (two x elements and then w elements)
    x = np.copy(v[:N//2])
    w = np.copy(v[N//2:N])  # LOL spent 3.5 hours here to figure out I need np.copy.....
    # Re-shuffle
    v[:N:2] = x
    v[1:N:2] = w
    # Calculate new column by inverting the transform
    # First two lines
    v_next[0] = c0 * v[0] + c3 * v[1] + c2 * v[L-2] + c1 * v[L-1]
    v_next[1] = c1 * v[0] - c2 * v[1] + c3 * v[L-2] - c0 * v[L-1]
    for i in range(2,N,2):
        # odd lines
        v_next[i] = c2 * v[i-2] + c1 * v[i-1] + c0 * v[i] + c3 * v[i+1]
        # even lines
        v_next[i+1] = c3 * v[i-2] - c0 * v[i-1] + c1 * v[i] - c2 * v[i+1]
    v[:N] = v_next[:N]
plt.plot(v)
plt.xlim([0, L])
plt.ylim([min(v)-.01, max(v)+.01])
plt.title('4th-order Daubechies Wavelet')
plt.show()

(12.3) Consider a measurement of a three-component vector $\vec{x}$, with $x_1$ and $x_2$ being drawn independently from a Gaussian distribution with zero mean and unit variance, and $x_3 = x_1 + x_2$.

(a) Analytically calculate the covariance matrix of $\vec{x}$.

We have $\mu_{\vec{x}} = \langle \vec{x} \rangle = (\begin{matrix}0 & 0 & 0\end{matrix}) $,

$\langle x_{\{1,2\}}^2 \rangle = Var(x_{\{1,2\}}) - \langle x_{\{1,2\}} \rangle = 1$,

$ \langle x_{3}^2 \rangle = Var(x_{3}) - \langle x_{3} \rangle = Var(x_1)+Var(x_2) = 2 $,

$\langle x_ix_j\rangle = \langle x_i\rangle\langle x_j\rangle = 0$ for $i,j$ in $\{1,2\}$ since $x_1$, $x_2$ independent and, lastly,

$\langle x_3x_{\{1,2\}}\rangle = \langle x_{\{1,2\}}^2\rangle + \langle x_{\{1,2\}}x_{\{2,1\}}\rangle = 1$.

Thus, $$\begin{align} \mathbf{C}_x &= \big\langle (\vec{x}-\mu_{\vec{x}})\cdot (\vec{x}-\mu_{\vec{x}})^T \big\rangle \\ &= \big\langle \vec{x}\cdot \vec{x}^T \big\rangle \\ &= \begin{bmatrix} \langle x_1^2 \rangle & \langle x_1x_2 \rangle & \langle x_1x_3 \rangle \\ \langle x_2x_1 \rangle & \langle x_2^2 \rangle & \langle x_2x_3 \rangle \\ \langle x_3x_1 \rangle & \langle x_3x_2 \rangle & \langle x_3^2 \rangle \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 2 \end{bmatrix} \end{align}$$

(b) What are the eigenvalues?

To calculate the eigenvalues we need to solve the characteristic equation. In illustration, $$\begin{align} |\mathbf{C}_x - \lambda \mathbf{I}| &= 0 \\ \begin{vmatrix} 1-\lambda & 0 & 1 \\ 0 & 1-\lambda & 1 \\ 1 & 1 & 2-\lambda \end{vmatrix} &= 0 \\ (1-\lambda)[(1-\lambda)(2-\lambda)-1]-(1-\lambda) &=0 \end{align}$$ which gives three solutions (i.e. eigenvalues): $\lambda = 0$, $\lambda = 1$, $\lambda= 3$.

(c) Numerically verify these results by drawing a data set from the distribution and computing the covariance matrix and eigenvalues.

In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="ticks")
import numpy as np
x1 = np.random.normal(0,1,20000)
x2 = np.random.normal(0,1,20000)
x = np.array([x1, x2, x1+x2])
C_x = np.cov(x)
eig = np.linalg.eigvals(C_x)
print('Covariance Matrix')
print(np.round(C_x,2))
print('Eigenvalues')
print(np.round(eig,2))

Covariance Matrix
[[ 0.99  0.01  1.  ]
 [ 0.01  1.    1.01]
 [ 1.    1.01  2.01]]
Eigenvalues
[ 3.02  0.99  0.  ]

(d) Numerically find the eigenvectors of the covariance matrix, and use them to construct a transformation to a new set of variables $\vec{y}$ that have a diagonal covariance matrix with no zero eigenvalues. Verify this on the data set.

In [3]:
# Calculate eigenvectors
eigvec = np.linalg.eig(C_x)[1]
# Select the ones with non-zero eigenvalues
new_eigvec = eigvec[eig > 1e-6]
# Transform the dataset
y = np.matmul(new_eigvec, x)
# Print new outputs
C_y = np.cov(y)
eigy = np.linalg.eigvals(C_y)
print('Covariance Matrix')
print(np.round(C_y,2))
print('Eigenvalues')
print(np.round(eigy,2))

Covariance Matrix
[[ 2.64  0.81]
 [ 0.81  0.98]]
Eigenvalues
[ 2.97  0.65]

(12.4) Generate pairs of uniform random variables $\{s_1, s_2\}$ with each component contained in $[0, 1]$.

(a) Plot these data.

In [4]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="white")
import numpy as np

s1 = np.random.uniform(0, 1, (5000,1))
s2 = np.random.uniform(0, 1, (5000,1))
s = np.hstack((s1, s2))
plt.scatter(s[:,0],s[:,1])
sns.despine()
plt.show()

(b) Mix them $(\vec{x} = \mathbf{A}\cdot\vec{s})$ with a square matrix $A = \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}$ and plot.

In [5]:
A = np.array([[1, 2],[3, 3]])
x = np.matmul(A,s.T).T
plt.scatter(x[:,0],x[:,1])
sns.despine()
plt.show()

(c) Make $\vec{x}$ zero mean, diagonalize with unit variance, and plot.

In [6]:
# Make zero mean
x_m = x - np.mean(x, 0, keepdims=True)

# Find covariance
C_x = np.cov(x_m.T)

# Calculate eigenvalues and vectors
eigval, eigvec = np.linalg.eig(C_x)

# Transformation diagonalize and remove variance
M = (eigvec / np.sqrt(eigval)).T

# # Apply transformation
y = np.matmul(M, x_m.T)

plt.scatter(y[0,:],y[1,:])
sns.despine()
plt.show()

(d) Find the independent components of $\vec{x}$ with the log cosh contrast function, and plot.

In [7]:
# Performing iterative ICA
# max number of iters
n_iters = 100

def f(x):
    # Contrast function
    return np.log(np.cosh(x))

def df(x):
    # Derivative of contrast function
    return np.tanh(x)

def sdf(x):
    # Second derivative of contrast function
    return (1/np.cosh(x))**2

def E(x):
    return np.mean(x, 1, keepdims=True)

# Initialize w (assuming that we want to find 2 components, perpendicular to each other)
w = np.random.random(size=(2,1))
w_next = np.zeros_like(w)
for i in range(n_iters):
    w_next = E(y*df(np.dot(w.T, y))) - E(sdf(np.matmul(w.T, y)))*w
    w = w_next / np.linalg.norm(w_next)
    
    if np.abs(np.dot(w.T, w_next) - 1) < 1e-1:
        break
        
plt.scatter(y[0,:],y[1,:])
plt.plot([0, w[0]], [0, w[1]], 'r')
plt.plot([0, -w[1]], [0, w[0]], 'r')
sns.despine()
plt.show()