<jupyter_start><jupyter_text>Principal Component Analysis (PCA)  Imports<jupyter_code># %matplotlib notebook
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Load the MNIST dataset.
data = sio.loadmat('mnist_digits.mat')
X = data['x']
Y = data['y']
n, dim = X.shape
print(f'The dataset has {n} examples of dimension {dim}.')

# Select only the digits in `digits`.
digits = [1, 9]  # 1 and 9 are cherry picked.
# digits = list(range(10))
indices = []
for digit in range(10):
    indices.append(np.argwhere(Y == digit)[:, 0])

X = np.concatenate([X[indices[digit]] for digit in digits], axis=0)
Y = np.concatenate([Y[indices[digit]] for digit in digits], axis=0)

indices = []
for digit in range(10):
    indices.append(np.argwhere(Y == digit)[:, 0])

n, dim = X.shape
print(f'The trimmed down dataset has {n} examples of dimension {dim}.')

# Add noise to the data.
sigma = 0.1
noisy_X = X + sigma * np.random.randn(*X.shape)<jupyter_output>The dataset has 66000 examples of dimension 784.
The trimmed down dataset has 13998 examples of dimension 784.<jupyter_text>PCA<jupyter_code>def pca(X, d=1, plot=False, label=''):
    # Assert that the data is centered.
    mean = np.mean(X)
    X_centered = np.array([x - mean for x in X])
    
    # Compute the eigenvectors and eigenvalues.
    eig_val, eig_vect = np.linalg.eig(X_centered.T @ X_centered)
    # Remove the imaginary part.
    eig_val = np.real(eig_val)
    eig_vect = np.real(eig_vect)
    
    # Sort by eigen value.
    sort_idx = np.argsort(eig_val)
    sorted_eig_val = eig_val[sort_idx]
    sorted_eig_vect = eig_vect[sort_idx]
    
    # Plot the eigenvalues.
    if plot:
        plt.plot(np.maximum(sorted_eig_val, 10 ** (-11)), label=label)
        plt.yscale('log')

    # Compute the projection on the subspace.
    sub_space = sorted_eig_vect[-d:]
    Z = np.array([np.array([np.dot(v, x) for v in sub_space]) for x in X_centered])
    
    return np.array([z + mean for z in Z])<jupyter_output><empty_output><jupyter_text>1D<jupyter_code>Z = np.sort(pca(X, plot=True, label='not noisy'), axis=None)
noisy_Z = np.sort(pca(noisy_X, plot=True, label='noisy'), axis=None)
plt.title('MNIST PCA 1D')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>2D<jupyter_code>def show_2d_pca(X, nb_points=100):
    Z = pca(X, d=2)
    for digit in digits:
        Z_digit = Z[indices[digit]][:nb_points]
        plt.scatter(Z_digit[:, 0], Z_digit[:, 1], label=f'{digit}')

    plt.title('MNIST PCA 2D')
    plt.legend()
    plt.show()
show_2d_pca(X)
show_2d_pca(noisy_X)<jupyter_output><empty_output><jupyter_text>3D<jupyter_code>def show_3d_pca(X, nb_points=100):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    Z = pca(X, d=3)
    for digit in digits:
        Z_digit = Z[indices[digit]][:nb_points]
        ax.scatter(Z_digit[:, 0], Z_digit[:, 1], Z_digit[:, 2], label=f'{digit}')

    plt.title('MNIST PCA 3D')
    plt.legend()
    plt.show()
show_3d_pca(X)
show_3d_pca(noisy_X)<jupyter_output><empty_output><jupyter_text>OLS<jupyter_code>def ols(X):
    X_powers = np.array([np.array([1, *x]) for x in X], dtype=np.float64)
    theta = np.linalg.lstsq(X, Y, rcond=None)[0]
    return theta
d = 30

# Not noisy.
Z = pca(X, d=d)

theta_ref = ols(X)
theta_pca = ols(Z)

# Noisy.
noisy_Z = pca(noisy_X, d=d)

noisy_theta_ref = ols(noisy_X)
noisy_theta_pca = ols(noisy_Z)
Y_pca = np.array([np.dot(theta_pca.T, z) for z in Z])
Y_ref = np.array([np.dot(theta_ref.T, x) for x in X])

noisy_Y_pca = np.array([np.dot(noisy_theta_pca.T, z) for z in noisy_Z])
noisy_Y_ref = np.array([np.dot(noisy_theta_ref.T, x) for x in noisy_X])
def show_Y(Y, digits, title=''):
    for digit in digits:
        ys = np.sort(Y[indices[digit]], axis=None)
        plt.plot(ys, label=f'{digit}')
        
    plt.title(title)
    plt.legend()
    plt.show()
show_Y(Y_ref, digits, 'ref')
show_Y(Y_pca, digits, 'pca')
show_Y(Y_ref, digits, 'noisy ref')
show_Y(Y_pca, digits, 'noisy pca')
def error(Y, y1, y2, threshold):
    n = len(indices[y1]) + len(indices[y2])
    return (np.sum(Y[indices[y1]] >= threshold) + np.sum(Y[indices[y2]] < threshold)) / n
y1 = 1
y2 = 9<jupyter_output><empty_output><jupyter_text>Optimize the thresholds by linear search<jupyter_code>Ys = [('ref', Y_ref), ('pca', Y_pca), ('noisy ref', noisy_Y_ref), ('noisy pca', noisy_Y_pca)]

print('error\tthreshold\tname\n' + '-' * 28)
for Y_name, Y in Ys:
    # -5 and 25 shouldn't be hardcoded :(
    e, t = min([(error(Y, y1, y2, threshold), threshold) for threshold in np.linspace(-5, 25, 10 ** 5)])
    print(f'{e * 100:.2f}%\t{t:.2f}\t\t{Y_name}')<jupyter_output>error	threshold	name
----------------------------
0.72%	3.92		ref
4.79%	4.29		pca
0.95%	3.71		noisy ref
3.16%	4.42		noisy pca
{'fraction_non_alphanumeric': 0.1637195756188891, 'fraction_numerical': 0.026211774495527355, 'mean_word_length': 3.9720785935884177, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '333515', 'license': 'MIT', 'n_tokens_mistral': 2149, 'n_tokens_neox': 2015, 'n_words': 452, 'path': 'PCA.ipynb', 'repo_name': 'gabriel-doriath-dohler/ML-experiments', 'size': '4807'}