notebooks/pca_from_scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Centralize and normalize $\\mathbf{X}$\n",
    "2. Compute $\\Sigma$, the covariance matrix of $\\mathbf{X}$: $\\Sigma = \\frac{1}{m-1}\\sum_{i=i}^{m} \\left(\\mathbf{x^{(i)}} \\cdot \\mathbf{x^{(i)}}^T\\right)$\n",
    "3. Decompose $\\Sigma$ using SVD: $\\Sigma = U S V^T$. The principal components are the columns of the matrix $U$.\n",
    "4. Select the $k$ columns of $U$ corresponding to the $k$ largest singular values of $S$. These columns are the $k$ principal components.\n",
    "5. Construct the Projection Matrix: $W = \\left[u_1, u_2, \\ldots, u_k\\right]$\n",
    "6. Project the data matrix over the basis formed by $U$: $\\mathbf{X}_{\\text{projected}} = \\mathbf{X}W$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %load ../src/pca.py\n",
    "import numpy as np\n",
    "from scipy import linalg as LA\n",
    "\n",
    "def standardize(X):\n",
    "    # Standardize the data matrix\n",
    "    mean_X = np.mean(X, axis=0)\n",
    "    std_X = np.std(X, axis=0)\n",
    "    standardized_X = (X - mean_X) / std_X\n",
    "    return standardized_X, mean_X, std_X\n",
    "\n",
    "def svd(Sigma):\n",
    "    U, S, V = LA.svd(Sigma)\n",
    "    return U, S, V\n",
    "\n",
    "def pca(X, k):\n",
    "    Sigma = np.cov(X, rowvar = False)\n",
    "    print(f'\\nCovariance matrix:\\n{Sigma}')\n",
    "\n",
    "    U, S, V = svd(Sigma)\n",
    "\n",
    "    print('\\nPercentages of variation:\\n', str(S/np.sum(S) * 100))\n",
    "\n",
    "    # If we wish to reduce dimensionality to k, we first need to compute W\n",
    "    W = U[:, 0:k]\n",
    "    print('\\nW (that is, truncated U):\\n', str(W))\n",
    "\n",
    "    # Now we compute X_projected, the matrix of projected points in k-dimensional space.\n",
    "    X_projected = np.matmul(X, W)\n",
    "    \n",
    "    return X_projected, W\n",
    "\n",
    "def reconstruct_data(X_projected, W):\n",
    "    return np.matmul(X_projected, W.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original Data Matrix:\n",
      "[[3.8700e-01 4.8780e+03 5.4200e+00]\n",
      " [7.2300e-01 1.2104e+04 5.2500e+00]\n",
      " [1.0000e+00 1.2756e+04 5.5200e+00]\n",
      " [1.5240e+00 6.7870e+03 3.9400e+00]]\n",
      "\n",
      "Standardized Data Matrix:\n",
      "[[-1.25237521 -1.2602333   0.60727372]\n",
      " [-0.44547575  0.88082256  0.34085686]\n",
      " [ 0.21973601  1.07400944  0.76398952]\n",
      " [ 1.47811494 -0.6945987  -1.7121201 ]]\n"
     ]
    }
   ],
   "source": [
    "# Data matrix X (four data points, each one is 3D)\n",
    "X = np.array([\n",
    "        [0.387,  4878, 5.42],\n",
    "        [0.723, 12104, 5.25],\n",
    "        [1,     12756, 5.52],\n",
    "        [1.524,  6787, 3.94],\n",
    "    ])\n",
    "\n",
    "# Display the standardized data matrix\n",
    "print(\"Original Data Matrix:\")\n",
    "print(X)\n",
    "\n",
    "standardized_X, _, _ = standardize(X)\n",
    "\n",
    "print(\"\\nStandardized Data Matrix:\")\n",
    "print(standardized_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Covariance matrix:\n",
      "[[ 1.33333333  0.1317339  -1.09173744]\n",
      " [ 0.1317339   1.33333333  0.51489873]\n",
      " [-1.09173744  0.51489873  1.33333333]]\n",
      "U:\n",
      "[[-0.64926351 -0.42018576  0.63395648]\n",
      " [ 0.24566452 -0.90471598 -0.34804876]\n",
      " [ 0.71979569 -0.07023474  0.69062381]]\n",
      "S:\n",
      "[2.49382602 1.43448827 0.07168571]\n",
      "V:\n",
      "[[-0.64926351  0.24566452  0.71979569]\n",
      " [-0.42018576 -0.90471598 -0.07023474]\n",
      " [ 0.63395648 -0.34804876  0.69062381]]\n",
      "\n",
      "Norms of eigenvectors (columns of U):\n",
      "0.9999999999999997 0.9999999999999999 0.9999999999999999\n",
      "dot:\n",
      "\n",
      "-4.0939474033052647e-16\n"
     ]
    }
   ],
   "source": [
    "# Data matrix X (four data points, each one is 3D)\n",
    "X = np.array([\n",
    "        [0.387,  4878, 5.42],\n",
    "        [0.723, 12104, 5.25],\n",
    "        [1,     12756, 5.52],\n",
    "        [1.524,  6787, 3.94],\n",
    "    ])\n",
    "\n",
    "X, _, _ = standardize(X)\n",
    "\n",
    "Sigma = np.cov(X, rowvar = False)\n",
    "print(f'\\nCovariance matrix:\\n{Sigma}')\n",
    "\n",
    "U, S, V = svd(Sigma)\n",
    "\n",
    "print('U:\\n' + str(U))\n",
    "print('S:\\n' + str(S))\n",
    "print('V:\\n' + str(V))\n",
    "\n",
    "# Sanity check: verifying that eivenvectors are indeed unit vectors\n",
    "print('\\nNorms of eigenvectors (columns of U):')\n",
    "print(np.linalg.norm(U[:,0]), np.linalg.norm(U[:,1]), np.linalg.norm(U[:,2]))\n",
    "\n",
    "print('dot:\\n')\n",
    "print(np.dot(U[:,0], U[:,1]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9820785725000001"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "( 2.49382602 + 1.43448827) / (2.49382602 + 1.43448827 + 0.07168571)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "X:\n",
      "[[-1.25237521 -1.2602333   0.60727372]\n",
      " [-0.44547575  0.88082256  0.34085686]\n",
      " [ 0.21973601  1.07400944  0.76398952]\n",
      " [ 1.47811494 -0.6945987  -1.7121201 ]]\n",
      "\n",
      "Covariance matrix:\n",
      "[[ 1.33333333  0.1317339  -1.09173744]\n",
      " [ 0.1317339   1.33333333  0.51489873]\n",
      " [-1.09173744  0.51489873  1.33333333]]\n",
      "\n",
      "Percentages of variation:\n",
      " [62.34565042 35.8622068   1.79214278]\n",
      "\n",
      "W (that is, truncated U):\n",
      " [[-0.64926351 -0.42018576]\n",
      " [ 0.24566452 -0.90471598]\n",
      " [ 0.71979569 -0.07023474]]\n",
      "\n",
      "X_projected (k=2):\n",
      "[[ 0.94063993  1.62373172]\n",
      " [ 0.7509653  -0.63365168]\n",
      " [ 0.6710958  -1.11766206]\n",
      " [-2.36270102  0.12758202]]\n"
     ]
    }
   ],
   "source": [
    "# Data matrix X (four data points, each one is 3D)\n",
    "X = np.array([\n",
    "        [0.387,  4878, 5.42],\n",
    "        [0.723, 12104, 5.25],\n",
    "        [1,     12756, 5.52],\n",
    "        [1.524,  6787, 3.94],\n",
    "    ])\n",
    "\n",
    "X, _, _ = standardize(X)    \n",
    "\n",
    "print('\\nX:\\n' + str(X))\n",
    "\n",
    "# Reduces the dimensionality from 3D to 2D.\n",
    "k=2\n",
    "X_projected, _ = pca(X, k)\n",
    "\n",
    "print('\\nX_projected (k=' + str(k) + '):\\n' + str(X_projected))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "X:\n",
      "[[-1.25237521 -1.2602333   0.60727372]\n",
      " [-0.44547575  0.88082256  0.34085686]\n",
      " [ 0.21973601  1.07400944  0.76398952]\n",
      " [ 1.47811494 -0.6945987  -1.7121201 ]]\n",
      "\n",
      "Covariance matrix:\n",
      "[[ 1.33333333  0.1317339  -1.09173744]\n",
      " [ 0.1317339   1.33333333  0.51489873]\n",
      " [-1.09173744  0.51489873  1.33333333]]\n",
      "\n",
      "Percentages of variation:\n",
      " [62.34565042 35.8622068   1.79214278]\n",
      "\n",
      "W (that is, truncated U):\n",
      " [[-0.64926351 -0.42018576]\n",
      " [ 0.24566452 -0.90471598]\n",
      " [ 0.71979569 -0.07023474]]\n",
      "\n",
      "X_reconstructed:\n",
      "[[-1.29299213 -1.23793419  0.56302618]\n",
      " [-0.22132296  0.75776033  0.58504595]\n",
      " [ 0.03390767  1.17603116  0.56155057]\n",
      " [ 1.48040742 -0.6958573  -1.7096227 ]]\n"
     ]
    }
   ],
   "source": [
    "X = np.array([\n",
    "        [0.387,  4878, 5.42],\n",
    "        [0.723, 12104, 5.25],\n",
    "        [1,     12756, 5.52],\n",
    "        [1.524,  6787, 3.94],\n",
    "    ])\n",
    "\n",
    "X, _, _ = standardize(X)    \n",
    "print('\\nX:\\n' + str(X))\n",
    "\n",
    "X_projected, W = pca(X, k=2)\n",
    "\n",
    "reconstruct_data(W, X_projected)\n",
    "\n",
    "print('\\nX_reconstructed:\\n' + str(np.matmul(X_projected, W.T)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you expect the projected datasets would be for the data matrix X? Explain the results for k=3 and k = 1. Compare with the results obtained in the cell above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Covariance matrix:\n",
      "[[ 9.58333333 19.16666667  0.        ]\n",
      " [19.16666667 38.33333333  0.        ]\n",
      " [ 0.          0.          0.        ]]\n",
      "\n",
      "Percentages of variation:\n",
      " [1.00000000e+02 1.32632086e-14 0.00000000e+00]\n",
      "\n",
      "W (that is, truncated U):\n",
      " [[-0.4472136  -0.89442719  0.        ]\n",
      " [-0.89442719  0.4472136   0.        ]\n",
      " [ 0.          0.          1.        ]]\n",
      "\n",
      "Covariance matrix:\n",
      "[[ 9.58333333 19.16666667  0.        ]\n",
      " [19.16666667 38.33333333  0.        ]\n",
      " [ 0.          0.          0.        ]]\n",
      "\n",
      "Percentages of variation:\n",
      " [1.00000000e+02 1.32632086e-14 0.00000000e+00]\n",
      "\n",
      "W (that is, truncated U):\n",
      " [[-0.4472136 ]\n",
      " [-0.89442719]\n",
      " [ 0.        ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([[ -2.23606798],\n",
       "        [ -4.47213595],\n",
       "        [ -8.94427191],\n",
       "        [-17.88854382]]),\n",
       " array([[-0.4472136 ],\n",
       "        [-0.89442719],\n",
       "        [ 0.        ]]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X1 = np.array([\n",
    "        [1,   2, -.05],\n",
    "        [2,   4, -.05],\n",
    "        [4,   8, -.05],\n",
    "        [8,  16, -.05],\n",
    "    ])\n",
    "pca(X=X1, k=3)\n",
    "pca(X=X1, k=1)"
   ]
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}