lectures/lecture7/lecture-7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 7: Matrix decompositions and how we compute them"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "- Eigenvectors and eigenvalues\n",
    "- Characteristic polynomial and why it is a bad idea\n",
    "- Power method: $x := Ax$\n",
    "- Gershgorin theorem\n",
    "- Schur theorem: $A = U T U^*$ \n",
    "- Normal matrices: $A^* A = A A^*$\n",
    "- Advanced topic: pseudospectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture\n",
    "\n",
    "- Today we will talk about **matrix factorizations** as general tool\n",
    "\n",
    "- Basic **matrix factorizations** in numerical linear algebra:\n",
    "\n",
    "    - LU decomposition and Gaussian elimination &mdash; already covered\n",
    "    - QR decomposition and Gram-Schmidt algorithm \n",
    "    - Schur decomposition and QR-algorithm\n",
    "    - Methods for computing SVD decomposition in the second part\n",
    "\n",
    "- We already introduced QR decomposition some time ago, but now we are going to discuss it in more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## General concept of matrix factorization\n",
    "\n",
    "- In numerical linear algebra we need to solve different tasks, for example:\n",
    "\n",
    "    - Solve linear systems $Ax = f$\n",
    "    - Compute eigenvalues / eigenvectors\n",
    "    - Compute singular values / singular vectors\n",
    "    - Compute inverses, even sometimes determinants \n",
    "    - Compute **matrix functions** like $\\exp(A), \\cos(A)$ (these are not elementwise functions)\n",
    "\n",
    "- In order to do this, we represent the matrix as a sum and/or product of matrices with **simpler structure**,   such that we can solve mentioned tasks faster / in a more stable form.\n",
    "\n",
    "- What is a **simpler structure**?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What is a simpler structure\n",
    "\n",
    "- We already encountered several classes of matrices **with structure**. \n",
    "\n",
    "- For dense matrices the most important classes are\n",
    "\n",
    "    - **unitary matrices**\n",
    "    - **upper/lower triangular matrices** \n",
    "    - **diagonal matrices**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other classes of structured matrices\n",
    "\n",
    "- For **sparse matrices** the sparse constraints are often included in the factorizations.\n",
    "- For **Toeplitz matrices** an important class of matrices is the class of matrices with **low displacement rank**, which is based on the low-rank matrices. \n",
    "- The class of **low-rank matrices** and **block low-rank matrices** appears in many applications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plan\n",
    "The plan for today's lecture is to discuss the decompositions one-by-one and point out:\n",
    "- How to compute a particular decomposition\n",
    "- When the decomposition exists \n",
    "- What is done in **real life** (LAPACK)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Decompositions we want to discuss today\n",
    "- LU factorization & Cholesky factorization &mdash; quick reminder, already done.\n",
    "- QR decomposition and Gram-Schmidt algorithm\n",
    "- One slide about the SVD (more details in the second part of today lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## PLU decomposition\n",
    "\n",
    "- Any nonsingular matrix can be factored as \n",
    "\n",
    "$$ A = P L U, $$\n",
    "\n",
    "where $P$ is a permutation matrix, $L$ is a **lower triangular matrix**, $U$ is an **upper triangular**\n",
    "\n",
    "- **Main goal** of the LU decomposition is to solve linear systems, because\n",
    "\n",
    "$$ A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, $$\n",
    "\n",
    "and this reduces to the solution of two linear systems \n",
    "\n",
    "$$ L y = f,  \\quad U x = y $$\n",
    "\n",
    "with lower and upper triangular matrices respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Positive definite matrices and Cholesky decomposition, reminder\n",
    "\n",
    "If the matrix is Hermitian positive definite, i.e. \n",
    "\n",
    "$$ A = A^*, \\quad (Ax, x) > 0, \\quad x \\ne 0, $$\n",
    "\n",
    "then it can be factored as \n",
    "\n",
    "$$ A = RR^*, $$\n",
    "\n",
    "where $R$ is lower triangular.\n",
    "\n",
    "We will need this for the QR decomposition.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition\n",
    "\n",
    "- The next decomposition: **QR** decomposition. \n",
    "- Again from the name it is clear that a matrix is represented as a product  \n",
    "\n",
    "$$\n",
    "    A = Q R, \n",
    "$$\n",
    "\n",
    "where $Q$ is an **column orthogonal (unitary)** matrix and $R$ is **upper triangular**.  \n",
    "\n",
    "- The matrix sizes: $Q$ is $n \\times m$, $R$ is $m \\times m$ if $n\\geq m$. See our [poster](../../decompositions.pdf) for visualization of QR decomposition\n",
    "\n",
    "- QR decomposition is defined for **any rectangular matrix**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition: applications\n",
    "\n",
    "This decomposition plays a crucial role in many problems:\n",
    "- Computing orthogonal bases in a linear space\n",
    "- Used in the preprocessing step for the SVD\n",
    "- QR-algorithm for the computation of eigenvectors and eigenvalues ([one of the 10 most important algorithms of the 20th century](https://archive.siam.org/pdf/news/637.pdf)) is based on the QR decomposition\n",
    "- Solving overdetermined systems of linear equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition and least squares, reminder\n",
    "\n",
    "- Suppose we need to solve\n",
    "\n",
    "$$ \\Vert A x - f \\Vert_2 \\rightarrow \\min_x, $$\n",
    "\n",
    "where $A$ is $n \\times m$, $n \\geq m$.\n",
    "\n",
    "- Then we factorize\n",
    "\n",
    "$$ A = Q R, $$\n",
    "\n",
    "and use equation for pseudo-inverse matrix in the case of the full rank matrix $A$:\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",
    "thus $x$ can be recovered from \n",
    "\n",
    "$$R x = Q^*b$$\n",
    "\n",
    "- Note that this is a square system of linear equations with lower triangular matrix. What is the complexity of solving this system?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Short side note: randomized least squares\n",
    "\n",
    "- One of the efficient ways to solve really overdetermined ($n\\gg m$) system of linear equations is to use **Kaczmarz method**.\n",
    "\n",
    "- Instead of solving all equations, pick one randomly, which reads\n",
    "\n",
    "$$a^{\\top}_i x = f_i,$$\n",
    "\n",
    "and given an approximation $x_k$ try to find $x_{k+1}$ as \n",
    "\n",
    "$$x_{k+1} = \\arg \\min_x \\Vert x - x_k \\Vert, \\quad \\mbox{s.t.} \\quad  a^{\\top}_i x = f_i.$$\n",
    "\n",
    "- A simple analysis gives \n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{(a_i, x_k) - f_i}{(a_i, a_i)} a_i. $$\n",
    "\n",
    "- A cheap update, but the analysis is quite complicated.\n",
    "- You can recognize in this method stochastic gradient descent with specific step size equal to $\\frac{1}{\\|a_i\\|_2^2}$ for every sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Existence of QR decomposition\n",
    "**Theorem.**\n",
    "Every rectangular $n \\times m$  matrix has a QR decomposition. \n",
    "\n",
    "\n",
    "There are several ways to prove it and compute it:\n",
    "\n",
    "- Theoretical: using the Gram matrices and Cholesky factorization\n",
    "- Geometrical: using the Gram-Schmidt orthogonalization\n",
    "- Practical: using Householder/Givens transformations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Proof using Cholesky decomposition\n",
    "\n",
    "If we have the representation of the form\n",
    "\n",
    "$$A = QR,$$\n",
    "\n",
    "then $A^* A = ( Q R)^* (QR)  = R^* (Q^* Q) R = R^* R$, the matrix $A^* A$ is called **Gram matrix**, and its elements are scalar products of the columns of $A$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Proof using Cholesky decomposition (full column rank)\n",
    "\n",
    "- Assume that $A$ has **full column rank**. Then, it is simple to show that $A^* A$ is positive definite:\n",
    "\n",
    "$$ (A^* A y, y) = (Ay, Ay) = \\Vert Ay \\Vert^2  > 0, \\quad y\\not = 0. $$\n",
    "\n",
    "- Therefore, $A^* A = R^* R$ always exists.\n",
    "\n",
    "- Then the matrix $Q = A R^{-1}$ is unitary:  \n",
    "\n",
    "$$ (A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Proof using Cholesky decomposition (rank-deficient case)\n",
    "\n",
    "- When an $n \\times m$ matrix does not have full column rank, it is said to be **rank-deficient**.  \n",
    "\n",
    "- The QR decomposition, however, also exists.\n",
    "\n",
    "- For any **rank-deficient matrix** there is a sequence of full-column rank matrices $A_k$ such that $A_k \\rightarrow A$ (why?).\n",
    "\n",
    "- Each $A_k$ can be decomposed as $A_k = Q_k R_k$. \n",
    "\n",
    "- The set of all unitary matrices is compact, thus there exists a converging subsequence $Q_{n_k} \\rightarrow Q$ (why?), and $Q^* A_k \\rightarrow Q^* A = R$, which is triangular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Stability of QR decomposition via Cholesky decomposition\n",
    "\n",
    "- So, the simplest way to compute QR decomposition is then\n",
    "\n",
    "$$A^* A = R^* R,$$\n",
    "\n",
    "and\n",
    "\n",
    "$$Q = A R^{-1}.$$\n",
    "\n",
    "- It is a **bad idea** for numerical stability. Let us do some demo (for a submatrix of the Hilbert matrix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Built-in QR orth 1.595746790565955e-15\n",
      "Via Cholesky: 0.8819197747599913\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DeviceArray(-2.1925951e-11, dtype=float64)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "from jax.config import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "n = 40\n",
    "r = 9\n",
    "a = [[1.0 / (i + j + 0.5) for i in range(r)] for j in range(n)]\n",
    "a = jnp.array(a)\n",
    "q, Rmat = jnp.linalg.qr(a)\n",
    "e = jnp.eye(r)\n",
    "print('Built-in QR orth', jnp.linalg.norm(jnp.dot(q.T, q) - e))\n",
    "gram_matrix = a.T.dot(a)\n",
    "Rmat1 = jnp.linalg.cholesky(gram_matrix)\n",
    "q1 = jnp.dot(a, jnp.linalg.inv(Rmat1.T))\n",
    "#q1 = jnp.linalg.solve(Rmat1, a.T).T\n",
    "print('Via Cholesky:', jnp.linalg.norm(jnp.dot(q1.T, q1) - e))\n",
    "\n",
    "\n",
    "q1[:, 0]@q1[:, 5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Second way: Gram-Schmidt orthogonalization\n",
    "\n",
    "- QR decomposition is a mathematical way of writing down the Gram-Schmidt orthogonalization process.  \n",
    "- Given a sequence of vectors $a_1, \\ldots, a_m$ we want to find orthogonal basis $q_1, \\ldots, q_m$ such that every $a_i$ is a linear combination of such vectors.  \n",
    "\n",
    "**Gram-Schmidt:**\n",
    "1. $q_1 := a_1/\\Vert a_1 \\Vert$\n",
    "2. $q_2 := a_2 - (a_2, q_1) q_1, \\quad q_2 := q_2/\\Vert q_2 \\Vert$\n",
    "3. $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \\quad q_3 := q_3/\\Vert q_3 \\Vert$\n",
    "4. And so on  \n",
    "\n",
    "Note that the transformation from $Q$ to $A$ has triangular structure, since from the $k$-th vector we subtract only the previous ones. It follows from the fact that the product of triangular matrices is a triangular matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Modified Gram-Schmidt\n",
    "- Gram-Schmidt can be very unstable (i.e., the produced vectors will be not orthogonal, especially if $q_k$ has small norm).  \n",
    "This is called **loss of orthogonality**.  \n",
    "\n",
    "- There is a remedy, called **modified Gram-Schmidt** method. Instead of doing \n",
    "\n",
    "$$q_k := a_k - (a_k, q_1) q_1 - \\ldots - (a_k, q_{k-1}) q_{k-1}$$\n",
    "\n",
    "we do it step-by-step. First we set $q_k := a_k$ and orthogonalize sequentially:\n",
    "\n",
    "$$\n",
    "   q_k := q_k - (q_k, q_1)q_1, \\quad q_k := q_{k} - (q_k,q_2)q_2, \\ldots\n",
    "$$\n",
    "\n",
    "- In exact arithmetic, it is the same. In floating point it is absolutely different!\n",
    "\n",
    "- Note that the complexity is $\\mathcal{O}(nm^2)$ operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition: the (almost) practical way\n",
    "\n",
    "- If $A = QR$, then  \n",
    "\n",
    "$$ R = Q^* A, $$\n",
    "\n",
    "and we need to find a certain orthogonal matrix $Q$ that brings a matrix into upper triangular form.  \n",
    "- For simplicity, we will look for an $n \\times n$ matrix such that\n",
    "\n",
    "$$ Q^* A = \\begin{bmatrix} * & * & *  \\\\ 0 & * & * \\\\ 0 & 0 & * \\\\ & 0_{(n-m) \\times m} \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- We will do it column-by-column.  \n",
    "- First, we find a Householder matrix $H_1 = (I - 2 uu^{\\top})$ such that (we illustrate on a $4 \\times 3$ matrix)\n",
    "\n",
    "$$ H_1 A = \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & * & * \\\\ 0 & * & * \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then, \n",
    "\n",
    "$$ H_2 H_1 A = \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & 0 & * \\\\ 0 & 0 & * \\end{bmatrix}, $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ H_2 = \\begin{bmatrix} 1 & 0 \\\\ 0 & H'_2, \\end{bmatrix} $$\n",
    "\n",
    "and $H'_2$ is a $3 \\times 3$ Householder matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Finally, \n",
    "\n",
    "$$ H_3 H_2 H_1 A = \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & 0 & * \\\\ 0 & 0 & 0 \\end{bmatrix}, $$\n",
    "\n",
    "You can try to implement it yourself, it is simple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition: real life\n",
    "\n",
    "- In reality, since this is a dense matrix factorization, you should implement the algorithm in terms of blocks (why?).  \n",
    "\n",
    "- Instead of using Householder transformation, we use **block Householder** transformation of the form\n",
    "\n",
    "$$H = I - 2UU^*, $$\n",
    "\n",
    "where $U^* U = I$.\n",
    "\n",
    "- This allows us to use BLAS-3 operations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Rank-revealing QR-decomposition\n",
    "\n",
    "- The QR-decomposition can be also used to compute the (numerical) rank of the matrix, see \n",
    "[Rank-Revealing QR Factorizations and the Singular Value Decomposition, Y. P. Hong; C.-T. Pan](http://scgroup.hpclab.ceid.upatras.gr/class/SCII/Various/HongPan.pdf)\n",
    "\n",
    "- It is done via so-called **rank-revealing** factorization.  \n",
    "\n",
    "- It is based on the representation\n",
    "\n",
    "$$P A = Q R,$$\n",
    "\n",
    "where $P$ is the permutation matrix (it permutes columns), and  $R$ has the block form\n",
    "\n",
    "\n",
    "$$R = \\begin{bmatrix} R_{11} & R_{12} \\\\ 0 & R_{22}\\end{bmatrix}.$$\n",
    "\n",
    "- The goal is to find $P$ such that the norm of $R_{22}$ is **small**, so you can find the numerical rank by looking at it.\n",
    "\n",
    "- An estimate is $\\sigma_{r+1} \\leq \\Vert R_{22} \\Vert_2$ (check why)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "- LU and QR decompositions can be computed using **direct methods** in finite amount of operations.\n",
    "\n",
    "- What about Schur form and SVD?  \n",
    "\n",
    "- They can not be computed by **direct methods** (why?) they can only be computed by **iterative methods**.\n",
    "\n",
    "- Although iterative methods still have the same $\\mathcal{O}(n^3)$ complexity in floating point arithmetic thanks to fast convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Schur form\n",
    "\n",
    "- Recall that every matrix can be written in the **Schur form**\n",
    "\n",
    "$$A = Q T Q^*,$$\n",
    "\n",
    "with upper triangular $T$ and unitary $Q$\n",
    "and this decomposition gives eigenvalues of the matrix (they are on the diagonal of $T$).\n",
    "\n",
    "- The first and the main algorithm for computing the Schur form is the **QR algorithm**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR algorithm\n",
    "\n",
    "- The QR algorithm was independently proposed in 1961 by Kublanovskaya and Francis.   \n",
    "\n",
    "- <font color='red'> Do not **mix** QR algorithm and QR decomposition! </font>\n",
    "\n",
    "- QR decomposition is the representation of a matrix, whereas QR algorithm uses QR decomposition to compute the eigenvalues!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Way to QR algorithm\n",
    "\n",
    "- Consider the equation\n",
    "\n",
    "$$A = Q T Q^*,$$\n",
    "\n",
    "and rewrite it in the form\n",
    "\n",
    "$$ Q T = A Q. $$\n",
    "\n",
    "- On the left we can see QR factorization of the matrix $AQ$.\n",
    "\n",
    "- We can use this to derive fixed-point iteration for the Schur form, also known as **QR algorithm.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Derivation of the QR algorithm as fixed-point iteration\n",
    "We can write down the iterative process \n",
    "\n",
    "$$\n",
    "    Q_{k+1} R_{k+1} = A Q_k, \\quad Q_{k+1}^* A = R_{k+1} Q^*_k\n",
    "$$\n",
    "\n",
    "Introduce\n",
    "\n",
    "$$A_k = Q^* _k A Q_k = Q^*_k Q_{k+1} R_{k+1} = \\widehat{Q}_k R_{k+1}$$\n",
    "\n",
    "and the new approximation reads\n",
    "\n",
    "$$A_{k+1} = Q^*_{k+1} A Q_{k+1} = ( Q_{k+1}^* A = R_{k+1} Q^*_k)  = R_{k+1} \\widehat{Q}_k.$$\n",
    "\n",
    "So we arrive at the standard form of the QR algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The final formulas are then written in the classical **QRRQ**-form:\n",
    "\n",
    "1. Start from $A_0 = A$.\n",
    "2. Compute QR factorization of $A_k = Q_k R_k$.\n",
    "3. Set $A_{k+1} = R_k Q_k$.\n",
    "\n",
    "Iterate until $A_k$ is **triangular enough** (e.g. norm of subdiagonal part is small enough).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What about the convergence and complexity\n",
    "\n",
    "**Statement**\n",
    "\n",
    "Matrices $A_k$ are unitary similar to $A$\n",
    "\n",
    "$$A_k = Q^*_{k-1} A_{k-1} Q_{k-1} = (Q_{k-1} \\ldots Q_1)^* A (Q_{k-1} \\ldots Q_1)$$\n",
    "\n",
    "and the product of unitary matrices is a unitary matrix.\n",
    "\n",
    "- Complexity of each step is $\\mathcal{O}(n^3)$, if a general QR decomposition is done.\n",
    "\n",
    "- Our hope is that $A_k$ will be **very close to the triangular  matrix** for suffiently large $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Leading 3x3 block of a:\n",
      "[[ 2.45414517e+00  1.95145959e-08  6.60861102e-18  5.76974836e-17]\n",
      " [ 1.95145961e-08  4.14159204e-01  2.66526168e-13  7.62304085e-18]\n",
      " [ 9.70054052e-21  2.66524048e-13  2.46972338e-02 -3.25466594e-18]\n",
      " [ 3.24776145e-36  9.72713876e-29  1.32450609e-17  7.05827401e-04]]\n"
     ]
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "n = 5\n",
    "a = [[1.0/(i + j + 0.5) for i in range(n)] for j in range(n)]\n",
    "niters = 10\n",
    "for k in range(niters):\n",
    "    q, rmat = jnp.linalg.qr(a)\n",
    "    a = rmat.dot(q)\n",
    "print('Leading 3x3 block of a:')\n",
    "print(a[:4, :4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence and complexity of the QR algorithm\n",
    "\n",
    "- The convergence of the QR algorithm is from the **largest eigenvalues** to the smallest.\n",
    "\n",
    "- At least 2-3 iterations is needed for an eigenvalue. \n",
    "\n",
    "- Each step is one QR factorization and one matrix-by-matrix product, as a result $\\mathcal{O}(n^3)$ complexity.\n",
    "\n",
    "**Q**: does it mean $\\mathcal{O}(n^4)$ complexity totally? \n",
    "\n",
    "**A**: fortunately, not. \n",
    "\n",
    "- We can speedup the QR algorithm by using **shifts**, since $A_k - \\lambda I$ has the same Schur vectors.\n",
    "\n",
    "- We will discuss these details later"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of QR algorithm\n",
    "\n",
    "Does QR algorithm always convergence? \n",
    "\n",
    "Provide an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Counterexample\n",
    "\n",
    "For a matrix $A = \\begin{bmatrix} 0 & 1 \\\\\n",
    "                                  1 & 0 \\end{bmatrix}$\n",
    "\n",
    "we have $A_k = A$. \n",
    "          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Connection to orthogonal iteration\n",
    "\n",
    "In the previous lecture, we considered **power iteration**, which is $A^k v$ -- approximation of the eigenvector. \n",
    "\n",
    "The QR algorithm computes (implicitly) QR-factorization of the matrix $A^k$:\n",
    "\n",
    "$$A^k = A \\cdot \\ldots \\cdot A = Q_1 R_1 Q_1 R_1 \\ldots = Q_1 Q_2 R_2 Q_2 R_2 \\ldots (R_2 R_1) = \\ldots (Q_1 Q_2 \\ldots Q_k) (R_k \\ldots R_1).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A few words about the SVD\n",
    "\n",
    "- Last but not least: the **singular value decomposition** of matrix.\n",
    "\n",
    "$$A = U \\Sigma V^*.$$\n",
    "\n",
    "- We can compute it via eigendecomposition of \n",
    "\n",
    "$$A^* A = V^* \\Sigma^2 V,$$\n",
    "\n",
    "and/or\n",
    "\n",
    "$$AA^* = U^* \\Sigma^2 U$$\n",
    "\n",
    "with QR algorithm, but it is a **bad idea** (c.f. Gram matrices).\n",
    "\n",
    "- We will discuss methods for computing SVD later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bonus: CS (Cosine-Sine) decomposition\n",
    "\n",
    "- Let $Q$ be square unitary matrix with even number of rows and it is splitted in four equal size blocks\n",
    "\n",
    "$$ Q = \\begin{bmatrix} Q_{11} & Q_{12} \\\\ Q_{21} & Q_{22} \\end{bmatrix} $$ \n",
    "\n",
    "- Then there exist unitary matrices $U_1, U_2, V_1, V_2$ such that\n",
    "\n",
    "$$ \\begin{bmatrix} U_1^* & 0 \\\\ 0 & U_2^* \\end{bmatrix} \\begin{bmatrix} Q_{11} & Q_{12} \\\\ Q_{21} & Q_{22} \\end{bmatrix} \\begin{bmatrix} V_1 & 0 \\\\ 0 & V_2 \\end{bmatrix} = \\begin{bmatrix} C & S \\\\ -S & C \\end{bmatrix}, $$\n",
    "\n",
    "where $C = \\mathrm{diag}(c)$ and $S = \\mathrm{diag}(s)$ such that $c_i \\geq 0, s_i \\geq 0$ and $c_i^2 + s_i^2 = 1$ \n",
    "\n",
    "\n",
    "- **Q**: how many SVD  do we have inside the CS decomposition?\n",
    "\n",
    "- The case of rectangular matrix with orthonormal columns\n",
    "\n",
    "$$ \\begin{bmatrix} U_1^* & 0 \\\\ 0 & U_2^* \\end{bmatrix} \\begin{bmatrix} Q_{1} \\\\ Q_{2} \\end{bmatrix} V = \\begin{bmatrix} C \\\\ S \\end{bmatrix} $$\n",
    "\n",
    "- The algorithm for computing this decomposition is presented [here](https://arxiv.org/pdf/0707.1838.pdf)\n",
    "\n",
    "- This decomposition naturally arises in the problem of finding distances and angles between subspaces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "- QR decomposition and Gram-Schmidt algorithm, reduction to a simpler form by Householder transformations\n",
    "- Schur decomposition and QR algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next steps\n",
    "\n",
    "- Efficient implementation of QR algorithm, convergence properties\n",
    "- Efficient computation of the SVD: 4 algorithms\n",
    "- More applications of the SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Слайд-шоу",
  "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.8.5"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}