lectures/lecture-2/lecture-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": [
     1
    ],
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 2. Matrix norms and unitary matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "\n",
    "- Floating point  (double, single, half precisions, number of bytes), rounding error\n",
    "- Vector norms are measures of smallness, used to compute the distance and accuracy\n",
    "- Forward/backward error (and stability of algorithms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Notations\n",
    "\n",
    "We use notation \n",
    "\n",
    "$$A= \\begin{bmatrix} a_{11} & \\dots & a_{1m} \\\\ \\vdots & \\ddots & \\vdots \\\\ a_{n1} & \\dots & a_{nm}  \\end{bmatrix} \\equiv \\{a_{ij}\\}_{i,j=1}^{n,m}\\in \\mathbb{C}^{n\\times m}.$$\n",
    "\n",
    "$A^*\\stackrel{\\mathrm{def}}{=}\\overline{A^\\top}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrices and norms\n",
    "\n",
    "- Recall vector norms that allow to evaluate distance between two vectors or how large are elements of a vector.\n",
    "\n",
    "- How to generalize this concept to matrices?\n",
    "\n",
    "- A trivial answer is that there is no big difference between matrices and vectors, and here comes the simplest matrix norm –– **Frobenius** norm:\n",
    "\n",
    "$$ \\Vert A \\Vert_F \\stackrel{\\mathrm{def}}{=} \\Big(\\sum_{i=1}^n \\sum_{j=1}^m |a_{ij}|^2\\Big)^{1/2} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix norms\n",
    "$\\Vert \\cdot \\Vert$ is called a **matrix norm** if it is a vector norm on the vector space of $n \\times m$ matrices:\n",
    "1. $\\|A\\| \\geq 0$ and if $\\|A\\| = 0$ then $A = O$\n",
    "3. $\\|\\alpha A\\| = |\\alpha| \\|A\\|$\n",
    "4. $\\|A+B\\| \\leq \\|A\\| + \\|B\\|$ (triangle inequality)\n",
    "\n",
    "Additionally some norms can satisfy the *submultiplicative property*\n",
    "\n",
    "* <span style=\"color:red\"> $\\Vert A B \\Vert \\leq \\Vert A \\Vert \\Vert B \\Vert$ </span>\n",
    "\n",
    "- These norms are called **submultiplicative norms**.\n",
    "\n",
    "- The submultiplicative property is needed in many places, for example in the estimates for the error of solution of linear systems (we will cover this topic later). \n",
    "\n",
    "- Example of a non-submultiplicative norm is Chebyshev norm \n",
    "\n",
    "$$ \\|A\\|_C = \\displaystyle{\\max_{i,j}}\\, |a_{ij}| $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Operator norms\n",
    "\n",
    "- The most important class of the matrix norms is the class of **operator norms**. They are defined as\n",
    "\n",
    "$$ \\Vert A \\Vert_{*,**} = \\sup_{x \\ne 0} \\frac{\\Vert A x \\Vert_*}{\\Vert x \\Vert_{**}}, $$\n",
    "\n",
    "where $\\Vert \\cdot \\Vert_*$ and $\\| \\cdot \\|_{**}$ are **vector norms**.\n",
    "\n",
    "- It is easy to check that operator norms are submultiplicative if $\\|\\cdot\\|_* = \\|\\cdot\\|_{**}$. Otherwise, it can be non-submultiplicative, think about example.\n",
    "\n",
    "- **Frobenius norm** is a matrix norm, but not an operator norm, i.e. you can not find $\\Vert \\cdot \\Vert_*$ and $\\| \\cdot \\|_{**}$ that induce it. \n",
    "- This is a nontrivial fact and the general criterion for matrix norm to be an operator norm can be found in [Theorem 6 and Corollary 4](http://www.sciencedirect.com/science/article/pii/S0024379515004346).\n",
    "For $\\Vert \\cdot \\Vert_* = \\| \\cdot \\|_{**}$ let us check on the blackboard!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix $p$-norms\n",
    "\n",
    "Important case of operator norms are matrix $p$-norms, which are defined for $\\|\\cdot\\|_* = \\|\\cdot\\|_{**} = \\|\\cdot\\|_p$. <br>\n",
    "\n",
    "Among all $p$-norms three norms are the most common ones:  \n",
    "\n",
    "- $p = 1, \\quad \\Vert A \\Vert_{1} = \\displaystyle{\\max_j \\sum_{i=1}^n} |a_{ij}|$.\n",
    "\n",
    "- $p = 2, \\quad$ spectral norm, denoted by $\\Vert A \\Vert_2$.\n",
    "\n",
    "- $p = \\infty, \\quad \\Vert A \\Vert_{\\infty} = \\displaystyle{\\max_i \\sum_{j=1}^m} |a_{ij}|$.\n",
    "\n",
    "Let us check it for $p=\\infty$ on a blackboard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Spectral norm\n",
    "\n",
    "- Spectral norm, $\\Vert A \\Vert_2$ is one of the most used matrix norms (along with the Frobenius norm). \n",
    "- It can not be computed directly from the entries using a simple formula, like the Frobenius norm, however, there are efficient algorithms to compute it.  \n",
    "- It is directly related to the **singular value decomposition** (SVD) of the matrix. It holds\n",
    "\n",
    "$$ \\Vert A \\Vert_2 = \\sigma_1(A) = \\sqrt{\\lambda_\\max(A^*A)} $$\n",
    "\n",
    "where $\\sigma_1(A)$ is the largest singular value of the matrix $A$ and $^*$ is a *conjugate transpose* of the matrix. \n",
    "\n",
    "- We will soon learn all about the SVD. Meanwhile, we can already compute the norm in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "import jax\n",
    "n = 100\n",
    "m = 2000\n",
    "a = jax.random.normal(jax.random.PRNGKey(0), (n, m)) #Random n x m matrix\n",
    "s1 = jnp.linalg.norm(a, 2) #Spectral\n",
    "s2 = jnp.linalg.norm(a, 'fro') #Frobenius\n",
    "s3 = jnp.linalg.norm(a, 1) #1-norm\n",
    "s4 = jnp.linalg.norm(a, jnp.inf) \n",
    "print('Spectral: {0:} \\nFrobenius: {1:} \\n1-norm: {2:} \\ninfinity: {3:}'.format(s1, s2, s3, s4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples\n",
    "\n",
    "Several examples of optimization problems where matrix norms arise:\n",
    "* <span style=\"color:red\"> $\\displaystyle{\\min_{\\mathrm{rank}(A_r)=r}}\\| A - A_r\\|$ </span> –– finding best rank-r approximation. SVD helps to solve this problem for $\\|\\cdot\\|_2$ and $\\|\\cdot\\|_F$.\n",
    "\n",
    "\n",
    "* $\\displaystyle{\\min_B}\\| P_\\Omega \\odot(A - B)\\| + \\mathrm{rank}(B)$ –– matrix completion. \n",
    "\n",
    "$$ (P_\\Omega)_{ij} = \\begin{cases} 1 & i,j\\in\\Omega \\\\ 0 & \\text{otherwise}, \\end{cases} $$\n",
    "\n",
    "where $\\odot$ denotes Hadamard product (elementwise)\n",
    "\n",
    "\n",
    "* $\\displaystyle{\\min_{B,C\\geq 0}} \\|A - BC\\|_F$ –– nonnegative matrix factorization. Symbol $B\\geq0$ here means that all elements of $B$ are nonnegative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Scalar product\n",
    "While norm is a measure of distance, the **scalar product** takes angle into account.  \n",
    "\n",
    "It is defined as\n",
    "\n",
    "* **For vectors:**\n",
    "$$ (x, y) =  x^* y = \\sum_{i=1}^n \\overline{x}_i y_i, $$\n",
    "where $\\overline{x}$ denotes the *complex conjugate* of $x$. The Euclidean norm is then\n",
    "\n",
    "$$ \\Vert x \\Vert_2 = \\sqrt{(x, x)}, $$\n",
    "\n",
    "or it is said the norm is **induced** by the scalar product.  \n",
    "\n",
    "\n",
    "* **For matrices** (Frobenius scalar product):\n",
    "\n",
    "$$ (A, B)_F = \\displaystyle{\\sum_{i=1}^{n}\\sum_{j=1}^{m}} \\overline{a}_{ij} b_{ij} \\equiv \\mathrm{trace}(A^* B), $$\n",
    "\n",
    "where $\\mathrm{trace}(A)$ denotes the sum of diagonal elements of $A$. One can check that $\\|A\\|_F = \\sqrt{(A, A)_F}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Remark**. The angle between two vectors is defined as\n",
    "\n",
    "$$ \\cos \\phi = \\frac{(x, y)}{\\Vert x \\Vert_2 \\Vert y \\Vert_2}. $$\n",
    "\n",
    "Similar expression holds for matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- An important property of the scalar product is the **Cauchy-Schwarz-Bunyakovski inequality**:\n",
    "\n",
    "$$|(x, y)| \\leq \\Vert x \\Vert_2 \\Vert y \\Vert_2,$$\n",
    "\n",
    "and thus the angle between two vectors is defined properly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrices preserving the norm\n",
    "\n",
    "- For stability it is really important that the error does not grow after we apply some transformations. \n",
    "\n",
    "- Suppose you are given $\\widehat{x}$ –– the approximation of $x$ such that,  \n",
    "\n",
    "$$ \\frac{\\Vert x - \\widehat{x} \\Vert}{\\Vert x \\Vert} \\leq \\varepsilon. $$\n",
    "\n",
    "- Let us calculate a linear transformation of $x$ and $\\widehat{x}$:  \n",
    "\n",
    "$$ y = U x, \\quad \\widehat{y} = U \\widehat{x}. $$\n",
    "\n",
    "- When building new algorithms, we want to use transformations that do not increase (or even preserve) the error:\n",
    "\n",
    "$$ \\frac{\\Vert y - \\widehat{y} \\Vert}{\\Vert y \\Vert } = \\frac{\\Vert U ( x - \\widehat{x}) \\Vert}{\\Vert U  x\\Vert}  \\leq \\varepsilon. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- The question is for which kind of matrices the norm of the vector **will not change**, so that\n",
    "\n",
    "$$ \\frac{\\Vert U ( x - \\widehat{x}) \\Vert}{\\Vert U  x\\Vert} = \\frac{ \\|x - \\widehat{x}\\|}{\\|x\\|}. $$\n",
    "\n",
    "- For the euclidean norm $\\|\\cdot\\|_2$ the answer is **unitary** (or orthogonal in the real case) matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Unitary (orthogonal) matrices\n",
    "\n",
    "- Let $U$ be complex $n \\times n$ matrix, and $\\Vert U z \\Vert_2 = \\Vert z \\Vert_2$ for all $z$. \n",
    "\n",
    "- This can happen **if and only if** (can be abbreviated as **iff**)\n",
    "\n",
    "$$ U^* U = I_n, $$\n",
    "\n",
    "where $I_n$ is an identity matrix $n\\times n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- Complex $n\\times n$ square matrix is called **unitary** if\n",
    "\n",
    "$$ U^*U = UU^* = I_n, $$\n",
    "\n",
    "which means that columns and rows of unitary matrices both form orthonormal basis in $\\mathbb{C}^{n}$.\n",
    "\n",
    "- For rectangular matrices of size $m\\times n$ ($n\\not= m$) only one of the equalities can hold\n",
    "\n",
    "    - $ U^*U = I_n$ –– left unitary for $m>n$\n",
    "    - $ UU^* = I_m$  –– right unitary for $m<n$\n",
    "\n",
    "- In the case of real matrices $U^* = U^T$ and matrices such that\n",
    "\n",
    "$$ U^TU = UU^T = I $$\n",
    "\n",
    "are called **orthogonal**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Unitary matrices\n",
    "\n",
    "Important property: **a product of two unitary matrices is a unitary matrix:**  \n",
    "\n",
    "$$(UV)^* UV = V^* (U^* U) V = V^* V = I,$$\n",
    "\n",
    "- Later we will show that there are types of matrices (**Householder reflections** and **Givens rotations**) composition of which is able to produce arbitrary unitary matrix\n",
    "- This idea is a core of some algorithms, e.g. QR decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Unitary invariance of $\\|\\cdot\\|_2$ and $\\|\\cdot\\|_F$ norms\n",
    "\n",
    "- For vector 2-norm we have already seen that $\\Vert U z \\Vert_2 = \\Vert z \\Vert_2$ for any unitary $U$.\n",
    "\n",
    "- One can show that unitary matrices also do not change matrix norms $\\|\\cdot\\|_2$ and $\\|\\cdot\\|_F$, i.e. for any square $A$ and unitary $U$,$V$: \n",
    "\n",
    "$$ \\| UAV\\|_2 = \\| A \\|_2 \\qquad \\| UAV\\|_F = \\| A \\|_F.$$\n",
    "\n",
    "- For $\\|\\cdot\\|_2$ it follows from the definition of an operator norm and the fact that vector 2-norm is unitary invariant.\n",
    "\n",
    "- For $\\|\\cdot\\|_F$ it follows from $\\|A\\|_F^2 = \\mathrm{trace}(A^*A)$ and the fact that $\\mathrm{trace}(BC) = \\mathrm{trace}(CB)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples of unitary matrices\n",
    "- There are two important classes of unitary matrices, using composition of which we can construct any unitary matrix:\n",
    "\n",
    "    1. Householder matrices\n",
    "    2. Givens (Jacobi) matrices\n",
    "\n",
    "- Other important examples are\n",
    "    * **Permutation matrix** $P$ whose rows (columns) are permutation of rows (columns) of the identity matrix.\n",
    "    * **Fourier matrix** $F_n = \\frac{1}{\\sqrt{n}} \\left\\{ e^{-i\\frac{2\\pi kl}{n}}\\right\\}_{k,l=0}^{n-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Householder matrices\n",
    "\n",
    "- Householder matrix is the matrix of the form \n",
    "\n",
    "$$H \\equiv H(v) = I - 2 vv^*,$$\n",
    "\n",
    "where $v$ is an $n \\times 1$ column and $v^* v = 1$. \n",
    "- Can you show that $H$ is unitary and Hermitian ($H^* = H$)?  \n",
    "- It is also a reflection:\n",
    "\n",
    "$$ Hx = x - 2(v^* x) v$$\n",
    "\n",
    "<img src=\"householder.jpeg\" width=500>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important property of Householder reflection\n",
    "\n",
    "- A nice property of Householder transformation is that it can zero all elements of a vector except for the first one:\n",
    "\n",
    "$$ H \\begin{bmatrix} \\times \\\\ \\times \\\\ \\times \\\\ \\times  \\end{bmatrix} =  \\begin{bmatrix} \\times \\\\ 0 \\\\ 0 \\\\ 0  \\end{bmatrix}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "*Proof (for real case).* Let $e_1 = (1,0,\\dots, 0)^T$, then we want to find $v$ such that\n",
    "\n",
    "$$ H x = x - 2(v^* x) v = \\alpha e_1, $$\n",
    "\n",
    "where $\\alpha$ is an unknown constant. Since $\\|\\cdot\\|_2$ is unitary invariant we get\n",
    "\n",
    "$$\\|x\\|_2 = \\|Hx\\|_2 = \\|\\alpha e_1\\|_2 = |\\alpha|.$$\n",
    "\n",
    "and $$\\alpha = \\pm \\|x\\|_2$$\n",
    "\n",
    "Also, we can express $v$ from $x - 2(v^* x) v = \\alpha e_1$:\n",
    "\n",
    "$$v = \\dfrac{x-\\alpha e_1}{2 v^* x}$$\n",
    "\n",
    "Multiplying the latter expression by $x^*$ we get\n",
    "\n",
    "$$x^* x - 2 (v^* x) x^* v = \\alpha x_1; $$\n",
    "\n",
    "or \n",
    "\n",
    "$$ \\|x\\|_2^2 - 2 (v^* x)^2 = \\alpha x_1. $$\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$ (v^* x)^2 = \\frac{\\|x\\|_2^2 - \\alpha x_1}{2}. $$\n",
    "\n",
    "So, $v$ exists and equals\n",
    "\n",
    "$$ v = \\dfrac{x \\mp \\|x\\|_2 e_1}{2v^* x} = \\dfrac{x \\mp \\|x\\|_2 e_1}{\\pm\\sqrt{2(\\|x\\|_2^2 \\mp \\|x\\|_2 x_1)}}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Householder algorithm for QR decomposition\n",
    "\n",
    "Using the obtained property we can make arbitrary matrix $A$ lower triangular:\n",
    "\n",
    "$$\n",
    "H_2 H_1 A = \\begin{bmatrix} \\times & \\times & \\times & \\times \\\\ 0 & \\times & \\times & \\times  \\\\  0 & 0 & \\boldsymbol{\\times} & \\times\\\\ 0 &0 & \\boldsymbol{\\times} & \\times  \\\\ 0 &0 & \\boldsymbol{\\times} & \\times \\end{bmatrix} $$\n",
    "\n",
    "then finding $H_3=\\begin{bmatrix}I_2 & \\\\ & {\\widetilde H}_3 \\end{bmatrix}$ such that\n",
    "\n",
    "$$ {\\widetilde H}_3 \\begin{bmatrix} \\boldsymbol{\\times}\\\\ \\boldsymbol{\\times} \\\\ \\boldsymbol{\\times}  \\end{bmatrix} = \\begin{bmatrix} \\times \\\\ 0 \\\\ 0  \\end{bmatrix}. $$\n",
    "\n",
    "we get\n",
    "\n",
    "$$ H_3 H_2 H_1 A =  \\begin{bmatrix} \\times & \\times & \\times & \\times \\\\  0 & \\times & \\times & \\times  \\\\  0 & 0 & {\\times} & \\times\\\\  0 &0 & 0 & \\times  \\\\  0 &0 & 0 & \\times  \\end{bmatrix} $$\n",
    "\n",
    "Finding $H_4$ by analogy we arrive at upper-triangular matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Since product and inverse of unitary matrices is a unitary matrix we get:\n",
    "\n",
    "**Corollary:** (QR decomposition) Every $A\\in \\mathbb{C}^{n\\times m}$ can be represented as\n",
    "\n",
    "$$ A = QR, $$\n",
    "\n",
    "where $Q$ is unitary and $R$ is upper triangular. \n",
    "\n",
    "See [poster](../../decompositions.pdf), what are the sizes of $Q$ and $R$ for $n>m$ and $n<m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Givens (Jacobi) matrix\n",
    "\n",
    "- A Givens matrix is a matrix  \n",
    "\n",
    "$$ G = \\begin{bmatrix} \\cos \\alpha & -\\sin \\alpha \\\\ \\sin \\alpha & \\cos \\alpha \\end{bmatrix},$$\n",
    "\n",
    "which is a rotation. \n",
    "\n",
    "- For a general case, we select two $(i, j)$ planes and rotate vector $x$  \n",
    "\n",
    "$$ x' = G x, $$\n",
    "\n",
    "only in the $i$-th and $j$-th positions:\n",
    "\n",
    "$$ x'_i =  x_i\\cos \\alpha - x_j\\sin \\alpha , \\quad x'_j = x_i \\sin \\alpha  +  x_j\\cos\\alpha, $$\n",
    "\n",
    "with all other $x_i$ remain unchanged.\n",
    "- Therefore, we can make elements in the $j$-th  position zero by choosing $\\alpha$ such that\n",
    "\n",
    "$$ \\cos \\alpha = \\frac{x_i}{\\sqrt{x_i^2 + x_j^2}}, \\quad \\sin \\alpha = -\\frac{x_j}{\\sqrt{x_i^2 + x_j^2}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "alpha = -3*jnp.pi / 4\n",
    "G = jnp.array([\n",
    "    [jnp.cos(alpha), -jnp.sin(alpha)],\n",
    "    [jnp.sin(alpha), jnp.cos(alpha)]\n",
    "])\n",
    "x = jnp.array([-1./jnp.sqrt(2), 1./jnp.sqrt(2)])\n",
    "y = G @ x\n",
    "\n",
    "plt.quiver([0, 0], [0, 0], [x[0], y[0]], [x[1], y[1]], angles='xy', scale_units='xy', scale=1)\n",
    "plt.xlim(-1., 1.)\n",
    "plt.ylim(-1., 1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR via Givens rotations\n",
    "\n",
    "Similarly we can make matrix upper-triangular using Givens rotations:\n",
    "\n",
    "$$\\begin{bmatrix} \\times & \\times & \\times \\\\ \\bf{*} & \\times & \\times \\\\ \\bf{*} & \\times & \\times \\end{bmatrix} \\to \\begin{bmatrix} * & \\times & \\times \\\\ * & \\times & \\times \\\\ 0 & \\times & \\times \\end{bmatrix} \\to \\begin{bmatrix} \\times & \\times & \\times \\\\ 0 & * & \\times \\\\ 0 & * & \\times \\end{bmatrix} \\to \\begin{bmatrix} \\times & \\times & \\times \\\\ 0 & \\times & \\times \\\\ 0 & 0 & \\times \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Givens vs. Householder transformations\n",
    "\n",
    "- Householder reflections are useful for dense matrices (complexity is $\\approx$ twice smaller than for Jacobi) and we need to zero large number of elements.\n",
    "- Givens rotations are more suitable for sparse matrice or parallel machines as they act locally on elements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular Value Decomposition\n",
    "\n",
    "SVD will be considered later in more details.\n",
    "\n",
    "**Theorem.** Any matrix $A\\in \\mathbb{C}^{n\\times m}$ can be written as a product of three matrices:  \n",
    "\n",
    "$$ A = U \\Sigma V^*, $$\n",
    "\n",
    "where \n",
    "- $U$ is an $n \\times n$ unitary matrix\n",
    "- $V$ is an $m \\times m$ unitary matrix\n",
    "- $\\Sigma$ is a diagonal matrix with non-negative elements $\\sigma_1 \\geq  \\ldots, \\geq \\sigma_{\\min (m,n)}$ on the diagonal.\n",
    "\n",
    "Moreover, if $\\text{rank}(A) = r$, then $\\sigma_{r+1} = \\dots = \\sigma_{\\min (m,n)} = 0$.\n",
    "\n",
    "See [poster](../../decompositions.pdf) for the visualization.\n",
    "\n",
    "- **Important note:** if one truncates (replace by $0$) all singular values except for $r$ first, then the resulting matrix yields best rank-$r$ approximation both in $\\|\\cdot\\|_2$ and $\\|\\cdot\\|_F$. \n",
    "- This is called Eckart-Young theorem and will be proved later in our course."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary \n",
    "\n",
    "- Most important matrix norms: Frobenius and spectral\n",
    "\n",
    "- Unitary matrices preserve these norms\n",
    "\n",
    "- There are two \"basic\" classes of unitary matrices: Householder and Givens matrices"
   ]
  },
  {
   "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": "Slideshow",
  "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"
  },
  "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
}