EOF_simple.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple EOF analysis based on SVD\n",
    "A is  3 $\\times$ 4 matrix. (location, time) or (n_features, n_samples)\n",
    "\n",
    "Singular value decomposition\n",
    "$A\\ = U\\ \\Sigma\\ V^T$\n",
    "\n",
    "$U$ (3 $\\times$ 3) = EOF (location, mode)   \n",
    "$\\Sigma$ (3 $\\times$ 4)  square roots of the nonzero eigenvalues  \n",
    "$V$ (4 $\\times$ 4)\n",
    "\n",
    "$ PC = \\Sigma\\ V^T$ (3 $\\times$ 4) (mode, time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ -3.          -4.66666667  -1.          -5.        ]\n",
      " [  3.          16.33333333  -1.           3.        ]\n",
      " [  0.         -11.66666667   2.           2.        ]] (3, 4)\n"
     ]
    }
   ],
   "source": [
    "a = np.array([[2, 4, 1, -6], [8, 25, 1, 2], [5, -3, 4, 1]])\n",
    "mu = a.mean(axis=0)\n",
    "\n",
    "# Remove sample mean\n",
    "a = a - mu\n",
    "\n",
    "print(a, a.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "eof, sigma, VT = np.linalg.svd(a, full_matrices=True)\n",
    "smat = np.zeros(a.shape)\n",
    "smat[:len(sigma), :len(sigma)] = np.diag(sigma)\n",
    "\n",
    "pc = smat @ VT\n",
    "pca_score = sigma**2 / (sigma**2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EOF1 = [ 0.26447791 -0.80122219  0.53674429]\n",
      "EOF2 = [-0.77247531  0.15719307  0.61528224]\n",
      "dot product of EOF1 and EOF2 = 0.000\n"
     ]
    }
   ],
   "source": [
    "print('EOF1 =', eof[:, 0])\n",
    "print('EOF2 =', eof[:, 1])\n",
    "print('dot product of EOF1 and EOF2 = {:.3f}'.format(eof[:, 0] @ eof[:, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstruction\n",
    "\n",
    "We can obtain the original matrix by multiplying eigen vectors and princple components!!\n",
    "\n",
    "X = EZ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.,  4.,  1., -6.],\n",
       "       [ 8., 25.,  1.,  2.],\n",
       "       [ 5., -3.,  4.,  1.]])"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eof @ pc + mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## check $E^T E = I$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00  3.37786034e-18 -2.33237467e-16]\n",
      " [ 3.37786034e-18  1.00000000e+00 -5.16664735e-17]\n",
      " [-2.33237467e-16 -5.16664735e-17  1.00000000e+00]]\n",
      "3.000000000000001 True\n"
     ]
    }
   ],
   "source": [
    "I = eof.T @ eof\n",
    "print(I)\n",
    "print(np.trace(I), np.allclose(I, np.eye(I.shape[0]), atol=1e-06))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projection\n",
    "\n",
    "Project arbitary x onto EOF space  \n",
    "Z = E^T X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 1)"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#x = np.array([[2, 4], [3, 3], [2, 6]])\n",
    "x = np.array([[1], [3], [2]])\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0657001 ],\n",
       "       [ 0.92966838],\n",
       "       [ 3.46410162]])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# projection coefficients\n",
    "z = eof.T @ x\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.],\n",
       "       [3.],\n",
       "       [2.]])"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eof @ z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check if the projection results can be obtained from the least squre fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0657001 ],\n",
       "       [ 0.92966838],\n",
       "       [ 3.46410162]])"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z1, res, rank, s = np.linalg.lstsq(eof, x, rcond=None)\n",
    "z1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The coefficients do not change even in the truncated EOFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0657001 ],\n",
       "       [ 0.92966838]])"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eof_tr = eof[:, :2]\n",
    "z1_tr, res, rank, s = np.linalg.lstsq(eof_tr, x, rcond=None)\n",
    "z1_tr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.00000000e+00],\n",
       "       [ 1.00000000e+00],\n",
       "       [-1.55431223e-15]])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eof_tr @ eof_tr.T @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.66666667, -0.33333333, -0.33333333],\n",
       "       [-0.33333333,  0.66666667, -0.33333333],\n",
       "       [-0.33333333, -0.33333333,  0.66666667]])"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EET = eof_tr @ eof_tr.T\n",
    "EET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9.00719925e+15, 9.00719925e+15, 9.00719925e+15],\n",
       "       [9.00719925e+15, 9.00719925e+15, 9.00719925e+15],\n",
       "       [9.00719925e+15, 9.00719925e+15, 9.00719925e+15]])"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(eof_tr @ eof_tr.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[5.40431955e+16],\n",
       "       [5.40431955e+16],\n",
       "       [5.40431955e+16]])"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_im = np.linalg.inv(eof_tr @ eof_tr.T)@ x\n",
    "x_im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5.33333333],\n",
       "       [-6.66666667],\n",
       "       [-0.66666667]])"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eof_tr @ eof_tr.T @ x_im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.0049636 , -0.53453177, -1.47043183],\n",
       "       [ 0.52886841,  0.9144099 , -1.44327831],\n",
       "       [ 0.26447791, -0.80122219,  0.53674429]])"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(eof_tr @ eof_tr.T) @ eof_tr @ eof_tr.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.66666667, -0.33333333, -0.33333333],\n",
       "       [-0.33333333,  0.66666667, -0.33333333],\n",
       "       [-0.33333333, -0.33333333,  0.66666667]])"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.pinv(EET) @ EET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00,  3.33333333e-01,  5.00000000e-01],\n",
       "       [-5.00000000e-01,  3.33333333e-01, -5.00000000e-01],\n",
       "       [ 2.33146835e-15, -6.66666667e-01,  0.00000000e+00]])"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EET @ np.linalg.inv(EET)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.7007434154172e-17"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(EET)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00, -9.28077013e-17, -7.56324913e-17],\n",
       "       [-9.28077013e-17,  1.00000000e+00, -2.80501886e-16],\n",
       "       [-7.56324913e-17, -2.80501886e-16,  1.00000000e+00]])"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eof @ eof.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}