lectures/overview_part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Total recall of the first part of NLA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The main topics that were covered\n",
    "\n",
    "- Floating point arithmetic and related issues\n",
    "- Main operations with matrices and vectors, their complexity and special cases\n",
    "- Matrix decompositions: SVD, skeleton decompositions, LU/Cholecky decomposition, QR, Schur decomposition\n",
    "- SVD applications\n",
    "- Computing eigenvalues and eigenvectors: QR algorithm and power method\n",
    "- How we compute SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular value decomposition (SVD)\n",
    "\n",
    "- Representing matrix in the form\n",
    "\n",
    "$$ A = U\\Sigma V^*,$$\n",
    "\n",
    "- $U$ and $V$ are unitary\n",
    "- $\\Sigma$ is diagonal\n",
    "\n",
    "- It exists for **any** matrix\n",
    "- It is used to get the best low rank approximation of the given matrix in Frobenius and spectral norms (Eckart–Young theorem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Problem 1: solving linear systems\n",
    "\n",
    "$$ Ax = b $$\n",
    "\n",
    "- LU/Cholesky decomposition\n",
    "- Least-squares problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Case 1: square matrix\n",
    "\n",
    "A solution to the linear system of equations with a square matrix $A$\n",
    "\n",
    "$$ Ax = b $$\n",
    "\n",
    "exists, iff\n",
    "\n",
    "- $\\det A \\neq 0$\n",
    "\n",
    "or\n",
    "\n",
    "- matrix $A$ has full rank."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How to solve this problem?\n",
    "\n",
    "- Gaussian elimination or LU decomposition\n",
    "\n",
    "**Definition**: LU decomposition of the matrix $A$ is representation in the form\n",
    "\n",
    "$$A =  LU,$$\n",
    "\n",
    "where $L$ – **lower triangular** and $U$ – **upper triangular** matrix.\n",
    "\n",
    "- Forward and backward steps in Gaussian elimination\n",
    "    - Forward step\n",
    "\n",
    "    $$ L y = b.$$\n",
    "\n",
    "    - Backward step\n",
    "    \n",
    "    $$ U x = y.$$\n",
    "\n",
    "Does LU decomposition always exist?\n",
    "- Complexity is $\\mathcal{O}(n^3)$\n",
    "- **Q:** can we reduce it?\n",
    "- Problems with stability\n",
    "- Strictly regular matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Pivoting\n",
    "\n",
    "- Any matrix can be factorized with PLU decomposition\n",
    "\n",
    "$$ A = PLU, $$\n",
    "\n",
    "where $P$ is a permutation matrix.\n",
    "\n",
    "- What about complexity?\n",
    "- How to solve linear system with this decomposition?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### If matrix is Hermitian and positive definite\n",
    "\n",
    "- Hermitian positive definite matrix $A$ is strictly regular (check Sylvester's criterion) and has Cholesky decomposition\n",
    "\n",
    "$$A = RR^*,$$\n",
    "\n",
    "where $R$ is lower triangular matrix.\n",
    "\n",
    "- Sometimes matrix $R$ is called \"square root\" of the matrix $A$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Stability of solving linear system\n",
    "\n",
    "- Algorithms in exact and floating-point arithmetics behave **differently**!\n",
    "- Concept of **condition number** is introduced to estimate this difference in solving linear systems\n",
    "\n",
    "$$\\mathrm{cond}(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert.$$\n",
    "\n",
    "- It leads to the following estimate\n",
    "\n",
    "$$ \\frac{\\Vert \\widehat{x} - x \\Vert}{\\Vert x \\Vert} \\leq \\mathrm{cond}(A) \\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big),$$\n",
    "\n",
    "where $(A + \\Delta A) \\widehat{x} = f + \\Delta f.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Case 2: rectangular matrix\n",
    "\n",
    "- the number of rows is **larger** than the number of columns – overdetermined linear system\n",
    "    - a solution may not exist\n",
    "- the number of rows is **smaller** than the number of columns – underdetermined linear system\n",
    "    - a solution is not unique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Least-squares problem \n",
    "\n",
    "$$\\Vert A x - b \\Vert^2_2 \\rightarrow \\min_x$$\n",
    "\n",
    "- Geometric interpretation\n",
    "- Standard problem: you need to call a proper function from the package and understand its output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gram matrix and normal equation\n",
    "\n",
    "- From the first order optimality condition we get normal equation\n",
    "\n",
    "$$\n",
    "A^* A x = A^* b\n",
    "$$\n",
    "\n",
    "Matrix $A^* A$ is called **Gram matrix**.\n",
    "\n",
    "Property of the normal equation:\n",
    "\n",
    "- Condition number of the matrix $A^* A$ is square of the condition number of the matrix $A$ (check this fact!).\n",
    "- Therefore solving normal equation in this form is not a good idea!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Pseudoinverse matrix\n",
    "\n",
    "- If matrix $A$ is rank-deficient, then Gram matrix is singular\n",
    "- Therefore, introduce concept of the pseudoinverse matrix $A^{\\dagger}$\n",
    "\n",
    "$$x = A^{\\dagger} b.$$\n",
    "\n",
    "- Matrix \n",
    "\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^*$$ \n",
    "\n",
    "is called pseudoinverse matrix for the matrix $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Pseudoinverse matrix property\n",
    "\n",
    "* If matrix $A$ is full-rank, then $A^* A$ is invertible, and we get\n",
    "\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^*.$$\n",
    "\n",
    "* If matrix $A$ is square and non-singular, then we get\n",
    "\n",
    "$$ A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^* = A^{-1} A^{-*} A^* = A^{-1}$$\n",
    "\n",
    "We have obtained standard inverion of the matrix $A$\n",
    "\n",
    "* If $A$ has linear dependent columns, then $A^\\dagger b$ gives solution with minimum Euclidean norm - remind about the geometry of the problem!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Computing pseudoinverse through SVD\n",
    "\n",
    "Let $A = U \\Sigma V^*$ be SVD of the matrix $A$. Then,\n",
    "\n",
    "$$A^{\\dagger} = V \\Sigma^{\\dagger} U^*,$$\n",
    "\n",
    "where $\\Sigma^{\\dagger}$ is diagonal matrix such that diagonal entries are inverse of non-zero singular values of the matrix $A$. Indeed,\n",
    "\n",
    "\\begin{align*}\n",
    "A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = \\lim_{\\alpha \\rightarrow 0}( \\alpha VV^* + V \\Sigma^2 V^*)^{-1} V \\Sigma U^* = \\lim_{\\alpha \\rightarrow 0}( V(\\alpha I + \\Sigma^2) V^*)^{-1} V \\Sigma U^* = V \\lim_{\\alpha \\rightarrow 0}(\\alpha I + \\Sigma^2)^{-1} \\Sigma U^* = V \\Sigma^{\\dagger} U^*,\n",
    "\\end{align*}\n",
    "\n",
    "* You can check that $\\Sigma^{\\dagger}$ is constructed with inversion of singular values\n",
    "* If some singular values are small, then you can skip them. This gives more stable solution with respect to the noise in the right-hand side.\n",
    "\n",
    "**Q:** what about condition number?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### To solve the problem we use QR decomposition\n",
    "\n",
    "Any matrix can be represented in the form\n",
    "\n",
    "$$\n",
    "A = Q R,\n",
    "$$\n",
    "\n",
    "where $Q$ is unitary matrix, and $R$ is upper tringular matrix.\n",
    "\n",
    "- Gram-Schmidt orthogonalization is unstable (cf. Problem 2 in PS2) \n",
    "- Modified version is more stable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Computing with Householder reflections. This transformation zeroes subvector\n",
    "\n",
    "$$\n",
    "    H \\begin{bmatrix} \\times \\\\ \\times \\\\ \\times \\\\ \\times  \\end{bmatrix} = \n",
    "      \\begin{bmatrix} \\times \\\\ 0 \\\\ 0 \\\\ 0  \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "- Householder matrix \n",
    "\n",
    "$$ H = I - 2 vv^*. $$\n",
    "\n",
    "Complexity of Householder matrix by vector product is $\\mathcal{O}(n)$!\n",
    "\n",
    "- Complexity of the computing QR decomposition is $\\mathcal{O}(n^3)$\n",
    "- What is complexity for the matrix $m \\times n,$ where $m > n$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Computing with Givens matrix (rotations)\n",
    "- Givens matrix in 2D has the form\n",
    "\n",
    "$$\n",
    "    G = \\begin{bmatrix}\n",
    "          \\cos \\alpha & -\\sin \\alpha \\\\\n",
    "          \\sin \\alpha & \\cos \\alpha\n",
    "        \\end{bmatrix},\n",
    "$$\n",
    "\n",
    "- It zeroes elements sequentially\n",
    "- It is appropriate for sparse matrices and parallel computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Come back to the least-squares problem\n",
    "\n",
    "- If matrix $A$ is full-rank, then\n",
    "\n",
    "$$\n",
    "x = A^{\\dagger}b = (A^*A)^{-1}A^*b = ((QR)^*(QR))^{-1}(QR)^*b = (R^*Q^*QR)^{-1}R^*Q^*b = R^{-1}Q^*b. \n",
    "$$ \n",
    "\n",
    "- Thus, we have to solve square linear system\n",
    "\n",
    "$$\n",
    "Rx = Q^* b.\n",
    "$$\n",
    "\n",
    "- $R$ is upper-triangular\n",
    "- The complexity of solving is $\\mathcal{O}(n^2)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Problem 2: eigenvalue problem\n",
    "\n",
    "$$ Ax = \\lambda x, \\qquad A = S \\Lambda S^{-1} $$\n",
    "\n",
    "- Power method\n",
    "- QR algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Theoretical approach\n",
    "\n",
    "- Characteristic equation\n",
    "\n",
    "$$ \\det(A - \\lambda I) = 0 $$\n",
    "\n",
    "- Solve it with polynomial root finding procedure\n",
    "- It is very unstable operation! \n",
    "- Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Practical approach: partial eigenvalue problem\n",
    "\n",
    "- Power method\n",
    "- One iteration reads as\n",
    "\n",
    "$$ x_{k+1} = A x_k, \\quad x_{k+1} := \\frac{x_{k+1}}{\\Vert x_{k+1} \\Vert_2}.$$\n",
    "\n",
    "- It converges (not always!) to the maximum by modulus eigenvalue\n",
    "- To perform iterations you need only procedure for computing matrix by vector product (cf. Problem 4 in PS2)\n",
    "\n",
    "**Q:** what about complexity?\n",
    "\n",
    "- Convergence is linear with factor $q = \\left|\\frac{\\lambda_{2}}{\\lambda_{1}}\\right| < 1$, where $\\lambda_1>\\lambda_2\\geq\\dots\\geq \\lambda_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Krylov subspace\n",
    "\n",
    "- Eigenvector computed with power method lies in **Krylov subspace** $\\{x_0, Ax_0,\\dots,A^{k}x_0\\}$ and has the form $\\mu A^k x_0$, where $\\mu$ is normalized constant. \n",
    "- Further Krylov subspace will be a key tool to solve large-scale problems!\n",
    "- More details about why Krylov subpace is so important will be later in the course"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Practical approach: full eigenvalue problem\n",
    "\n",
    "- We want to find all eigenvalues!\n",
    "- Schur decomposition\n",
    "\n",
    "$$A = UTU^*,$$ \n",
    "\n",
    "where $U$ is unitary, and $T$ is upper triangular.\n",
    "\n",
    "- Spectra of $T$ and $A$ are the same since the matrix $U$ is unitary\n",
    "- How to compute Schur form?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  QR algorithm\n",
    "1. Initialize $A_0 = A$.\n",
    "2. Compute QR decomposition of the matrix $A_k$: $A_k = Q_k R_k$.\n",
    "3. Update approximation $A_{k+1} = R_k Q_k$.\n",
    "\n",
    "\n",
    "- How to check convergence?\n",
    "- What about complexity?\n",
    "- How to find matrices $U$ and $T$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Converegence and acceleration\n",
    "\n",
    "- Accelerate QR algorithm up to $\\mathcal{O}(n^3)$\n",
    "- Use property of the matrix in upper-hessenberg form\n",
    "- Matrix $A$ is in upper-hessenberg form, if\n",
    "\n",
    "$$a_{ij} = 0, \\quad \\mbox{if } i \\geq j+2.$$\n",
    "\n",
    "or \n",
    "\n",
    "$$A = \\begin{bmatrix} * & * & * & * & * \\\\ * & * & * & * & * \\\\ 0 & * & * & * & *\\\\ 0 & 0 & * & * & *\\\\ 0 & 0 & 0 & * & * \\\\ \\end{bmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Where are we moving further?\n",
    "\n",
    "- Problems are the same, but the scale is much larger!\n",
    "- Exploiting the structure of matrices\n",
    "    - sparsity\n",
    "    - low rank\n",
    "    - Toeplitz and circulant matrices\n",
    "- Iterative methods for\n",
    "    - solving linear systems\n",
    "    - computing eigenvalues and eigenvectors \n",
    "    - computing **matrix functions**\n",
    "- Tensor decompoisitions and their applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important dates \n",
    "\n",
    "- November, 21 –– midterm\n",
    "- November, 24 –– deadline for PS2\n",
    "- November, 24 –– deadline for submiting project proposals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plans for next classes\n",
    "\n",
    "- Tomorrow Q&A session\n",
    "- Midterm is on Thursday\n",
    "- Next regular lecture is on Monday"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}