Session 3 - SVD, Differential Equations, Quantum Mechanics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Singular Value Decomposition\n",
    "\n",
    "\n",
    "\n",
    "# Toy example: Artificial Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "%matplotlib notebook\n",
    "plt.rcParams['figure.dpi'] = 150\n",
    "\n",
    "# pca takes centered data matrix\n",
    "def pca(data):\n",
    "    covariance_matrix = (data.transpose() @ data) / data.shape[0]\n",
    "    \n",
    "    eigvals, eigvecs = np.linalg.eig(covariance_matrix)\n",
    "    \n",
    "    sorted_idx = np.argsort(eigvals)[::-1]\n",
    "    \n",
    "    return eigvals[sorted_idx], eigvecs[:,sorted_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of data points\n",
    "N = 100\n",
    "\n",
    "# Data dimension\n",
    "D = 3\n",
    "\n",
    "# Generate data\n",
    "data = np.random.normal(size=(N, D))\n",
    "center = np.mean(data, axis=0)\n",
    "data = data - center\n",
    "\n",
    "eigvals, eigvecs = pca(data)\n",
    "# Eigenvalues are real and positive\n",
    "print(eigvals)\n",
    "\n",
    "# Eigenvectors are orthogonal\n",
    "products = eigvecs.transpose() @ eigvecs\n",
    "print(products)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.axes(projection ='3d')\n",
    "ax.set_xlim(-3, 3)\n",
    "ax.set_ylim(-3, 3)\n",
    "ax.set_zlim(-3, 3)\n",
    "\n",
    "ax.scatter(data[:, 0], data[:, 1], data[:, 2])\n",
    "for v in np.transpose(eigvecs):\n",
    "    ax.plot3D([0, v[0]*3], [0, v[1]*3], [0, v[2]*3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of data points\n",
    "N = 100\n",
    "\n",
    "# Random Embedding Matrix\n",
    "embedding = np.random.normal(size=(2, 3))\n",
    "\n",
    "# Generate and center data\n",
    "data = np.random.normal(size=(N, 2)) @ embedding + np.random.normal(size=(N, 3))*0.1\n",
    "center = np.mean(data, axis=0)\n",
    "data = data - center\n",
    "\n",
    "eigvals, eigvecs = pca(data)\n",
    "print(f'Eigenvalues are {eigvals}')\n",
    "\n",
    "(U, singvals, V) = np.linalg.svd(data)\n",
    "print(f'Singular Values are {singvals}')\n",
    "print(f'Squared Singular Values are {singvals*singvals}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.axes(projection ='3d')\n",
    "ax.set_xlim(-3, 3)\n",
    "ax.set_ylim(-3, 3)\n",
    "ax.set_zlim(-3, 3)\n",
    "\n",
    "ax.scatter(data[:, 0], data[:, 1], data[:, 2])\n",
    "for v in np.transpose(eigvecs):\n",
    "    ax.plot3D([0, v[0]*3], [0, v[1]*3], [0, v[2]*3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.axes(projection ='3d')\n",
    "ax.set_xlim(-3, 3)\n",
    "ax.set_ylim(-3, 3)\n",
    "ax.set_zlim(-3, 3)\n",
    "\n",
    "ax.scatter(data[:, 0], data[:, 1], data[:, 2])\n",
    "for v in np.transpose(eigvecs):\n",
    "    ax.plot3D([0, v[0]*3], [0, v[1]*3], [0, v[2]*3], 'b')\n",
    "for v in np.transpose(V):\n",
    "    ax.plot3D([0, v[0]*3], [0, v[1]*3], [0, v[2]*3], 'r')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform data into eigenvectors coordinate system\n",
    "transformed_data = data @ eigvecs\n",
    "print(f'Variance along first eigenvector: {np.var(transformed_data[:, 0])}')\n",
    "print(f'Variance along second eigenvector: {np.var(transformed_data[:, 1])}')\n",
    "print(f'Veriance along third eigenvector: {np.var(transformed_data[:, 2])}')\n",
    "\n",
    "print(f'Eigenvalues: {eigvals}')\n",
    "\n",
    "transformed_data = data @ V\n",
    "print(f'Variance along first eigenvector: {np.var(transformed_data[:, 0])}')\n",
    "print(f'Variance along second eigenvector: {np.var(transformed_data[:, 1])}')\n",
    "print(f'Veriance along third eigenvector: {np.var(transformed_data[:, 2])}')\n",
    "\n",
    "print(f'Squared Singular values: {singvals*singvals}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Neglect the dimension which corresponds to the lowest eigenvalue\n",
    "transformed_data[:, 2] = 0.\n",
    "\n",
    "# Transform to the original coordinate system\n",
    "transformed_back = transformed_data @ eigvecs.transpose()\n",
    "\n",
    "average_squared_distance = np.sum((data - transformed_back)**2) / N\n",
    "print(f'Average Squared Distance from original data = {average_squared_distance}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EigenFaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('face_data.csv')\n",
    "\n",
    "# 400 Samples, 64x64 images in shades of gray + 1 target label = 4097\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_faces(pixels):\n",
    "    fig, axes = plt.subplots(5, 5, figsize=(6, 6))\n",
    "    for i, ax in enumerate(axes.flat):\n",
    "        ax.imshow(np.array(pixels)[i].reshape(64, 64), cmap='gray')\n",
    "    plt.show()\n",
    "    \n",
    "X = df.drop('target', axis=1)\n",
    "y = df['target']\n",
    "\n",
    "plot_faces(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring the data to the format for our pca function\n",
    "data = np.array(X)\n",
    "\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute SVD.\n",
    "U, singvals, V = np.linalg.svd(data)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot(singvals[:50])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot eigenfaces\n",
    "plot_faces(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transformed_data = data @ V[:25, :].T\n",
    "\n",
    "print(f'Shape of the transformed data is {transformed_data.shape}')\n",
    "\n",
    "transformed_back = transformed_data @ V[:25, :]\n",
    "\n",
    "plot_faces(transformed_back)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Compression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.image as mpimg\n",
    "\n",
    "img = mpimg.imread('cat.png')\n",
    "\n",
    "print(img.shape) \n",
    "\n",
    "fig = plt.figure()\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Treat every column as a vector\n",
    "data = np.reshape(img, (900, 900*4)) \n",
    "print(data.shape) \n",
    "\n",
    "# Center the data\n",
    "center = np.mean(data, axis=0)\n",
    "\n",
    "centered_data = data - center\n",
    "\n",
    "# Compute SVD\n",
    "U, singvals, V = np.linalg.svd(data)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot(singvals[:50])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 50\n",
    "# Project to first n eigenvectors\n",
    "transformed_data = data @ V[:n, :].T\n",
    "\n",
    "print(f'Shape of the transformed data is {transformed_data.shape}')\n",
    "\n",
    "transformed_back = transformed_data @ V[:n, :] \n",
    "\n",
    "# reshape the transformed back data\n",
    "rec_img = np.reshape(transformed_back, (900, 900, 4))\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.imshow(rec_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import resize\n",
    "from skimage.color import rgb2gray, rgba2rgb\n",
    "\n",
    "img_resized = resize(img, (64, 64))\n",
    "img_bw = rgb2gray(rgba2rgb(img_resized))\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "plt.imshow(img_bw, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compressing a cat with EigenFaces?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring the data to the format from eigenfaces \n",
    "data = np.array(X)\n",
    "\n",
    "# Compute SVD.\n",
    "U, singvals, V = np.linalg.svd(data)\n",
    "\n",
    "# Pretend the cat is a face\n",
    "cat = np.reshape(img_bw, (1, 4096))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = np.reshape(img_bw, (1, 4096))\n",
    "\n",
    "fig, axes = plt.subplots(4, 4, figsize=(6, 6))\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    num_vectors = 50 * (i + 1)\n",
    "    reconstruction = data @ V[:num_vectors, :].T @ V[:num_vectors, :]\n",
    "    ax.imshow(reconstruction.reshape(64, 64), cmap='gray')\n",
    "    ax.set_title(f\"#{num_vectors} vectors\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving differential equations\n",
    "\n",
    "## Linear system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a,b,c,d=1, 0,0,1 # Identity\n",
    "# a,b,c,d=-4,1,2,3 # Compression - expansion\n",
    "# a,b,c,d=1,-1,1,1 # Complex eigenvalues\n",
    "# a,b,c,d=1,1,0,1 # Non-diagonalizable\n",
    "transformation = np.array([[a, b], [c, d]])\n",
    "print(\"Coefficients:\", transformation)\n",
    "\n",
    "(vals, vecs) = np.linalg.eig(transformation)\n",
    "print(\"Eigenvalues:\",vals )\n",
    "print(\"Eigenvectors:\",vecs[:, 0], vecs[:, 1] )\n",
    "\n",
    "# Solving the ODE\n",
    "def dX_dt(X, t):\n",
    "    return [a * X[0] + b * X[1], c * X[0] + d * X[1]]\n",
    "\n",
    "ts = np.linspace(0, 0.7, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "for px, py in circle.T:\n",
    "    Ps = odeint(dX_dt, [px, py], ts)\n",
    "    plt.plot(Ps[:,0],  Ps[:,1])\n",
    "\n",
    "r = 5\n",
    "plt.arrow(0,0,r * vals[0] * vecs[0,0], r * vals[0] * vecs[1,0],head_width=1,head_length=2)\n",
    "plt.arrow(0,0,r * vals[1] * vecs[0,1], r * vals[1] * vecs[1,1],head_width=1,head_length=2)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analytical solution\n",
    "\n",
    "def time(P, t):\n",
    "    return vecs @  np.diag(np.exp(vals * t)) @ np.linalg.inv(vecs) @ P\n",
    "\n",
    "fig = plt.figure()\n",
    "for P in circle.T:\n",
    "    Ps = np.array([time(P, t) for t in ts])\n",
    "    plt.plot(Ps[:,0],  Ps[:,1])\n",
    "\n",
    "r = 5\n",
    "plt.arrow(0,0,r * vals[0] * vecs[0,0], r * vals[0] * vecs[1,0],head_width=1,head_length=2)\n",
    "plt.arrow(0,0,r * vals[1] * vecs[0,1], r * vals[1] * vecs[1,1],head_width=1,head_length=2)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analytic solution, second order, for Non-diagonalizable\n",
    "# Only for a,b,c,d=1,1,0,1 # Non-diagonalizable\n",
    "\n",
    "w = np.array([0,1])\n",
    "\n",
    "assert ((transformation - vals[0] * np.identity(2)) @ w == vecs[:, 0]).all()\n",
    "\n",
    "def time(P, t):\n",
    "    return [P[0] * np.exp(t) + P[1] * t * np.exp(t) , P[1] * np.exp(t)]\n",
    "\n",
    "ts = np.linspace(0, 0.7, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "for P in circle.T:\n",
    "    Ps = np.array([time(P, t) for t in ts])\n",
    "    plt.plot(Ps[:,0],  Ps[:,1])\n",
    "\n",
    "r = 5\n",
    "plt.arrow(0,0,r * vals[0] * vecs[0,0], r * vals[0] * vecs[1,0],head_width=1,head_length=2)\n",
    "plt.arrow(0,0,r * vals[1] * vecs[0,1], r * vals[1] * vecs[1,1],head_width=1,head_length=2)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([vecs[:, 0], w])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Eigendecomposition and Quantum Mechanics\n",
    "\n",
    "## Eigendecomposition of a continuous system\n",
    "\n",
    "\n",
    "$$ H = - \\frac{\\hbar^2}{2m} \\nabla^2 + V(x) $$ \n",
    "\n",
    "Trick to solve for $\\nabla^2$\n",
    "\n",
    "$$\\frac{\\partial^2}{\\partial x} \\psi(x) \\simeq \\frac 1 {dx^2} \\left(\\psi(x+dx) - 2*\\psi(x) + \\psi(x-dx)  \\right) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "xmax = 10\n",
    "x = np.linspace(-xmax,xmax,num=n)\n",
    "dx = 2*xmax/(n-1)\n",
    "\n",
    "# Potential\n",
    "V = np.diag(x*x) / 2\n",
    "# Kinetic term\n",
    "T =  -2 * np.identity(n) + np.diag([1]*(n-1), k=1) + np.diag([1]*(n-1), k=-1)\n",
    "T =  - T / 2 / dx**2 \n",
    "# Hamiltonian\n",
    "H = T + V \n",
    "\n",
    "(vals, vecs) = np.linalg.eigh(H)\n",
    "\n",
    "sorted(vals.round(3))[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.plot(x, np.diag(V))\n",
    "for i in range(5):\n",
    "    plt.plot(x, vals[i] + 5 * vecs[:, i])\n",
    "\n",
    "plt.ylim([0, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenbases of a Discrete hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eigs(interaction):\n",
    "    H = np.array([[interaction,1,1,0],[1,0,0,1],[1,0,0,1],[0,1,1,0]])\n",
    "    (vals, vecs) = np.linalg.eigh(H)\n",
    "    return vals\n",
    "    \n",
    "deltas = range(10)\n",
    "l = np.array([eigs(delta) for delta in deltas])\n",
    "    \n",
    "fig = plt.figure()\n",
    "plt.plot(deltas, l)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}