15_Probabilistische_Lineare_Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probabilistische Lineare Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grundlagen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aus der linearen Algebra sind viele Matrix-Zerlegungen bekannt, z.B. LU, QR, SVD.\n",
    "Wir betrachten hier die SVD, also\n",
    "  \\begin{equation*} \n",
    "  A = U\\Sigma V^T\n",
    "  \\end{equation*}\n",
    "mit $U,V$ orthonormal und $\\Sigma = \\mathrm{diag}_i(\\sigma_i)$.\n",
    "  \n",
    "Hat $A$ den Rang $r$, dann gilt\n",
    "  \\begin{align*}\n",
    "   \\sigma_i &> 0 \\quad i=1,\\ldots,r \\\\\n",
    "   \\sigma_i &= 0 \\quad i=r+1,\\ldots\n",
    "  \\end{align*}\n",
    "und man erhält\n",
    "  \\begin{equation*} \n",
    "  \\underbrace{A}_{\\mathbb{R}^{m \\times n}}=\n",
    "  \\underbrace{U_r}_{\\mathbb{R}^{m \\times r}}\n",
    "  \\ \\underbrace{\\Sigma_r}_{\\mathbb{R}^{r \\times r}}\n",
    "  \\ \\underbrace{V_r^T}_{\\mathbb{R}^{r \\times n}}\n",
    "  \\end{equation*}\n",
    "mit $r \\leq \\min(m,n)$.\n",
    "Ist $r \\ll \\min(m,n)$, dann können Operationen mit $A$, wie z.B.\n",
    "  Matrix-Vektor-Produkte, über die Zerlegung mit deutlich weniger\n",
    "  Aufwand berechnet werden.\n",
    "\n",
    "Die Spalten von $U_r$ bilden eine Orthonormalbasis des Bildraums\n",
    "  $R(A)$, d.h. sie enthalten die komplette Information über  $R(A)$.\n",
    "Gilt nun\n",
    "  \\begin{equation*} \n",
    "  \\sigma_1 \\geq \\ldots \\geq \\sigma_k \\geq \\sigma_{k+1}\\geq \\ldots \\geq \\sigma_r,\n",
    "  \\quad\n",
    "  \\sigma_1 \\gg \\sigma_{k+1}\n",
    "  \\end{equation*}\n",
    "dann hatten wir bei TSVD einfach $\\sigma_{k+1} = \\ldots = \\sigma_r =0$ gesetzt\n",
    "und damit eine Matrix-Approximation durchgeführt\n",
    "  \\begin{equation*} \n",
    "  A = U\\Sigma V^T \\approx U_k \\Sigma_k V_k^T.\n",
    "  \\end{equation*}\n",
    "Die Spalten von $U_k$ sind eine Orthonormalbasis eines Unterraums von $R(A)$.\n",
    "Wegen $\\sigma_1 \\gg \\sigma_{k+1}$ \"hoffen\" wir, das wir damit den \"wesentlichen Teil\" \n",
    "  von $R(A)$ erwischt haben und damit $U_k \\Sigma_k V_k^T$ eine gute Approximation von\n",
    "  $A$ ist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diese Vorgehensweise lässt sich verallgemeinern:\n",
    "\n",
    "- für $A\\in \\mathbb{R}^{m\\times n}$ sucht man eine Rang-$k$-Approximation\n",
    "  \\begin{equation*} \n",
    "  A \\approx QC, \\quad Q \\in \\mathbb{R}^{m\\times k}, \\quad C \\in \\mathbb{R}^{k\\times n}\n",
    "  \\end{equation*}\n",
    "  \n",
    "- dazu konstruiert man zunächst einen Unterraum von $R(A)$ der möglichst kleine \n",
    "  Dimension $k$ haben sollte, aber gleichzeitig die wesentlichen Abbildungseigenschaften \n",
    "  von $A$ wiedergeben soll\n",
    "  \n",
    "- eine Orthonormalbasis dieses Unterraums legen wir als Spalten in der Matrix\n",
    "  \\begin{equation*} \n",
    "  Q\\in\\mathbb{R}^{m\\times k}\n",
    "  \\end{equation*}\n",
    "  ab und approximieren $A$ durch Orthogonalprojektion in diesen Unterraum, d.h.\n",
    "  \\begin{equation*} \n",
    "  A \\approx QQ^T A =  QC = B, \\quad C = Q^T A \\in \\mathbb{R}^{k\\times n}\n",
    "  \\end{equation*}\n",
    "  \n",
    "- mit $Q$ erzeugen wir also eine niedrig dimensionale \n",
    "  Matrix $C = Q^TA$ so dass $QC$ die Ausgangsmatrix $A$\n",
    "  approximiert\n",
    "    \n",
    "- auf $C$ wenden wir dann Standard-Zerlegungen (LU, QR, SVD) \n",
    "  an, was aufgrund der niedrigen Dimension wenig Aufwand verursacht,\n",
    "  und konstruieren daraus approximative Zerlegungen für \n",
    "  die Ausgangs-Matrix $A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir gehen also zweistufig vor:\n",
    "\n",
    "  - bestimme $Q$\n",
    "  \n",
    "  - berechne die Zerlegung von $C = Q^T A$ und darüber\n",
    "    eine Zerlegung von $B = QC$, was dann eine \n",
    "    approximative Zerlegung von $A$ darstellt\n",
    "  \n",
    "Für den ersten Teil werden wir probabilistische Methoden benutzen.\n",
    "\n",
    "\n",
    "**Beispiel:** Approximative SVD für $A\\in \\mathbb{R}^{m \\times n}$\n",
    "\n",
    "- $Q_k\\in \\mathbb{R}^{m \\times k}$ orthonormal sei gegeben, $k\\leq\\min(m,n)$\n",
    "  \n",
    "- ist $C_k = Q_k^TA \\in \\mathbb{R}^{k \\times n}$, dann ist\n",
    "  \\begin{equation*} \n",
    "  B_k = Q_kC_k = Q_k Q_k^TA \\approx A\n",
    "  \\end{equation*}\n",
    "  eine Rang-$k$-Approximation von $A$\n",
    "  \n",
    "- wir berechnen die SVD von $C_k$ und erhalten\n",
    "  \\begin{equation*} \n",
    "  C_k = \\tilde{U}_k \\Sigma_k V_k^T,\n",
    "  \\quad\n",
    "  \\tilde{U}_k\\in \\mathbb{R}^{k \\times k},\n",
    "  \\quad\n",
    "  \\Sigma_k\\in \\mathbb{R}^{k \\times k},\n",
    "  \\quad\n",
    "  V_k^T \\in \\mathbb{R}^{k \\times n}\n",
    "  \\end{equation*}\n",
    "  \n",
    "- mit $A\\approx B_k = Q_k C_k$ folgt dann\n",
    "  \\begin{equation*} \n",
    "  A\\approx Q_k C_k  = \\underbrace{Q_k \\tilde{U}_k}_{U_k} \\Sigma_k V_k^T,\n",
    "  \\end{equation*}\n",
    "  mit $U_k\\in \\mathbb{R}^{m \\times k}$ orthonormal\n",
    "  \n",
    "  \n",
    "Die offene Frage ist jetzt, wie wir geeignete Matrizen $Q_k$ finden können."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probabilistische Verfahren"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In diesem Abschnitt betrachten wir die Konstruktion von $Q_k$.\n",
    "\n",
    "Ab jetzt lassen wir den Index $k$ weg.\n",
    "Wir müssen folgendes Problem lösen:\n",
    "\n",
    "- gegeben ist $A\\in \\mathbb{R}^{m \\times n}$, $\\varepsilon > 0$\n",
    "  \n",
    "- finde eine orthonormale Matrix $Q\\in \\mathbb{R}^{m \\times k}$,\n",
    "    $k = k(\\varepsilon) \\leq \\min(m,n)$, mit\n",
    "    \\begin{equation*} \n",
    "    \\|A -QQ^TA\\| \\leq \\varepsilon,\n",
    "    \\end{equation*}\n",
    "  d.h.\n",
    "    \\begin{equation*} \n",
    "    \\|A - X\\| \\leq \\varepsilon\n",
    "    \\end{equation*}\n",
    "  mit $X=QQ^TA$ mit $\\mathrm{rang}(X) \\leq k$\n",
    "    \n",
    "$k(\\varepsilon)$ sollte so klein wie möglich sein.\n",
    "Als Matrix-Norm benutzen wir in der Regel\n",
    "  \\begin{align*} \n",
    "  \\|A\\|_2 \n",
    "  &= \\sqrt{\\rho(A^TA)},\n",
    "  \\\\\n",
    "  \\|A\\|_F \n",
    "  &= \\sqrt{\\mathrm{tr}(A^TA)} = \\sqrt{\\sum_{i=1}^m\\sum_{j=1}^n a_{ij}^2}.\n",
    "  \\end{align*}\n",
    "Für orthonormal-invariante Matrix-Normen, d.h.\n",
    "  \\begin{equation*} \n",
    "  \\|A\\| = \\|QA\\| = \\|A\\tilde{Q}\\|\n",
    "  \\quad\n",
    "  \\forall Q,\\tilde{Q} \\ \\text{orthonormal}\n",
    "  \\end{equation*}\n",
    "kennt man eine Optimal-Lösung dieses Problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Satz (Eckart, Young, Mirsky):** Ist $A\\in \\mathbb{R}^{m \\times n}$ \n",
    "mit SVD $A=U\\Sigma V^T$ und die Matrixnorm $\\|\\cdot\\|$ orthonormal invariant, dann\n",
    "gilt\n",
    "  \\begin{equation*} \n",
    "  \\min_{X\\in \\mathbb{R}^{m \\times n}, \\mathrm{rang}(X)\\leq k} \\|A-X\\| = \\sigma_{k+1}.\n",
    "  \\end{equation*}\n",
    "Das Minimum wird für\n",
    "  \\begin{equation*} \n",
    "  X_{\\min} = U_k \\Sigma_k V_k^T,\n",
    "  \\end{equation*}\n",
    "also TSVD, angenommen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Normen $\\|\\cdot\\|_2$ und $\\|\\cdot\\|_F$ sind orthonormal \n",
    "invariant (Übung), somit können wir das Ergebnis hier in anwenden.\n",
    "In unserem Set-Up bedeutet das:\n",
    "\n",
    "- in $X_{\\min}$ kommen nur die Spalten $u_1,\\ldots,u_k$ von $U$ vor, d.h.\n",
    "    \\begin{equation*} \n",
    "    Q = U_k = (u_1,\\ldots,u_k) \\in \\mathbb{R}^{m\\times k}\n",
    "    \\end{equation*}\n",
    "    \n",
    "- somit ist\n",
    "  \\begin{align*}\n",
    "  QQ^TA \n",
    "  &= U_k U_k^T U\\Sigma V^T \\\\\n",
    "  &= U_k (I_k,0)\\Sigma V^T \\\\\n",
    "  &= U_k (\\Sigma_k,0) V^T \\\\\n",
    "  &= U_k \\Sigma_k V_k^T \n",
    "  \\end{align*}\n",
    "  und wir könnten wie folgt vorgehen:\n",
    "  \n",
    "    - gebe $\\varepsilon > 0$ vor\n",
    "    \n",
    "    - suche $k$ mit $\\sigma_{k+1}<\\varepsilon$\n",
    "    \n",
    "    - benutze $Q=U_k$ und \n",
    "    \\begin{equation*} \n",
    "    B = QQ^TA = U_k \\Sigma_k V_k^T \n",
    "    \\end{equation*}\n",
    "    als Rang-$k$-Approximation\n",
    "    \n",
    "Das Problem dabei ist, dass man dazu eine komplette SVD von $A$ benötigt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir müssen also überlegen, wie wir \"billiger\" an $Q$ kommen können.\n",
    "Dazu betrachten wir für $A\\in \\mathbb{R}^{m \\times n}$ zwei Fälle:\n",
    "\n",
    "- fixed-precision: \n",
    "  gegeben $\\varepsilon > 0$, \n",
    "  suche $Q$ orthonormal mit \n",
    "  $\\mathrm{rang}(Q)\\leq k(\\varepsilon)$,  \n",
    "  wobei\n",
    "  $k(\\varepsilon)$ möglichst klein sein soll,\n",
    "  und \n",
    "  \\begin{equation*} \n",
    "  \\|A -QQ^TA\\| \\leq \\varepsilon\n",
    "  \\end{equation*}\n",
    "  \n",
    "- fixed-rank: $k,p$ ist gegeben, suche $Q \\in \\mathbb{R}^{m\\times (k+p)}$ orthonormal mit \n",
    "  $\\mathrm{rang}(Q)\\leq k$\n",
    "  und \n",
    "  \\begin{equation*} \n",
    "  \\|A -QQ^TA\\| \n",
    "  \\approx\n",
    "  \\min_{X\\in \\mathbb{R}^{m \\times n}, \\mathrm{rang}(X)\\leq k} \\|A-X\\| \n",
    "  = \n",
    "  \\sigma_{k+1}\n",
    "  \\end{equation*}\n",
    "  \n",
    "Bei fixed-rank lassen wir bei $Q$ mehr Spalten zu als für eine Rang-$k$-Approximation \n",
    "  unbedingt nötig wären ($p$).\n",
    "Dieser Freiheitsgrad wird später eine wichtige Rolle bei der Genauigkeit\n",
    "  probabilistischer Methoden spielen.\n",
    "  \n",
    "Wir untersuchen zunächst fixed-rank und knacken damit dann später auch fixed-precision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed-Rank"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wie kann der \"Zufall\" bei der Konstruktion von $Q$ helfen?\n",
    "\n",
    "Wir betrachten dazu zunächst $A\\in \\mathbb{R}^{m \\times n}$ mit $\\mathrm{rang}(A)=k$:\n",
    "\n",
    "- ziehe $k$ Zufallsvektoren $\\omega^{(i)}$ (Verteilung zunächst uninteressant)\n",
    "  und berechne \n",
    "  \\begin{equation*} \n",
    "  y^{(i)} = A \\omega^{(i)}, \\quad i=1,\\ldots,k\n",
    "  \\end{equation*}\n",
    "  \n",
    "- mit hoher Wahrscheinlichkeit gilt:\n",
    "  \n",
    "  - die $\\omega^{(i)}$ sind linear unabhängig\n",
    "    \n",
    "  - $\\omega^{(i)} \\not\\in N(A)$\n",
    "    \n",
    "- damit sind dann auch die $y^{(i)}$ mit hoher Wahrscheinlichkeit linear unabhängig\n",
    "  \n",
    "- wegen $\\mathrm{rang}(A)=k$ gilt dann\n",
    "  \\begin{equation*} \n",
    "  R(A)=\\mathrm{span}\\big(y^{(1)},\\ldots,y^{(k)}\\big),\n",
    "  \\end{equation*}\n",
    "  d.h. $Y = \\big(y^{(1)},\\ldots,y^{(k)}\\big)$ spannt $R(A)$ auf\n",
    "  \n",
    "- orthonormalisieren wir jetzt $Y$ (modifizierter Gram-Schmidt, QR, etc.), so erhalten wir\n",
    "  $Q\\in \\mathbb{R}^{m \\times k}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jetzt sei $A = A_0 + E$, $\\mathrm{rang}(A_0)=k$, wobei $E$ sei eine \"kleine\" Störung ist.\n",
    "  \n",
    "  - wir suchen einen Unterraum in $R(A)$, der einen möglichst \"großen Teil\" von\n",
    "    $R(A_0)$ abdeckt, aber nicht unbedingt die kleinstmögliche Dimension $k$ hat\n",
    "    sondern Dimension $k+p$\n",
    "    \n",
    "  - wir betrachten wieder\n",
    "    \\begin{equation*} \n",
    "    y^{(i)} = A \\omega^{(i)} = A_0 \\omega^{(i)} + E  \\omega^{(i)}\n",
    "    \\end{equation*}\n",
    "  \n",
    "  - der \"Störanteil\" $E  \\omega^{(i)}$ \"schiebt\" $A_0 \\omega^{(i)}$ eventuell\n",
    "    aus $R(A_0)$ hinaus\n",
    "  \n",
    "  - wenn wir genau $k$ verschiedene $y^{(i)}$ benutzen, dann ist die Wahrscheinlichkeit \n",
    "    hoch, dass wir \"wichtige Teile\" von $R(A_0)$ nicht erfassen\n",
    "    \n",
    "  - deshalb benutzen wir $k+p$ Vektoren $y^{(i)}$, $p\\in\\mathbb{N}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fassen wir diese Schritte zusammen, so erhalten wir das **Fixed-Rank-Verfahren**:\n",
    "\n",
    "  - gegeben ist $A\\in \\mathbb{R}^{m \\times n}$, $k,p\\in \\mathbb{N}$\n",
    "  \n",
    "  - bestimme eine \"zufällige\" Matrix $\\Omega \\in \\mathbb{R}^{n \\times (k+p)}$\n",
    "  \n",
    "  - berechne damit\n",
    "  \\begin{equation*} \n",
    "  Y = A\\Omega \\in \\mathbb{R}^{m \\times (k+p)}\n",
    "  \\end{equation*}\n",
    "  \n",
    "  - orthonormiere die Spalten von $Y$ und erzeuge damit die orthonormale Matrix\n",
    "  \\begin{equation*} \n",
    "  Q \\in \\mathbb{R}^{m \\times (k+p)}\n",
    "  \\end{equation*}\n",
    "  \n",
    "  \n",
    "Dabei sind folgende Fragen noch offen:\n",
    "\n",
    "- welches $\\Omega$ und welches $p$ soll benutzt werden?\n",
    "  \n",
    "- wie orthonormalisiert man $Y$, das in der Regel schlechte Kondition hat\n",
    "  (Spalten sind \"fast\" linear abhängig)?\n",
    "  \n",
    "- wie ist der Aufwand, die Genauigkeit?\n",
    "  \n",
    "- wie geht fixed-precision?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Satz:** $A \\in \\mathbb{R}^{m \\times n}$, $k\\geq 2$, $p\\geq 2$, $k+p \\leq \\min(m,n)$,\n",
    "\\begin{equation*} \n",
    "\\Omega = \\big(\\omega_{ij}\\big)_{i,j} \\in  \\mathbb{R}^{n \\times (k+p)},\n",
    "\\quad\n",
    "\\omega_{ij} \\sim \\mathcal{N}(0,1) \\ \\text{iid},  \n",
    "\\end{equation*}\n",
    "dann gilt\n",
    "\\begin{align*} \n",
    "\\mathbb{E}_\\Omega\\big(\\|A - QQ^TA\\|_2\\big) \n",
    "&\\leq (1+\\delta)\\sigma_{k+1},\n",
    "\\\\\n",
    "\\delta \n",
    "&= \\frac{4\\sqrt{k+p}}{p-1}\\sqrt{\\min(m,n)}\n",
    "\\end{align*}\n",
    "bzw.\n",
    "\\begin{align*} \n",
    "\\mathbb{P}\\big(\\|A - QQ^TA\\|_2 \n",
    "\\leq (1+\\tilde{\\delta})\\sigma_{k+1}\\big)\n",
    "&\\geq\n",
    "1 - 3p^{-p},\n",
    "\\\\\n",
    "\\tilde{\\delta} \n",
    "&= 9 \\sqrt{k+p} \\sqrt{\\min(m,n)}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bemerkung:**\n",
    "\n",
    "- $\\min_{X\\in \\mathbb{R}^{m \\times n}, \\mathrm{rang}(X)\\leq k} \\|A-X\\|_2 = \\sigma_{k+1}$ \n",
    "  so dass $\\delta$, $\\tilde{\\delta}$ die Abweichung vom Optimum angibt\n",
    "    \n",
    "- $\\delta$, $\\tilde{\\delta}$ wächst nur schwach in $k,m,n$\n",
    "  \n",
    "- es ist (unabhängig von $k,m,n$)\n",
    "  \\begin{equation*} \n",
    "  3 p^{-p}\n",
    "  \\approx\n",
    "  \\begin{cases}\n",
    "  10^{-3} & p=5 \\\\\n",
    "  3\\cdot 10^{-10} & p=10 \n",
    "  \\end{cases},\n",
    "  \\end{equation*}\n",
    "  so dass für recht kleine $p$-Werte\n",
    "  $\\mathbb{P}\\big(\\|A - QQ^TA\\|_2 \\leq (1+\\tilde{\\delta})\\sigma_{k+1}\\big)\\approx 1$\n",
    "  \n",
    "- fallen die Singulärwerte schnell ab, so dass $\\sigma_{k+1}$ klein ist, dann\n",
    "  sind die Approximationen sehr gut"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomized SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir wenden jetzt die Methoden des letzten Abschnitts, um eine approximative SVD\n",
    "  von $A$ zu berechnen.\n",
    "  \n",
    "Eine direkte Anwendung liefert oft schlechte Ergebnisse, insbesondere wenn die \n",
    "  Singulärwerte von $A$ nur langsam abfallen.\n",
    "Statt $Y=A\\Omega$ benutzt man dann die Power-Methode \n",
    "  \\begin{equation*} \n",
    "  Y = (AA^T)^q A\\Omega,\n",
    "  \\end{equation*}\n",
    "wobei in der Praxis meist $q=1$ oder $q=2$ ausreicht.\n",
    "  \n",
    "Es gilt $R\\big((AA^T)^q A\\big)\\subset R(A)$ und mit $A=U\\Sigma V^T$ folgt\n",
    "  \\begin{equation*} \n",
    "  AA^T = U\\Sigma\\Sigma^T U^T\n",
    "  \\end{equation*}\n",
    "  also\n",
    "  \\begin{align*} \n",
    "  (AA^T)^q A \n",
    "  &= (U\\Sigma\\Sigma^T U^T)^q U\\Sigma V^T\\\\\n",
    "  &= U(\\Sigma\\Sigma^T)^q\\Sigma V^T\\\\\n",
    "  &= U\\Sigma_{q+1} V^T\n",
    "  \\end{align*}\n",
    "  mit\n",
    "  \\begin{equation*} \n",
    "  \\quad\n",
    "  \\Sigma_{q+1} = \\mathrm{diag}(\\sigma_i^{q+1}),\n",
    "  \\end{equation*}\n",
    "  wobei $\\sigma_i^{q+1}$ für $q>0$ schneller abfällt als $\\sigma_i$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Damit erhalten wir schließlich folgenden Algorithmus für RSVD:\n",
    "\n",
    "- $A \\in \\mathbb{R}^{m \\times n}$, $k\\in\\mathbb{N}$, $p,q\\in\\mathbb{N}_0$\n",
    "  \n",
    "- Schritt 1:\n",
    "  \n",
    "  - $\\Omega = \\big(\\omega_{ij}\\big)_{i,j} \\in \\mathbb{R}^{n \\times (k+p)}$,\n",
    "    $\\omega_{ij} \\sim \\mathcal{N}(0,1)$ iid  \n",
    "\n",
    "  - $Y = (AA^T)^q A\\Omega \\in \\mathbb{R}^{m \\times (k+p)}$\n",
    "    \n",
    "  - orthonormalisiere die Spalten von $Y$ und lege das Ergebnis in\n",
    "    $Q\\in \\mathbb{R}^{m \\times (k+p)}$ ab\n",
    "  \n",
    "- Schritt 2:\n",
    "    \n",
    "  - $C = Q^T A\\in \\mathbb{R}^{(k+p) \\times n}$\n",
    "    \n",
    "  - berechne SVD von $C$, $C = \\tilde{U}\\Sigma V^T$\n",
    "    \n",
    "  - setze $U = Q\\tilde{U}$\n",
    "    \n",
    "- insgesamt ist dann $A \\approx U\\Sigma V^T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bemerkung:**\n",
    "\n",
    "- man berechnet eine Rang-$(k+p)$-Approximation um eine (gute) Näherung für die\n",
    "  ersten $k$ Singulärwerte und -vektoren zu erhalten\n",
    "   \n",
    "- die Produkte $AA^TA\\Omega$ sind anfällig gegen Rundungsfehler, weshalb man nach\n",
    "  jedem Schritt orthonormalisieren sollte\n",
    "   \n",
    "   \n",
    "**Satz:** $A \\in \\mathbb{R}^{m \\times n}$, $q\\in\\mathbb{N}_0$, \n",
    "$2 \\leq k \\leq \\frac{1}{2}\\min(m,n)$, dann gilt für $U,\\Sigma,V$ aus der RSVD\n",
    "mit $p=k$\n",
    "  \\begin{align*} \n",
    "  \\mathbb{E}_\\Omega\\big(\\|A - U\\Sigma V^T\\|_2\\big) \n",
    "  &\\leq \n",
    "  (1+\\delta)^{\\frac{1}{2q+1}}\\sigma_{k+1},\n",
    "  \\\\\n",
    "  \\delta &= 4 \\sqrt{\\frac{2\\min(m,n)}{k-1}}\n",
    "  \\end{align*}  \n",
    "\n",
    "\n",
    "Da wir nur an den ersten $k$ Approximationen der Singulärwerte und -vektoren\n",
    "interessiert sind, ist es naheliegend, die $\\sigma_i$ für $i>k$ abzuschneiden.\n",
    "Wir ersetzen oben $\\Sigma$ durch \n",
    "\\begin{equation*} \n",
    "\\Sigma_k =  \\mathrm{diag}(\\sigma_1,\\ldots,\\sigma_k,0,\\ldots,0),\n",
    "\\end{equation*}\n",
    "d.h. wir benutzen\n",
    "eine TSVD von $C$ in Schritt 2.\n",
    "In diesem Fall ändert sich die Fehlerabschätzung des letzten Satzes geringfügig.\n",
    "\n",
    "\n",
    "**Satz:** $A \\in \\mathbb{R}^{m \\times n}$, $q\\in\\mathbb{N}_0$, \n",
    "$2 \\leq k \\leq \\frac{1}{2}\\min(m,n)$, dann gilt für $U,\\Sigma_k,V$ aus der RSVD\n",
    "mit $p=k$\n",
    "\\begin{align*} \n",
    "\\mathbb{E}_\\Omega\\big(\\|A - U\\Sigma_k V^T\\|_2\\big) \n",
    "&\\leq \n",
    "(1+\\delta)^{\\frac{1}{2q+1}}\\sigma_{k+1} + \\sigma_{k+1},\n",
    "\\\\\n",
    "\\delta &= 4 \\sqrt{\\frac{2\\min(m,n)}{k-1}}\n",
    "\\end{align*}\n",
    "  \n",
    "  \n",
    "Damit ist die fixed-rank Variante der RSVD soweit erledigt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed-Precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Das fixed-precision Problem könnte wie folgt auf die fixed-rank Variante\n",
    "zurückgeführt werden:\n",
    "  \n",
    "- starte mit fixed-rank mit kleinem $k$\n",
    "  \n",
    "- teste, ob $\\|A - QQ^TA\\|_2 \\leq \\varepsilon$\n",
    "  \n",
    "- wenn nicht, dann erhöhe $k$ und wiederhole den vorherigen Schritt\n",
    "  \n",
    "  \n",
    "Das Berechnen der Matrixnorm $\\|A - QQ^TA\\|_2$ ist in der Regel zu aufwendig,\n",
    "deshalb benötigt man dafür einen einfachen a-posteriori Fehlerindikator.\n",
    "Auch hier hilft uns wieder der \"Zufall\" weiter.\n",
    "\n",
    "\n",
    "**Lemma:** Sei $B\\in\\mathbb{R}^{m\\times n}$, $\\omega^{(j)}\\in\\mathbb{R}^n$,\n",
    "$j=1,\\ldots,r$,\n",
    "$\\omega_i^{(j)} \\sim \\mathcal{N}(0,1)$  iid und $\\alpha > 1$.\n",
    "Dann gilt\n",
    "\\begin{equation*} \n",
    "\\mathbb{P}\n",
    "\\big(\n",
    "\\|B\\|_2 \n",
    "\\leq \n",
    "\\alpha\\sqrt{\\frac{2}{\\pi}} \\ \\max_{j=1,\\ldots,r}\\|B\\omega^{(j)}\\|_2\n",
    "\\big)\n",
    "\\geq\n",
    "1 - \\alpha^{-r}.\n",
    "\\end{equation*}\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Für die Praxis bedeutet das:\n",
    "\n",
    "- erzeuge $\\omega^{(j)}$,  $j=1,\\ldots,r$ für $r$ klein\n",
    "  \n",
    "- mit $m_r = \\max_{j=1,\\ldots,r}\\|(A-QQ^TA)\\omega^{(j)}\\|_2$ und $\\alpha=10$ folgt\n",
    "  \\begin{equation*} \n",
    "  \\mathbb{P}\n",
    "  \\big(\n",
    "  \\|A-QQ^TA\\|_2 \n",
    "  \\leq \n",
    "  10\\sqrt{\\frac{2}{\\pi}} \\ m_r\n",
    "  \\big)\n",
    "  \\geq\n",
    "  1 - 10^{-r},\n",
    "  \\end{equation*}\n",
    "  so dass $r\\leq 10$ ausreichend ist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prinzipiell wenden wir beim diesem Fehlerschätzer die gleichen Operationen\n",
    "wie bei der Gaussian-Projection\n",
    "\\begin{equation*} \n",
    "Y = A\\Omega,\n",
    "\\quad\n",
    "\\Omega = \\big(\\omega_{ij}\\big)_{i,j} \\in  \\mathbb{R}^{n \\times (k+p)},\n",
    "\\quad\n",
    "\\omega_{ij} \\sim \\mathcal{N}(0,1) \\ \\text{iid}\n",
    "\\end{equation*}\n",
    "an, nur dass wir statt der Matrix $\\Omega$ einzelne Vektoren $\\omega^{(j)}$\n",
    "benutzen.\n",
    "  \n",
    "Ordnen wir die Berechnungsschritte etwas um, dann können wir den Fehlerschätzer\n",
    "praktisch ohne Mehraufwand in unser fixed-rank-Verfahren integrieren:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- gegeben sind $A\\in\\mathbb{R}^{m\\times n}$, $r\\in\\mathbb{N}$, $\\varepsilon>0$:\n",
    "\n",
    "  - erzeuge $\\omega^{(j)}\\in\\mathbb{R}^n$, $j=1,\\ldots,r$,\n",
    "    $\\omega_i^{(j)} \\sim \\mathcal{N}(0,1)$  iid \n",
    "    \n",
    "  - $y^{(j)} = A \\omega^{(j)}$, $j=1,\\ldots,r$\n",
    "\n",
    "  - $k=0$, $Q^{(0)} = [\\ ]$ (leer)\n",
    "  \n",
    "  - wiederhole bis \n",
    "    $\\displaystyle \n",
    "    \\max_{j=1,\\ldots,r}\\|y^{(k+j)}\\|_2 \\leq \\frac{\\varepsilon}{10\\sqrt{\\frac{2}{\\pi}}}$:\n",
    "    \n",
    "    - $k := k+1$\n",
    "    \n",
    "    - $y^{(k)} := \\big(I - Q^{(k-1)}{Q^{(k-1)}}^T\\big) y^{(k)}$\n",
    "    \n",
    "    - $q^{(k)} := \\frac{y^{(k)}}{\\|y^{(k)}\\|_2}$\n",
    "    \n",
    "    - $Q^{(k)} = [Q^{(k-1)}, q^{(k)}]$\n",
    "    \n",
    "    - erzeuge $\\omega^{(k+r)} \\in \\mathbb{R}^n$, \n",
    "      $\\omega^{(k+r)} \\sim \\mathcal{N}(0,1)$ iid\n",
    "      \n",
    "    - $y^{(k+r)} := \\big(I - Q^{(k)}{Q^{(k)}}^T\\big) A \\omega^{(k+r)}$\n",
    "    \n",
    "    - $y^{(l)} := y^{(l)} - {q^{(k)}}^T y^{(l)} {q^{(k)}}$ für $l = k+1,k+2,\\ldots,k+r-1$\n",
    "    \n",
    "  - mit $Q = Q^{(j)}$ gilt dann \n",
    "  \\begin{equation*} \n",
    "  \\mathbb{P}\\big(\\|A-QQ^TA\\|_2 \\leq \\varepsilon\\big) \\geq 1- 10^{-r}\\min(m,n)\n",
    "  \\end{equation*}\n",
    "    \n",
    "Damit ist $Q$ bestimmt.\n",
    "Der Rest der RSVD ist identisch zum fixed-rank-Fall, insbesondere\n",
    "kann dabei auch wieder die Power-Methode benutzt werden."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zusammenfassung"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir haben probabilistische Methoden für approximative Berechnung von TSVD\n",
    "betrachtet und Algorithmen für fixed-rank und fixed-precision Aufgabenstellungen\n",
    "vorgestellt."
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.9.6"
  },
  "nbTranslate": {
   "displayLangs": [
    "fr",
    "en"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": false,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}