# COMPUTATIONAL METHODS FOR DATA SCIENCE HW1
郭訓佑, R11946025

## 0
#### code

# Q0: A Practice on the Randomization
import random

# read tsv file
def read_tsv(file_name) -> list:
    l = []
    with open(file_name, 'r') as f:
        for line in f:
            a = line.strip().split(' ')
            # remove empty string
            a = [x for x in a if x != '']
            l.append(a)
    return l

# Q1a
"""
Write a simple code (any programming language of your choice) to randomly pick 6
objects out of 12, where all 12 possible objects have the same probability to be picked
"""
# read tsv file
a = read_tsv("E:\\GitHub\\NTU_Computational_Methods_for_Data_Science_2023\\hw1\\staff.stat.sinica.edu.tw_fredphoa_HW_HW1_CAmaxTemp.txt")
# remove column 0,1, 14
a = [x[2:14] for x in a]
# check it is a 12 * 12 matrix
assert len(a) == 12, "the matrix must be a 12 * 12 matrix"
for i in range(12):
    assert len(a[1]) == 12, "the matrix must be a 12 * 12 matrix"
    
# randomly pick(not replace) 6 rows and 6 columns and sort them
# set seed
random.seed(0)
columns = random.sample(range(12), 6)
rows = random.sample(range(12), 6)
columns.sort()
rows.sort()

# Q0b
"""
Run the code on both the stations and the months, and state which stations and months
you obtain as the reduced data X.
"""

print("Pikced columns: ", columns)
print("Pikced rows: ", rows)
X = [[a[i][j] for j in columns] for i in rows]
# %%
# Q0c
"""
Do a quick check if X has at least four REAL eigenvalues. You can do this question by
using existing library or packages implemented in the software of your choice.
"""
# check if X has at least four REAL eigenvalues
import numpy as np
from numpy import linalg as LA
X_np = np.array(X, dtype=float)
eigenvalues, eigenvectors = LA.eig(X_np)
real_eigenvalues = []
for i in range(len(eigenvalues)):
    if np.isreal(eigenvalues[i]):
        real_eigenvalues.append(eigenvalues[i])
print("# of real eigenvalues: ", len(real_eigenvalues))
assert len(real_eigenvalues) >= 4, "X must have at least four REAL eigenvalues"
# %%
print("X: ")
print(X_np)

#### Ans:
Picked Stations:
23225BLUE_CANYON,CA
93193FRESNO,CA
23129LONG_BEACH,CA
93134LOS_ANGELES,CA
23232SACRAMENTO,CA
23190SANTA_BARBARA,CA
Picked Months:
1,4,5,7,8,12

Pikced columns:  [0, 3, 4, 6, 7, 11]
Pikced rows:  [2, 3, 4, 5, 7, 10]
# of real eigenvalues:  4

X: 
[[ 76.  82.  88.  95.  97.  75.]
 [ 78. 100. 107. 113. 112.  77.]
 [ 93. 105. 104. 109. 105.  92.]
 [ 95. 106. 102. 108. 105.  92.]
 [ 76.  95. 105. 114. 112.  74.]
 [ 86.  96. 101. 109. 105.  89.]]

## 1
### a
Run the LU factorization on X to obtain a lower triangular matrix L and an upper triangular matrix U.
#### code:

import numpy as np
import copy
np.set_printoptions(precision=3)
np.set_printoptions(suppress=True)

# Q1.a
"""
Run the LU factorization on X to obtain a lower triangular matrix L and an upper
triangular matrix U.
"""
def LU(X: np.array)-> tuple:
    # check if X is a square matrix
    assert len(X) == len(X[0]), "X must be a square matrix"
    n = len(X)
    # initialize L and U
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    # calculate L and U
    for i in range(n):
        L[i][i] = 1
        for j in range(i, n):
            U[i][j] = X[i][j] - sum([L[i][k] * U[k][j] for k in range(i)])
        for j in range(i+1, n):
            L[j][i] = (X[j][i] - sum([L[j][k] * U[k][i] for k in range(i)])) / U[i][i]
    return L, U

L, U = LU(X_np)
print("Q1.a")
print("L: ")
print(L)
print("U: ")
print(U)

# %%
# check if LU decomposition is correct
print("Check if LU decomposition is correct")
print("L * U: ")
print(L.dot(U))

#### Ans:

L: 
[[ 1.     0.     0.     0.     0.     0.   ]
 [ 1.026  1.     0.     0.     0.     0.   ]
 [ 1.224  0.294  1.     0.     0.     0.   ]
 [ 1.25   0.221  1.36   1.     0.     0.   ]
 [ 1.     0.821 -0.385  0.917  1.     0.   ]
 [ 1.132  0.203  0.228  0.558  0.962  1.   ]]
U: 
[[ 76.     82.     88.     95.     97.     75.   ]
 [  0.     15.842  16.684  15.5    12.447   0.026]
 [  0.      0.     -8.59  -11.807 -17.357   0.216]
 [  0.      0.      0.      1.889   4.614  -2.05 ]
 [  0.      0.      0.      0.     -6.132   0.941]
 [  0.      0.      0.      0.      0.      4.315]]
 
 Check if LU decomposition is correct
L * U: 
[[ 76.  82.  88.  95.  97.  75.]
 [ 78. 100. 107. 113. 112.  77.]
 [ 93. 105. 104. 109. 105.  92.]
 [ 95. 106. 102. 108. 105.  92.]
 [ 76.  95. 105. 114. 112.  74.]
 [ 86.  96. 101. 109. 105.  89.]]

### b
Use Power Iteration method to find the largest eigenvalue-eigenvector pair of X.
#### code:

def largest_eigen_pair(X: np.array, epsilon: float = 1e-6, max_iter=1000000)-> tuple:
    # power method
    assert len(X) == len(X[0]), "X must be a square matrix"
    n = len(X)
    # initialize x
    x = np.random.rand(n)
    v = x / np.linalg.norm(x)
    
    # eignevalue = v.T * X * v
    eigen = np.dot(v.T, np.dot(X, v))
    num_iter = 0
    while True:
        num_iter += 1
        # update v
        v = np.dot(X, v)
        v = v / np.linalg.norm(v)
        # update eigenvalue
        eigen_new = np.dot(v.T, np.dot(X, v))
        # check if converge
        if abs(eigen_new - eigen) < epsilon or num_iter > max_iter:
            if num_iter > max_iter:
                print("Not converge")
            break
        if num_iter % 100 == 0:
            print("Iteration: ", num_iter)
            
        eigen = eigen_new
        num_iter += 1
    return eigen, v

def QR(X: np.array) -> tuple:
    # QR decomposition using Gram-Schmidt process
    Q = np.zeros(X.shape)
    R = np.zeros((X.shape[1], X.shape[1]))
    for k in range(X.shape[1]):
        for i in range(k):
            R[i][k] = np.dot(Q[:, i].T, X[:, k])
        qk = copy.deepcopy(X[:, k])
        for i in range(k):
            qk = qk - R[i][k] * Q[:, i]
        R[k][k] = np.linalg.norm(qk)
        Q[:, k] = qk / R[k][k]
    return Q, R

largest_value, largest_vector = largest_eigen_pair(X_np)
print("Q1.b")
print("Largest eigenvalue: ")
print(largest_value)
print("Largest eigenvector: ")
print(largest_vector)
# %%
# check 
print("Check if largest eigenvalue and eigenvector is correct")
print("X * v: ")
print(X_np.dot(largest_vector))
print("lambda * v: ")
print(largest_value * largest_vector)

#### Ans:

Largest eigenvalue: 
581.7047038772532
Largest eigenvector: 
[0.361 0.414 0.427 0.427 0.406 0.412]

Check if largest eigenvalue and eigenvector is correct
X * v: 
[209.707 240.581 248.406 248.26  236.129 239.643]
lambda * v: 
[209.707 240.581 248.406 248.26  236.129 239.643]

### c
Use QR factorization to find all eigenvectors with REAL eigenvalues of X.
#### code:

def real_eigen_pairs(X: np.array, max_iter=500000)-> tuple:
    # QR iteration
    """
    if matrix X is symmetric, then the eigenvalues are real and the eigenvectors is almost correct
    if matrix X is not symmetric, then the eigenvalues are complex and the eigenvectors is not correct
    """
    assert len(X) == len(X[0]), "X must be a square matrix"
    # check if X is symmetric
    if not np.allclose(X, X.T):
        print("X is not symmetric, the eigenvalues are complex and the eigenvectors is not correct")
    
    n = len(X)
    # initialize Q and R
    Q = np.identity(n)
    R = X
    num_iter = 0
    Q_eigen = np.identity(n)
    while True:
        num_iter += 1
        # update Q and R
        Q, R = QR(R.dot(Q))
        Q_eigen = Q_eigen.dot(Q)
        # check if converge
        if num_iter > max_iter:
            break
        if num_iter % 100 == 0:
            print("Iteration: ", num_iter)

    eigenvalues = []
    for i in range(n):
        eigenvalues.append(R[i][i])
    
    eigenvectors = []
    for i in range(n):
        eigenvectors.append(Q_eigen[:, i])
    return eigenvalues, eigenvectors

eigenvalues, eigenvectors = real_eigen_pairs(X_np)
print("Q1.c")
print("Eigenvalues: ")
print(eigenvalues)
print("Eigenvectors: ")
print(eigenvectors)
# %%
# check
print("Check if eigenvalues and eigenvectors are correct")
for i in range(len(eigenvalues)):
    print("X * v: ")
    print(X_np.dot(eigenvectors[i]))
    print("lambda * v: ")
    print(eigenvalues[i] * eigenvectors[i])

#### Ans:
I found some cases will break this method which the eigenvectors are difficult to converge. e.f. QR = RQ or some asymmetric.

Eigenvalues: 
[581.7047038564525, 9.393027168397762, 3.4800347037085526, 8.627255256768713, 3.915093807405027, 0.8047459721312528]
Eigenvectors: 
[array([0.361, 0.414, 0.427, 0.427, 0.406, 0.412]), array([-0.109,  0.585, -0.311, -0.185,  0.522, -0.492]), array([-0.533, -0.173,  0.063, -0.414,  0.476,  0.533]), array([-0.646, -0.101,  0.378,  0.548,  0.064, -0.355]), array([ 0.342, -0.663, -0.171,  0.189,  0.574, -0.219]), array([-0.201,  0.079, -0.738,  0.526, -0.044,  0.36 ])]

### d
Find the inverse of X.

#### code:

# Q1.d
"""
Find the inverse of X.
"""
def inverse_LU_version(X: np.array)-> np.array:
    # find inverse of X using LU decomposition
    assert len(X) == len(X[0]), "X must be a square matrix"
    n = len(X)
    L, U = LU(X)
    # initialize inverse of X
    X_inv = np.zeros((n, n))
    # solve for each column of X_inv
    for i in range(n):
        b = np.zeros(n)
        b[i] = 1
        # solve Ly = b
        y = np.zeros(n)
        for j in range(n):
            y[j] = b[j] - sum([L[j][k] * y[k] for k in range(j)])
        # solve Ux = y
        x = np.zeros(n)
        for j in range(n-1, -1, -1):
            x[j] = (y[j] - sum([U[j][k] * x[k] for k in range(j+1, n)])) / U[j][j]
        X_inv[:, i] = x
    return X_inv
X_inv = inverse_LU_version(X_np)
print("Q1.d")
print("Inverse of X: ")
print(X_inv)
# %%
# check
print("Check if inverse of X is correct")
print("X * X_inv: ")
print(X_np.dot(X_inv))

#### Ans:
I use LU decomposition to do invertion.

Q1.d
Inverse of X: 
[[ 0.102 -0.302  0.397 -0.172  0.256 -0.269]
 [-0.114  0.204 -0.391  0.316 -0.146  0.118]
 [ 0.14  -0.097  0.822 -0.623  0.063 -0.292]
 [-0.19  -0.208 -0.241  0.217  0.24   0.165]
 [ 0.113  0.208 -0.303  0.161 -0.197  0.036]
 [-0.034  0.19  -0.241  0.075 -0.223  0.232]]
 
 Check if inverse of X is correct
X * X_inv: 
[[ 1.  0.  0.  0. -0.  0.]
 [-0.  1.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0. -0.]
 [-0.  0. -0.  1.  0.  0.]
 [-0.  0.  0.  0.  1.  0.]
 [-0.  0.  0.  0.  0.  1.]]

## 2
### a
Center the data and compute the variance-covariance matrix.
#### code:

def cov_matrix(X: np.array)-> np.array:
    # Center the data and compute the variance-covariance matrix.
    # center X
    X = copy.deepcopy(X - np.mean(X, axis=0))
    # calculate variance-covariance matrix of X
    n = len(X[0])
    # initialize variance-covariance matrix
    var_cov = np.zeros((n, n))
    # calculate variance-covariance matrix
    for i in range(n):
        for j in range(i, n):
            var_cov[i][j] = np.dot(X[:, i].T, X[:, j]) / (n-1)
            var_cov[j][i] = var_cov[i][j]
    return var_cov

X_mean = X_np - np.mean(X_np, axis=0)
var_cov = cov_matrix(X_np)
print("Q2.a")
print("Centered data: ")
print(X_mean)
print("Variance-covariance matrix: ")
print(var_cov)

#### Ans:

Q2.a
Centered data: 
[[ -8.    -15.333 -13.167 -13.     -9.     -8.167]
 [ -6.      2.667   5.833   5.      6.     -6.167]
 [  9.      7.667   2.833   1.     -1.      8.833]
 [ 11.      8.667   0.833   0.     -1.      8.833]
 [ -8.     -2.333   3.833   6.      6.     -9.167]
 [  2.     -1.333  -0.167   1.     -1.      5.833]]
Variance-covariance matrix: 
[[74.    57.4   14.8    7.4   -6.8   72.8  ]
 [57.4   76.667 47.533 41.    25.    53.333]
 [14.8   47.533 46.167 45.2   34.6   13.567]
 [ 7.4   41.    45.2   46.4   36.2    7.   ]
 [-6.8   25.    34.6   36.2   31.2   -8.4  ]
 [72.8   53.333 13.567  7.    -8.4   75.767]]

### b
Find the top three principal components using power iteration. Calculate the cumulative
percentage of the total eigenvalues that these three principal components cover.
#### code:

def deflated_power_method(X: np.array, n: int,epsilon: float = 1e-6, max_iter=1000000)-> tuple:
    # power method with deflation to find the first n largest eigenvalues and eigenvectors
    assert len(X) == len(X[0]), "X must be a square matrix"
    # initialize eigenvalues and eigenvectors
    eigenvalues = []
    eigenvectors = []
    for i in range(n):
        eigen, vec = largest_eigen_pair(X, epsilon, max_iter)
        eigenvalues.append(eigen)
        eigenvectors.append(vec)
        # deflation
        X = X - eigen * np.outer(vec, vec)
    return eigenvalues, eigenvectors

def find_topn_principal_components(X: np.array, n: int)-> np.array:
    # find the top n principal components of X
    # calculate covariance matrix of X
    C = cov_matrix(X)
    # find the top n eigenvalues and eigenvectors
    eigenvalues, eigenvectors = deflated_power_method(C, n)
    # sort eigenvalues and eigenvectors
    eigenvalues = np.array(eigenvalues)
    eigenvectors = np.array(eigenvectors)
    idx = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[idx]
    # return the top n principal components and eigenvalues
    return eigenvectors, eigenvalues

# find the top 3 principal components of X
PCn, eigvals = find_topn_principal_components(X_np, 3)
cov_X = cov_matrix(X_np)
all_eigenvalues, _ = deflated_power_method(cov_X, 10)
print("Q2.b")
print("Top 3 principal components: ")
print(PCn)
print("Cumulative percentage of the total eigenvalues that these three principal components cover: ", sum(eigvals) / sum(all_eigenvalues))

#### Ans:

Q2.b
Top 3 principal components: 
[[ 0.487  0.57   0.329  0.283  0.148  0.477]
 [-0.412  0.111  0.426  0.477  0.469 -0.435]
 [ 0.176  0.632 -0.083 -0.443 -0.029 -0.605]]
Cumulative percentage:  0.9973567855222676

### c
Plot the data on a 3D space with three principal component axes. Provide the coordinates of the recast data

#### code:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# plot original data (first three columns) with three principal components as axes
fig = plt.figure()
ax = Axes3D(fig)
# plot three principal component times eigenvalues as axes
ax.plot([0, PCn[0][0] * np.sqrt(eigvals[0])*2], [0, PCn[0][1] * np.sqrt(eigvals[0])*2], [0, PCn[0][2] * np.sqrt(eigvals[0])*2], color='r')
ax.plot([0, PCn[1][0] * np.sqrt(eigvals[1])*2], [0, PCn[1][1] * np.sqrt(eigvals[1])*2], [0, PCn[1][2] * np.sqrt(eigvals[1])*2], color='g')
ax.plot([0, PCn[2][0] * np.sqrt(eigvals[2])*2], [0, PCn[2][1] * np.sqrt(eigvals[2])*2], [0, PCn[2][2] * np.sqrt(eigvals[2])*2], color='b')
# center X
X_mean = copy.deepcopy(X_np - np.mean(X_np, axis=0))
# plot data
ax.scatter(X_mean[:, 0], X_mean[:, 1], X_mean[:, 2])
plt.show()
# %%
def transform_data(X: np.array, PCs: np.array, eigenvalues: np.array)-> np.array:
    # center X
    X_mean = copy.deepcopy(X - np.mean(X, axis=0))
    # transform X to the new coordinate system
    return np.dot(X_mean, PCs.T)
# plot transformed data (first three columns) with three principal components as axes
PCn_transformed = transform_data(X_np, PCn, eigvals)
fig = plt.figure()
ax = Axes3D(fig)
# plot eigenvalues as axes
ax.plot([0, np.sqrt(eigvals[0])], [0, 0], [0, 0], color='r')
ax.plot([0, 0], [0, np.sqrt(eigvals[1])], [0, 0], color='g')
ax.plot([0, 0], [0, 0],[0, np.sqrt(eigvals[2])], color='b')
# plot data
ax.scatter(PCn_transformed[:, 0], PCn_transformed[:, 1], PCn_transformed[:, 2])
plt.show()

#### Ans:
Original data with three principal compnents
After transformed

### d
Find all principal components with their eigenvalues using SVD.
#### code:

"""
Find all principal components with their eigenvalues using SVD.
"""
def SVD(X: np.array)-> tuple:
    # Step1: find eigenvalues and eigenvectors
    # I use numpy.linalg.eig to find eigenvalues and eigenvectors because power method can not work on non square matrix and QR decomposition can not grarantee the eigenvectors are correct
    eigenvalues, eigenvectors = np.linalg.eig(np.dot(X.T, X))
    # Step2: obtain S, which is the diagonal matrix whose diagonal entries are the square roots of the eigenvalues of X
    idx = np.argsort(eigenvalues)[::-1]
    e = eigenvalues[idx]
    s = np.sqrt(e)
    S = np.diag(s)
    S1 = np.diag(s)
    for i in range(len(S1)):
        if S1[i][i] != 0:
            S1[i][i] = 1 / S1[i][i]
    # print(s)
    # Step3: obtain V, which is the matrix whose columns are the eigenvectors of X.T * X
    V = eigenvectors[:, idx]
    # Step4: obtain U, which is the matrix whose columns are the eigenvectors of X * X.T
    U = X.dot(V.dot(S1))
    return U, S, V

def find_all_principal_components_svd(X: np.array)-> np.array:
    #  Find all principal components with their eigenvalues using SVD.
    # center X
    X_mean = copy.deepcopy(X - np.mean(X, axis=0))
    # calculate SVD of X
    U, s, V = SVD(X_mean)
    s = s ** 2 / (len(X) - 1)
    return V, s

PCn, eigvals = find_all_principal_components_svd(X_np)

print("Q2.d")
print("Principal components: ")
print(PCn.T)
print("Eigenvalues: ")
print(eigvals)
# %%
U, s, V = SVD(X_np)
print("Principal components: ")
print(V)
print("Eigenvalues: ")
print(s)

#### Ans:

Q2.d
Principal components: 
[[ 0.487  0.57   0.329  0.283  0.148  0.477]
 [ 0.412 -0.111 -0.426 -0.477 -0.469  0.435]
 [ 0.176  0.632 -0.083 -0.443 -0.029 -0.605]
 [-0.599  0.235  0.462 -0.45  -0.221  0.348]
 [-0.297  0.231 -0.573 -0.147  0.642  0.311]
 [-0.339  0.394 -0.403  0.521 -0.544  0.012]]
Eigenvalues: 
[[224.681   0.      0.      0.      0.      0.   ]
 [  0.    119.155   0.      0.      0.      0.   ]
 [  0.      0.      5.439   0.      0.      0.   ]
 [  0.      0.      0.      0.801   0.      0.   ]
 [  0.      0.      0.      0.      0.125   0.   ]
 [  0.      0.      0.      0.      0.      0.   ]]

### e
SVD provides an extra information on U that PCA does not usually have. Is it any interpretation of this U matrix? If yes, please state it
U is the eigenvectors of X * X.T
#### Ans:
C is the covariance matrix, and each matrix can be decomposed through SVD --> $X=USV^T$. Therefore, $U$ will be deleted since it is an orthogonal set. So in this case of PCA it will be eliminated and has no additional use.
$$
C=\frac{X^TX}{n-1}=\frac{VSU^TUSV^T}{n-1}=V\frac{S^2}{n-1}V^T
$$
But if $X$ is regarded as a group of data and if $XX^T$ is needed, $U$ will be its eigenvectors. because:
$$
XX^T = USV^TVSU^T=US^2U^T
$$

### f
Conduct a rank-3 approximation (SVD version) of X
#### code:

PCn, eigvals = find_all_principal_components_svd(X_np)
# get top 3 principal components
PCn = PCn[:3]
# get top 3 eigenvalues
eigvals = eigvals[:3]
# transform data
PCn_transformed = transform_data(X_np, PCn, eigvals)
print("Q2.f")
print("Top 3 principal components: ")
print(PCn_transformed)
print("Top 3 eigenvalues: ")
print(eigvals)

#### Ans:

Q2.f
Top 3 principal components: 
[[  0.697 -19.525   7.44 ]
 [ -3.488   0.106  -2.241]
 [  4.748   9.546  -3.065]
 [  6.381   9.075  -3.129]
 [ -6.457  -2.685   1.068]
 [ -1.881   3.483  -0.074]]
Top 3 eigenvalues: 
[[224.681   0.      0.      0.      0.      0.   ]
 [  0.    119.155   0.      0.      0.      0.   ]
 [  0.      0.      5.439   0.      0.      0.   ]]

## 3
#### code:

rows = range(12)
columns = [1, 5, 9]
X = [[a[i][j] for j in columns] for i in rows]
X_np = np.array(X, dtype=float)
print("Matrix of Q3:")
print(X_np)
print()

#### results:

Matrix of Q3:
[[ 87. 114. 103.]
 [ 81. 109.  97.]
 [ 73.  92.  88.]
 [ 83. 110. 102.]
 [ 92. 109. 111.]
 [ 95. 112. 108.]
 [ 83. 117. 106.]
 [ 78. 115. 104.]
 [ 90. 101. 107.]
 [ 81. 103. 102.]
 [ 87. 109. 103.]
 [ 89. 110. 108.]]

### a
Explain why the monthly data is unlikely to be Gaussian.
#### code:

def guassianaity_of_a_column(X: np.array)-> float:
    # check Guassianaity of X (each column of X) by Kurtosis
    # calculate mean and variance of X
    mean = np.mean(X)
    var = np.var(X)
    # calculate kurtosis
    kurtosis = np.mean((X - mean) ** 4) / var ** 2 - 3
    return kurtosis
print("Q3.a")
for i in range(len(X[0])):
    print("Kurtosis of row {}: {}".format(i, guassianaity_of_a_column(X_np[:,i])))

#### Ans:
All Kurtosis of each row is not close to 0. Forthermore, monthly temperature data typically exhibit seasonal patterns these variations lead to non-Gaussian behavior, as temperatures tend to cluster around certain ranges during specific months (e.g., colder in winter, warmer in summer). For the long term view, climate change and other long-term trends impact temperature data. Non-stationarity (changing statistical properties over time) violates the assumptions of Gaussianity.

Q3.a
Kurtosis of row 0: -0.6224382359992942
Kurtosis of row 1: 0.7578802109307694
Kurtosis of row 2: 1.4866569765083826

### b
Run the three preprocessing steps of ICA on X
#### code

def ICA_preprocess(X: np.array)-> np.array:
    # Run the three preprocessing steps of ICA on X
    # Step 1: Center the data
    X_mean = copy.deepcopy(X - np.mean(X, axis=0))
    # print(X_mean)
    # Step 2: Decorrelate the data
    # calculate covariance matrix of X
    C = np.cov(X_mean.T)
    # find the top n eigenvalues and eigenvectors
    eigenvalues, eigenvectors = deflated_power_method(C, len(C))
    eigenvalues = np.array(eigenvalues)
    eigenvectors = np.array(eigenvectors)
    # sort eigenvalues
    indx = eigenvalues.argsort()
    eigenvalues = eigenvalues[indx]
    eigenvectors = eigenvectors[indx]
    # project X_mean to PCA space
    X_PCA = np.dot(eigenvectors, X_mean.T)
    # Step 3: Whiten the data
    l = np.sqrt(eigenvalues)
    l = 1 / l
    # diagonal matrix
    L = np.diag(l)
    # whiten data
    X_whiten = np.dot(L, X_PCA)
    return X_whiten.T

# X_111 = np.array([[1, 1, 2, 0, 5, 4, 5, 3],
#                   [3, 2, 3, 3, 4, 5, 5, 4]]).T
# X_whiten = ICA_preprocess(X_111)
W_whiten = ICA_preprocess(X_np)
print("Q3.b")
print("Whiten data: ")
print(W_whiten)

#### Ans:

Whiten data: 
[[ 1.283  0.591  0.444]
 [ 1.227  0.745 -0.594]
 [ 0.338 -0.512 -2.669]
 [ 0.104  0.486 -0.096]
 [-0.771 -0.978  0.946]
 [ 1.446 -0.759  1.112]
 [-0.526  1.352  0.583]
 [-1.394  1.683  0.041]
 [-0.795 -1.771  0.089]
 [-1.228 -0.311 -0.638]
 [ 0.718 -0.134  0.142]
 [-0.404 -0.393  0.639]]

### c
Provide a graphical illustration on the transformation, like the one in Lecture 03-2.
#### code

# plot original data
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(X_np[:, 0], X_np[:, 1], X_np[:, 2])
plt.show()

# plot transformed data
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(W_whiten[:, 0], W_whiten[:, 1], W_whiten[:, 2])
plt.show()

#### Ans:
Before:
After:

### d
Run the Fast ICA on Kurtosis Maximization to find the three independent components.
#### code

def fast_ICA(X: np.array,  alpha=1, thresh=1e-8, iterations=5000):
    m, n = X.shape
    # Initialize random weights
    W = np.random.rand(m, m)

    for c in range(m):
            w = W[c, :].copy().reshape(m, 1)
            w = w / np.sqrt((w ** 2).sum())

            i = 0
            lim = 100
            while ((lim > thresh) & (i < iterations)):
                # Dot product of weight and X
                ws = np.dot(w.T, X)
                # Pass w*s into contrast function g
                wg = np.tanh(ws * alpha).T
                # Pass w*s into g prime 
                wg_ = (1 - np.square(np.tanh(ws))) * alpha
                # Update weights
                w_new = (X * wg.T).mean(axis=1) - wg_.mean() * w.squeeze()
                # Decorrelate weights              
                w_new = w_new - np.dot(np.dot(w_new, W[:c].T), W[:c])
                w_new = w_new / np.sqrt((w_new ** 2).sum())
                # Calculate limit condition
                lim = np.abs(np.abs((w_new * w).sum()) - 1)
                # Update weights
                w = w_new
                # Update counter
                i += 1

            W[c, :] = w.T
    return W

X_whiten = ICA_preprocess(X_np).T
W = fast_ICA(X_whiten)
# transform data
X_ICA = np.dot(W, X_whiten)
print("Q3.d")
print("W")
print(W)
print("independent components: ")
print(X_ICA)
# %%
# check independent components are independent
print("Check independent components are independent")
print("Covariance matrix of independent components: ")
print(np.cov(X_ICA))

##### Ans:

Q3.d
W
[[ 0.141 -0.166 -0.976]
 [ 0.698 -0.682  0.217]
 [ 0.702  0.712 -0.019]]
independent components: 
[[-0.35   0.63   2.737  0.028 -0.871 -0.756 -0.868 -0.516  0.094  0.5
  -0.014 -0.616]
 [ 0.589  0.219  0.007 -0.279  0.334  1.769 -1.164 -2.113  0.673 -0.783
   0.624  0.125]
 [ 1.313  1.402 -0.076  0.421 -1.255  0.453  0.583  0.219 -1.82  -1.071
   0.406 -0.575]]
   
Check independent components are independent
Covariance matrix of independent components: 
[[ 1.  0.  0.]
 [ 0.  1. -0.]
 [ 0. -0.  1.]]

## B
determinant is equal to the area of the parallelogram. Please prove it.
#### Ans: