{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vector Space\n",
    "\n",
    "A ***vector space*** $V$ defined on the field $K$ is ...\n",
    "\n",
    "It satisfies the 10 following properties:\n",
    "1.\n",
    "\n",
    "A ***vector*** is an element of a vector space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Transformation\n",
    "\n",
    "Mapping between 2 vector spaces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Null and range space\n",
    "\n",
    "Assume two vector spaces $V$ and $W$, and a linear mapping $\\pmb{A}: V \\rightarrow W$.\n",
    "\n",
    "* The null (aka kernel) space of a matrix $\\pmb{A}$ is defined as the subspace of $V$ such that:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{N}(\\pmb{A}) = \\{\\pmb{x} \\in V : \\pmb{Ax} = \\pmb{0}\\} \\subseteq V\n",
    "\\end{equation}\n",
    "\n",
    "* The range (aka column or image) space of a matrix $\\pmb{A}$ is defined as the subspace of $W$ such that:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{R}(\\pmb{A}) = \\{\\pmb{y} \\in W : \\pmb{Ax} = \\pmb{y}, \\; \\forall \\pmb{x} \\in V \\} \\subseteq W\n",
    "\\end{equation}\n",
    "\n",
    "By the ***Rank–nullity theorem***:\n",
    "\\begin{equation}\n",
    "    dim(\\mathcal{N}(\\pmb{A})) + dim(\\mathcal{R}(\\pmb{A})) = dim(V)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rank of a matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transpose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinant\n",
    "\n",
    "Determinant of a square matrix.\n",
    "\n",
    "\\begin{equation}\n",
    "    det : \\mathbb{K}^{N \\times N} \\rightarrow \\mathbb{K}: det(\\pmb{A}) = |\\pmb{A}|  \n",
    "\\end{equation}\n",
    "\n",
    "Properties:\n",
    "* $det(\\pmb{A}) = det(\\pmb{A}^T)$\n",
    "* $det(\\pmb{A}^{-1}) = det(\\pmb{A})^{-1}$\n",
    "* If $\\pmb{A}$ and $\\pmb{B}$ are square matrices of same size then $det(\\pmb{A}\\pmb{B}) = det(\\pmb{A}) det(\\pmb{B}) = det(\\pmb{B}) det(\\pmb{A}) = det(\\pmb{B}\\pmb{A})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Invertible matrix\n",
    "\n",
    "A square matrix $\\pmb{A} \\in \\mathbb{R}^{N \\times N}$ is invertible, and if $\\pmb{A}\\pmb{A}^{-1} = \\pmb{A}^{-1}\\pmb{A} = \\pmb{I}$.\n",
    "\n",
    "* $\\pmb{A}$ is invertible.\n",
    "* $\\pmb{A}$ is full-rank, i.e. \n",
    "* The number 0 is not an eigenvalue of $\\pmb{A}$.\n",
    "\n",
    "Time complexity: $O(N^3)$\n",
    "\n",
    "Notes:\n",
    "* If you want to know if a matrix $\\pmb{A}$ (which is at least PSD) is PD, you can check if it is invertible. If that's the case, then 0 can not be an eigenvalue of $\\pmb{A}$ and thus the matrix is PD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Orthogonal Matrices\n",
    "\n",
    "$\\pmb{A}$ is an orthogonal (square) matrix if $\\pmb{A}^{-1} = \\pmb{A}^T$, and thus $\\pmb{A}^T\\pmb{A} = \\pmb{A}\\pmb{A}^T = \\pmb{I}$. This implies that the columns (resp. rows) of $\\pmb{A}$ are othonormal to each other. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Symmetric Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Positive (Semi) Definite\n",
    "\n",
    "* PSD $\\rightarrow$ symmetric.\n",
    "\n",
    "* All the evals of a PD matrix are positive, and 0 is not one of these.\n",
    "* All the evals of a PSD matrix are non-negative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trace\n",
    "\n",
    "\\begin{equation}\n",
    "    tr(\\pmb{A}) = \\sum_i a_{ii}\n",
    "\\end{equation}\n",
    "\n",
    "Properties:\n",
    "* The trace is a linear operator: $tr(c_1 \\pmb{A}+ c_2 \\pmb{B}) = c_1 tr(\\pmb{A}) + c_2 tr(\\pmb{B})$\n",
    "* The transpose has the same trace: $tr(\\pmb{A}) = tr(\\pmb{A}^T)$\n",
    "* Invariance under cyclic permutation: $tr(\\pmb{A}\\pmb{B}\\pmb{C}) = tr(\\pmb{B}\\pmb{C}\\pmb{A}) = tr(\\pmb{C}\\pmb{A}\\pmb{B})$\n",
    "* Product: $tr(\\pmb{A}\\pmb{B}) \\neq tr(\\pmb{A})tr(\\pmb{B})$\n",
    "* Frobenius Norm: $||\\pmb{A}||_F = \\sqrt{tr(\\pmb{A}^T\\pmb{A})} = \\sqrt{tr(\\pmb{A}\\pmb{A}^T)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Calculus\n",
    "\n",
    "More information can be found on the [matrix cookbook](https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf).\n",
    "\n",
    "First order:\n",
    "\\begin{align}\n",
    "    \\nabla_{\\pmb{X}} tr(\\pmb{AX}) = \\nabla_{\\pmb{X}} tr(\\pmb{XA}) &= \\pmb{A}^T \\\\\n",
    "    \\nabla_{\\pmb{X}} tr(\\pmb{AX}^T) = \\nabla_{\\pmb{X}} tr(\\pmb{X}^T\\pmb{A}) &= \\pmb{A}\n",
    "\\end{align}\n",
    "\n",
    "Second order:\n",
    "\\begin{align}\n",
    "    \\nabla_{\\pmb{X}} tr(\\pmb{A}\\pmb{X}\\pmb{B}\\pmb{X}) = \\nabla_{\\pmb{X}} tr(\\pmb{X}\\pmb{B}\\pmb{X}\\pmb{A}) &= (\\pmb{BXA})^T + (\\pmb{AXB})^T \\\\\n",
    "    \\nabla_{\\pmb{X}} tr(\\pmb{A}\\pmb{X}^T\\pmb{B}\\pmb{X}) = \\nabla_{\\pmb{X}} tr(\\pmb{X}^T\\pmb{B}\\pmb{X}\\pmb{A}) = \\nabla_{\\pmb{X}} tr(\\pmb{B}\\pmb{X}\\pmb{A}\\pmb{X}^T) &= \\pmb{BXA} + (\\pmb{A}\\pmb{X}^T\\pmb{B})^T \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pseudo inverse\n",
    "\n",
    "This is the generalization of the inverse, in the sense that it can be applied to rectangular matrices. Assume $\\pmb{A} \\in \\mathbb{R}^{N \\times N}$.\n",
    "\n",
    "Right pseudo-inverse:\n",
    "\\begin{equation}\n",
    "    \\pmb{A}^\\dagger = (\\pmb{A}^T\\pmb{A})^{-1}\\pmb{A}^T\n",
    "\\end{equation}\n",
    "\n",
    "Left pseudo-inverse:\n",
    "\\begin{equation}\n",
    "    ^\\dagger\\pmb{A} = \\pmb{A}^T(\\pmb{A}\\pmb{A}^T)^{-1}\n",
    "\\end{equation}\n",
    "\n",
    "#### Linear Regression (LR)\n",
    "\n",
    "Assume the input is given by $\\pmb{X} \\in \\mathbb{R}^{N \\times D_x}$, the output by $\\pmb{Y} \\in \\mathbb{R}^{N \\times D_y}$, and the weight matrix by $\\pmb{W} \\in \\mathbb{R}^{D_x \\times D_y}$.\n",
    "\n",
    "The MSE loss is defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{L} = ||\\pmb{Y} - \\pmb{XW}||^2_F\n",
    "\\end{equation}\n",
    "\n",
    "Taking the gradient of this loss and setting it to zero gives us the minimum:\n",
    "\n",
    "\\begin{align}\n",
    "    \\nabla_{\\pmb{W}} \\mathcal{L} &= \\nabla_{\\pmb{W}} ||\\pmb{Y} - \\pmb{XW}||^2_F \\\\\n",
    "    &= \\nabla_{\\pmb{W}} tr((\\pmb{Y} - \\pmb{XW})^T(\\pmb{Y} - \\pmb{XW})) \\\\\n",
    "    &= \\nabla_{\\pmb{W}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr((\\pmb{XW})^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW}) + tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W})] \\\\\n",
    "    &= \\nabla_{\\pmb{W}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW}) - tr(\\pmb{Y}^T\\pmb{XW}) + tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W})] \\\\\n",
    "    &= -2 \\nabla_{\\pmb{W}} tr(\\pmb{Y}^T\\pmb{XW}) + \\nabla_{\\pmb{W}} tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W}) \\\\\n",
    "    &= -2 (\\pmb{Y}^T\\pmb{X})^T + \\pmb{X}^T\\pmb{X} \\pmb{W} + (\\pmb{X}^T\\pmb{X})^T \\pmb{W} \\\\\n",
    "    &= -2 \\pmb{X}^T\\pmb{Y} + 2 \\pmb{X}^T\\pmb{X} \\pmb{W} \\\\\n",
    "    &= 0 \\\\\n",
    "    \\Leftrightarrow & \\quad (\\pmb{X}^T\\pmb{X}) \\pmb{W} = \\pmb{X}^T\\pmb{Y}\n",
    "\\end{align}\n",
    "\n",
    "If the covariance $(\\pmb{X}^T\\pmb{X})$ is invertible i.e. is PD, then the best set of weights are given by:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\pmb{W}^* = (\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T\\pmb{Y} = \\pmb{X}^\\dagger \\pmb{Y} = \\pmb{\\Sigma_{XX}}^{-1}\\pmb{\\Sigma_{XY}}\n",
    "\\end{equation}\n",
    "\n",
    "#### Linear Weighted Regression (LWR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Covariance\n",
    "\n",
    "The covariance $\\pmb{C}_{XY}$ is a positive semi-definite (PSD) matrix which captures linear correlation between 2 random variables $X$ and $Y$. The PSD implies that it is symmetric by definition.\n",
    "\n",
    "* If the cov is invertible, then 0 can not be an eigenvalue of $\\pmb{C}$. This means that it is positive definite (PD)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenvalue and Eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diagonalizable and Eigendecomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Idempotence\n",
    "\n",
    "An ***idempotent*** matrix $\\pmb{P}$ is a square 'matrix which, when multiplied by itself, yields itself'.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\pmb{P} = \\pmb{PP} = \\pmb{P}^2\n",
    "\\end{equation}\n",
    "\n",
    "Properties:\n",
    "* An idempotent matrix (except the identity) is singular (i.e. not full rank).\n",
    "* $\\pmb{I} - \\pmb{P}$ is also idempotent.\n",
    "* An idempotent matrix is always diagonalizable and its eigenvalues are either 0 or 1.\n",
    "* The trace of an idempotent matrix equals the rank of the matrix and thus is always an integer.\n",
    "\n",
    "In linear regression, the optimal solution of $\\mathcal{L} = ||\\pmb{Y} - \\pmb{XW}||^2_F$ with respect to $\\pmb{W}$ is $\\pmb{W}^* = (\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T\\pmb{Y} = \\pmb{X}^\\dagger \\pmb{Y}$.\n",
    "\n",
    "The residual error is then given by:\n",
    "\n",
    "\\begin{equation}\n",
    "    E = (\\pmb{Y} - \\pmb{XW}^*) = (\\pmb{Y} - \\pmb{X}(\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T\\pmb{Y}) = [\\pmb{I} - \\pmb{X}(\\pmb{X}^T\\pmb{X})^{-1} \\pmb{X}^T] \\pmb{Y} = [\\pmb{I} - \\pmb{X}\\pmb{X}^\\dagger] \\pmb{Y} = \\pmb{Q}\\pmb{Y}\n",
    "\\end{equation}\n",
    "\n",
    "The matrices $\\pmb{P} = \\pmb{X}\\pmb{X}^\\dagger$ and $\\pmb{Q} = [\\pmb{I} - \\pmb{X}\\pmb{X}^\\dagger] = [\\pmb{I} - \\pmb{P}]$ are symmetric and idempotent matrices.\n",
    "\n",
    "* An idempotent linear operator $\\pmb{P}$ is a projection operator on the range space $\\mathcal{R}(\\pmb{P})$ along its null space $\\mathcal{N}(\\pmb{P})$.\n",
    "* $\\pmb{P}$ is an orthogonal projection operator $\\Leftrightarrow$ it is idempotent and symmetric.\n",
    "Ex: $\\pmb{P} = \\left[\\begin{array}{cc} \\pmb{I} & \\pmb{0} \\\\ \\pmb{0} & \\pmb{0} \\end{array} \\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Orthogonal Complement and Projection Matrix\n",
    "\n",
    "\"A ***projection*** is a linear transformation $\\pmb{P}$ from a vector space $V$ to itself such that $\\pmb{P}^2 = \\pmb{P}$ (i.e. $\\pmb{P}$ is idempotent). That is, whenever $\\pmb{P}$ is applied twice to any value, it gives the same result as if it were applied once (idempotent). It leaves its image unchanged.\n",
    "\n",
    "Let $V$ be a finite dimensional vector space and $\\pmb{P}$ be a projection on $V$. Suppose the subspaces $\\mathcal{R}$ and $\\mathcal{N}$ are the range and null space (aka kernel) of $\\pmb{P}$ respectively. Then $\\pmb{P}$ has the following properties:\n",
    "1. $\\pmb{P}$ is idempotent by def (i.e. $\\pmb{P}^2 = \\pmb{P}$)\n",
    "2. $\\pmb{P}$ is the identity operator $\\pmb{I}$ on $\\mathcal{R}$ (i.e. $\\forall \\pmb{x} \\in \\mathcal{R}: \\pmb{Px} = \\pmb{x}$)\n",
    "3. We have a direct sum $V = \\mathcal{R} \\oplus \\mathcal{N}$. Every vector $\\pmb{x} \\in V$ may be decomposed uniquely as $\\pmb{x} = \\pmb{r} + \\pmb{n}$ with $\\pmb{r} = \\pmb{Px} \\in \\mathcal{R}$ and $\\pmb{n} = \\pmb{x} − \\pmb{Px} = (\\pmb{I} - \\pmb{P}) \\pmb{x} \\in \\mathcal{N}$.\n",
    "\n",
    "The range and kernel of a projection are complementary, as are $\\pmb{P}$ and $\\pmb{Q} = \\pmb{I} − \\pmb{P}$. The operator $\\pmb{Q}$ is also a projection, and the range and null spaces of $\\pmb{P}$ become the null and range spaces of $\\pmb{Q}$ and vice versa. We say $\\pmb{P}$ is a projection along $\\mathcal{N}$ onto $\\mathcal{R}$, and $\\pmb{Q}$ is a projection along $\\mathcal{R}$ onto $\\mathcal{N}$.\"\n",
    "\n",
    "##### Orthogonal Projection\n",
    "\n",
    "When the vector space $V$ has an inner product and is complete (i.e. it is a Hilbert space), the concept of orthogonality can be used. An ***orthogonal projection*** is a projection for which the range $\\mathcal{R}$ and the null space $\\mathcal{N}$ are orthogonal subspaces. That is, $\\forall \\pmb{r} \\in \\mathcal{R}, \\forall \\pmb{n} \\in \\mathcal{N}: \\pmb{r} . \\pmb{n} = 0$.\n",
    "\n",
    "The ***orthogonal complement*** of a subspace $W$ of a vector space $V$ equipped with a bilinear form $B$ is the set $W^\\perp$ of all vectors in $V$ that are orthogonal to every vector in $W$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SVD\n",
    "\n",
    "Any rectangular matrices $\\pmb{A} \\in \\mathbb{R}^{N \\times M}$ can be decomposed into a product of 3 matrices (which can be seen as the building blocks of $\\pmb{A}$):\n",
    "\n",
    "\\begin{equation}\n",
    "    \\pmb{A} = \\pmb{U_A}\\pmb{\\Sigma_A}\\pmb{V_A}^T\n",
    "\\end{equation}\n",
    "where:\n",
    "* $\\pmb{U_A} \\in \\mathbb{R}^{N \\times N}$ is an orthogonal matrix (i.e. $\\pmb{U_A}^T = \\pmb{U_A}^{-1}$ thus $\\pmb{U_A}\\pmb{U_A}^T = \\pmb{U_A}^T\\pmb{U_A} = \\pmb{I}$). The columns of $\\pmb{U_A}$ contains the eigenvectors of the PSD (thus symmetric) $\\pmb{A}\\pmb{A}^T$, and are also known as the left-singular vectors of $\\pmb{A}$. Its first columns associated with non-zero singular values span the range of $\\pmb{A}$. Note that because $\\pmb{U_A}$ is an orthogonal matrix, it is invertible and thus has full rank, i.e. its columns form a basis and span $\\mathbb{R}^N$.\n",
    "* $\\pmb{\\Sigma_A} \\in \\mathbb{R}^{N \\times M}$ is a rectangular matrix. The upper left submatrix is a diagonal matrix of size $r \\times r$ with $r = \\min(N,M)$ while the rest is filled with zeros. The diagonal elements are the singular values ordered by descending order. The singular values (SVs) $\\sigma_i$ are equals to the square roots of the eigenvalues $\\lambda_i$ of $\\pmb{A}\\pmb{A}^T$ and $\\pmb{A}^T\\pmb{A}$. The rank of $\\pmb{A}$ is given by the number of SVs different from 0.\n",
    "* $\\pmb{V_A} \\in \\mathbb{R}^{M \\times M}$ is an orthogonal matrix (i.e. $\\pmb{V_A}^T = \\pmb{V_A}^{-1}$ thus $\\pmb{V_A}\\pmb{V_A}^T = \\pmb{V_A}^T\\pmb{V_A} = \\pmb{I}$). The columns of $\\pmb{V_A}$ contains the eigenvectors of the PSD (thus symmetric) $\\pmb{A}^T\\pmb{A}$, and are also known as the right-singular vectors of $\\pmb{A}$. Its last columns associated with vanishing singular values $\\sigma_i = 0$ span the null space of $\\pmb{A}$. Note that because $\\pmb{V_A}$ is an orthogonal matrix, it is invertible and thus has full rank, i.e. its columns form a basis and span $\\mathbb{R}^M$.\n",
    "\n",
    "Few notes:\n",
    "* SVD can be seen as a generalization of eigendecomposition.\n",
    "* SVD can be compressed such that $\\pmb{A} = \\pmb{\\tilde{U}_A}\\pmb{\\tilde{\\Sigma}_A}\\pmb{\\tilde{V_A}^T}$ with $\\pmb{\\tilde{U}_A} \\in \\mathbb{R}^{N \\times r}$, $\\pmb{\\tilde{\\Sigma}_A} \\in \\mathbb{R}^{r \\times r}$, $\\pmb{V_A} \\in \\mathbb{R}^{M \\times r}$, and $r=\\min(N,M)$.\n",
    "* \n",
    "\n",
    "Some properties:\n",
    "* **Existence**:\n",
    "* **Uniqueness**:\n",
    "* **Transpose**: $\\pmb{A}^T = (\\pmb{U_A}\\pmb{\\Sigma_A}\\pmb{V_A}^T)^T = \\pmb{V_A}\\pmb{\\Sigma_A}^T\\pmb{U_A}^T $\n",
    "* If the matrix A is symmetric, then it has real eigenvalues.\n",
    "* **Pseudo-inverse**: The pseudo-inverse of $\\pmb{A}$ is given by $\\pmb{A^\\dagger} = \\pmb{V_A}\\pmb{\\Sigma_A}^{-1}\\pmb{U_A}^T$ with $\\pmb{\\Sigma_A}^{-1}$ containing the inverse of singular values on its diagonal.\n",
    "* **Null and range space** of $\\pmb{A}$: $\\mathcal{N}(\\pmb{A}) \\equiv \\mathcal{R}^\\perp(\\pmb{A}^T)$ and $\\mathcal{R}(\\pmb{A}) \\equiv \\mathcal{N}^\\perp(\\pmb{A}^T)$. This can be seen by applying SVD on $\\pmb{A}$ and $\\pmb{A}^T$. The last columns of $\\pmb{V_A}$ associated with vanishing SVs spans the null-space of $\\pmb{A}$, and the first columns of $\\pmb{V_A}$ associated with non-zero SVs spans the range space of $\\pmb{A}^T$. Because of the orthogonality property of $\\pmb{V_A}$, we have the $r$ first columns are orthogonals to the $M-r$ last columns. $dim(\\mathcal{R}(\\pmb{A})) + dim(\\mathcal{N}(\\pmb{A})) = M$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Block Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA\n",
    "\n",
    "The (MSE) loss minimized by PCA is:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{L} = ||\\pmb{X} - \\pmb{XWW}^T||^2_F\n",
    "\\end{equation}\n",
    "\n",
    "where is the $\\pmb{\\tilde{X}} = \\pmb{XW}$ is the projected data on the lower dimensional space, and $\\pmb{\\hat{X}} = \\pmb{\\tilde{X}W}^T$ is the data projected back to the original space.\n",
    "\n",
    "Taking the gradient of this loss with respect to $\\pmb{W}$ and setting it to $0$ gives us the minimum:\n",
    "\n",
    "\\begin{align}\n",
    "    \\nabla_{\\pmb{W}} \\mathcal{L} &= \\nabla_{\\pmb{W}} tr((\\pmb{X} - \\pmb{XWW}^T)^T(\\pmb{X} - \\pmb{XWW}^T)) \\\\\n",
    "    &= \\nabla_{\\pmb{W}} tr((\\pmb{X}^T - \\pmb{WW}^T \\pmb{X}^T) (\\pmb{X} - \\pmb{XWW}^T)) \\\\\n",
    "    &= \\nabla_{\\pmb{W}} [tr(\\pmb{X}^T\\pmb{X}) - tr(\\pmb{WW}^T \\pmb{X}^T\\pmb{X}) - tr(\\pmb{X}^T\\pmb{X}\\pmb{WW}^T) + tr(\\pmb{WW}^T\\pmb{X}^T\\pmb{XWW}^T)] \\\\\n",
    "    &= \\nabla_{\\pmb{W}} [- 2 tr(\\pmb{W}^T \\pmb{\\Sigma_{XX}} \\pmb{W}) + tr(\\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{WW}^T)] \\\\\n",
    "    &= - 2 (\\pmb{\\Sigma_{XX}} \\pmb{W} + \\pmb{\\Sigma_{XX}}^T \\pmb{W}) + (\\pmb{W}^T \\pmb{\\Sigma_{XX}} \\pmb{WW}^T)^T + (\\pmb{\\Sigma_{XX}} \\pmb{WW}^T)\\pmb{W} + (\\pmb{WW}^T \\pmb{\\Sigma_{XX}})^T (\\pmb{W}^T)^T + (\\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W}) \\\\\n",
    "    &= - 4 \\pmb{\\Sigma_{XX}} \\pmb{W} + \\pmb{WW}^T \\pmb{\\Sigma_{XX}}^T \\pmb{W} + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T\\pmb{W} + \\pmb{\\Sigma_{XX}}^T\\pmb{WW}^T \\pmb{W} + \\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W} \\\\\n",
    "    &= - 4 \\pmb{\\Sigma_{XX}} \\pmb{W} + 2 \\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W} + 2 \\pmb{\\Sigma_{XX}} \\pmb{WW}^T\\pmb{W} \\\\\n",
    "    &= 0 \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{WW}^T \\pmb{\\Sigma_{XX}} \\pmb{W} + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T\\pmb{W} = 2  \\pmb{\\Sigma_{XX}} \\pmb{W} \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{WW}^T \\pmb{\\Sigma_{XX}} + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T = 2 \\pmb{\\Sigma_{XX}} \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{WW}^T + \\pmb{\\Sigma_{XX}} \\pmb{WW}^T \\pmb{\\Sigma_{XX}}^{-1} = 2 \\pmb{I} \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{\\Sigma_{XX}} = (2 \\pmb{I} - \\pmb{WW}^T) \\pmb{\\Sigma_{XX}} (\\pmb{WW}^T)^{-1} \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{Q\\Lambda Q}^T = (2 \\pmb{I} - \\pmb{WW}^T) \\pmb{Q\\Lambda Q}^T (\\pmb{WW}^T)^{-1} \\\\\n",
    "    \\Leftrightarrow & \\quad \\left\\{ \\begin{array}{l} (2 \\pmb{I} - \\pmb{WW}^T) \\pmb{Q} = \\pmb{Q} \\\\ \\pmb{Q}^T (\\pmb{WW}^T)^{-1} = \\pmb{Q}^T \\: \\Leftrightarrow \\: \\pmb{QQ}^T = \\pmb{WW}^T \\end{array} \\right.  \\\\\n",
    "\\end{align}\n",
    "\n",
    "Algo using the covariance:\n",
    "1. subtract the mean $\\pmb{X}$\n",
    "2. compute the covariance matrix $\\pmb{\\Sigma_{XX}} = \\pmb{X}^T\\pmb{X}$\n",
    "3. compute the eigendecomposition of $\\pmb{\\Sigma_{XX}}$, i.e. compute $\\pmb{Q}$ and $\\pmb{\\Lambda}$ such that $\\pmb{\\Sigma_{XX}} = \\pmb{Q} \\pmb{\\Lambda} \\pmb{Q}^T$.\n",
    "4. return the sorted evals and the corresponding evecs\n",
    "\n",
    "Algo using SVD:\n",
    "1. substract the mean from $\\pmb{X}$\n",
    "2. compute SVD of $\\pmb{X}$, i.e. $\\pmb{X} = \\pmb{U_X \\Sigma_X V_X}^T$. The eigenvectors are given by $\\pmb{Q} = \\pmb{V_X}$ and the evals by $\\pmb{\\Lambda} = \\pmb{\\Sigma_X}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recursive PCA/SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hierarchical PCA/SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear FNN with 1 hidden layer\n",
    "\n",
    "Assume a linear Feedforward Neural Network (FNN) with 1 input, 1 hidden, and 1 output layer. Assume the input data is given by $\\pmb{X} \\in \\mathbb{R}^{N \\times D_x}$, the output data by $\\pmb{Y} \\in \\mathbb{R}^{N \\times D_y}$, the hidden data by $\\pmb{H} \\in \\mathbb{R}^{N \\times D_h}$, the input-hidden weight matrix by $\\pmb{W_1} \\in \\mathbb{R}^{D_x \\times D_h}$ and the hidden-output weight matrix $\\pmb{W_2} \\in \\mathbb{R}^{D_h \\times D_y}$. These variables are related by the following relationships:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\pmb{H} = \\pmb{X} \\pmb{W_1} \\qquad \\mbox{and} \\qquad \\pmb{Y} = \\pmb{H} \\pmb{W_2}\n",
    "\\end{equation}\n",
    "\n",
    "The MSE loss is thus defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{L} = ||\\pmb{Y} - \\pmb{XW_1W_2}||^2_F\n",
    "\\end{equation}\n",
    "\n",
    "Taking the gradient of this loss with respect to $\\pmb{W_1}$ and $\\pmb{W_2}$, and setting these to zero gives us the minimum:\n",
    "\n",
    "\\begin{align}\n",
    "    \\nabla_{\\pmb{W_1}} \\mathcal{L} &= \\nabla_{\\pmb{W_1}} ||\\pmb{Y} - \\pmb{XW_1W_2}||^2_F \\\\\n",
    "    &= \\nabla_{\\pmb{W_1}} tr((\\pmb{Y} - \\pmb{XW_1W_2})^T(\\pmb{Y} - \\pmb{XW_1W_2})) \\\\\n",
    "    &= \\nabla_{\\pmb{W_1}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr((\\pmb{XW_1W_2})^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW_1W_2}) + tr(\\pmb{W}^T\\pmb{X}^T\\pmb{X}\\pmb{W})] \\\\\n",
    "    &= \\nabla_{\\pmb{W_1}} [tr(\\pmb{Y}^T\\pmb{Y}) - tr(\\pmb{Y}^T\\pmb{XW_1W_2}) - tr(\\pmb{Y}^T\\pmb{XW_1W_2}) + tr(\\pmb{(W_1W_2)}^T\\pmb{X}^T\\pmb{X}\\pmb{(W_1W_2)})] \\\\\n",
    "    &= \\nabla_{\\pmb{W_1}} [-2 tr(\\pmb{Y}^T\\pmb{XW_1W_2}) + tr(\\pmb{W_2}^T\\pmb{W_1}^T\\pmb{X}^T\\pmb{X}\\pmb{W_1}\\pmb{W_2})] \\\\\n",
    "    &= \\nabla_{\\pmb{W_1}} [-2 tr(\\pmb{\\Sigma_{YX}}\\pmb{W_1W_2}) + tr(\\pmb{W_2}^T\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2})] \\\\\n",
    "    &= -2 (\\pmb{W_2}\\pmb{\\Sigma_{YX}})^T + \\pmb{\\Sigma_{XX}}\\pmb{W_1W_2W_2}^T + \\pmb{\\Sigma_{XX}}^T\\pmb{W_1}(\\pmb{W_2W_2}^T)^T \\\\\n",
    "    &= -2 \\pmb{\\Sigma_{XY}} \\pmb{W_2}^T + 2 \\pmb{\\Sigma_{XX}}\\pmb{W_1W_2W_2}^T \\\\\n",
    "    &= 0 \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{\\Sigma_{XX}}\\pmb{W_1W_2W_2}^T = \\pmb{\\Sigma_{XY}} \\pmb{W_2}^T \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{W_1} = \\pmb{\\Sigma_{XX}}^{-1} \\pmb{\\Sigma_{XY}} \\pmb{W_2}^T (\\pmb{W_2W_2}^T)^{-1}\n",
    "\\end{align}\n",
    "\n",
    "and \n",
    "\n",
    "\\begin{align}\n",
    "    \\nabla_{\\pmb{W_2}} \\mathcal{L} &= \\nabla_{\\pmb{W_2}} [-2 tr(\\pmb{\\Sigma_{YX}}\\pmb{W_1W_2}) + tr(\\pmb{W_2}^T\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2})] \\\\\n",
    "    &= -2 (\\pmb{\\Sigma_{YX}}\\pmb{W_1})^T + \\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2} + (\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1})^T\\pmb{W_2} \\\\\n",
    "    &= -2 \\pmb{W_1}^T\\pmb{\\Sigma_{XY}} + 2 \\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2} \\\\\n",
    "    &= 0 \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1}\\pmb{W_2} = \\pmb{W_1}^T \\pmb{\\Sigma_{XY}} \\\\\n",
    "    \\Leftrightarrow & \\quad \\pmb{W_2} = (\\pmb{W_1}^T\\pmb{\\Sigma_{XX}}\\pmb{W_1})^{-1} \\pmb{W_1}^T \\pmb{\\Sigma_{XY}} \\\\\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}