{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 1.1.0\n", "Commit 80516ca202 (2019-01-21 21:24 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin14.5.0)\n", " CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n", "Environment:\n", " JULIA_EDITOR = code\n" ] } ], "source": [ "versioninfo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Conjugate Gradient and Krylov Subspace Methods\n", "\n", "## Introduction\n", "\n", "* Conjugate gradient is the top-notch iterative method for solving large, **structured** linear systems $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{A}$ is pd. \n", "Earlier we talked about Jacobi, Gauss-Seidel, and successive over-relaxation (SOR) as the classical iterative solvers. They are rarely used in practice due to slow convergence. \n", "\n", " [Kershaw's results](http://www.sciencedirect.com/science/article/pii/0021999178900980?via%3Dihub) for a fusion problem.\n", "\n", "| Method | Number of Iterations |\n", "|----------------------------------------|----------------------|\n", "| Gauss Seidel | 208,000 |\n", "| Block SOR methods | 765 |\n", "| Incomplete Cholesky **conjugate gradient** | 25 |\n", "\n", "\n", "* History: Hestenes (**UCLA** professor!) and Stiefel proposed conjugate gradient method in 1950s.\n", "\n", "Hestenes and Stiefel (1952), [Methods of conjugate gradients for solving linear systems](http://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_A1b.pdf), _Jounral of Research of the National Bureau of Standards_.\n", "\n", "* Solve linear equation $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ is **pd**, is equivalent to \n", "$$\n", "\\begin{eqnarray*}\n", "\t\\text{minimize} \\,\\, f(\\mathbf{x}) = \\frac 12 \\mathbf{x}^T \\mathbf{A} \\mathbf{x} - \\mathbf{b}^T \\mathbf{x}.\n", "\\end{eqnarray*}\n", "$$\n", "Denote $\\nabla f(\\mathbf{x}) = \\mathbf{A} \\mathbf{x} - \\mathbf{b} =: r(\\mathbf{x})$.\n", "\n", "## Conjugate gradient (CG) method\n", "\n", "* Consider a simple idea: coordinate descent, that is to update components $x_j$ alternatingly. Same as the Gauss-Seidel iteration. Usually it takes too many iterations.\n", "\n", "\n", "\n", "* A set of vectors $\\{\\mathbf{x}^{(0)},\\ldots,\\mathbf{x}^{(l)}\\}$ is said to be **conjugate with respect to $\\mathbf{A}$** if\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\mathbf{p}_i^T \\mathbf{A} \\mathbf{p}_j = 0, \\quad \\text{for all } i \\ne j.\n", "\\end{eqnarray*}\n", "$$\n", "For example, eigen-vectors of $\\mathbf{A}$ are conjugate to each other. Why?\n", "\n", "* **Conjugate direction** method: Given a set of conjugate vectors $\\{\\mathbf{p}^{(0)},\\ldots,\\mathbf{p}^{(l)}\\}$, at iteration $t$, we search along the conjugate direction $\\mathbf{p}^{(t)}$\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\mathbf{x}^{(t+1)} = \\mathbf{x}^{(t)} + \\alpha^{(t)} \\mathbf{p}^{(t)},\n", "\\end{eqnarray*}\n", "$$\n", "where\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\alpha^{(t)} = - \\frac{\\mathbf{r}^{(t)T} \\mathbf{p}^{(t)}}{\\mathbf{p}^{(t)T} \\mathbf{A} \\mathbf{p}^{(t)}}\n", "\\end{eqnarray*}\n", "$$\n", "is the optimal step length.\n", "\n", "* Theorem: In conjugate direction method, $\\mathbf{x}^{(t)}$ converges to the solution in **at most** $n$ steps.\n", "\n", " Intuition: Look at graph.\n", " \n", "\n", "\n", "* **Conjugate gradient** method. Idea: generate $\\mathbf{p}^{(t)}$ using only $\\mathbf{p}^{(t-1)}$\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\mathbf{p}^{(t)} = - \\mathbf{r}^{(t)} + \\beta^{(t)} \\mathbf{p}^{(t-1)},\n", "\\end{eqnarray*}\n", "$$\n", "where $\\beta^{(t)}$ is determined by the conjugacy condition $\\mathbf{p}^{(t-1)T} \\mathbf{A} \\mathbf{p}^{(t)} = 0$\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\beta^{(t)} = \\frac{\\mathbf{r}^{(t)T} \\mathbf{A} \\mathbf{p}^{(t-1)}}{\\mathbf{p}^{(t-1)T} \\mathbf{A} \\mathbf{p}^{(t-1)}}.\n", "\\end{eqnarray*}\n", "$$\n", "\n", "* **CG algorithm (preliminary version)**: \n", "\n", " 0. Given $\\mathbf{x}^{(0)}$\n", " 0. Initialize: $\\mathbf{r}^{(0)} \\gets \\mathbf{A} \\mathbf{x}^{(0)} - \\mathbf{b}$, $\\mathbf{p}^{(0)} \\gets - \\mathbf{r}^{(0)}$, $t=0$\n", " 0. While $\\mathbf{r}^{(0)} \\ne \\mathbf{0}$\n", " 1. $\\alpha^{(t)} \\gets - \\frac{\\mathbf{r}^{(t)T} \\mathbf{p}^{(t)}}{\\mathbf{p}^{(t)T} \\mathbf{A} \\mathbf{p}^{(t)}}$\n", " 2. $\\mathbf{x}^{(t+1)} \\gets \\mathbf{x}^{(t)} + \\alpha^{(t)} \\mathbf{p}^{(t)}$\n", " 3. $\\mathbf{r}^{(t+1)} \\gets \\mathbf{A} \\mathbf{x}^{(t+1)} - \\mathbf{b}$\n", " 4. $\\beta^{(t+1)} \\gets \\frac{\\mathbf{r}^{(t+1)T} \\mathbf{A} \\mathbf{p}^{(t)}}{\\mathbf{p}^{(t)T} \\mathbf{A} \\mathbf{p}^{(t)}}$\n", " 5. $\\mathbf{p}^{(t+1)} \\gets - \\mathbf{r}^{(t+1)} + \\beta^{(t+1)} \\mathbf{p}^{(t)}$\n", " 6. $t \\gets t+1$\n", " \n", " Remark: The initial conjugate direction $\\mathbf{p}^{(0)} \\gets - \\mathbf{r}^{(0)}$ is crucial.\n", " \n", "* Theorem: With CG algorithm\n", " 0. $\\mathbf{r}^{(t)}$ are mutually orthogonal. \n", " 0. $\\{\\mathbf{r}^{(0)},\\ldots,\\mathbf{r}^{(t)}\\}$ is contained in the **Krylov subspace** of degree $t$ for $\\mathbf{r}^{(0)}$, denoted by\n", " $$\n", " \\begin{eqnarray*}\n", " {\\cal K}(\\mathbf{r}^{(0)}; t) = \\text{span} \\{\\mathbf{r}^{(0)},\\mathbf{A} \\mathbf{r}^{(0)}, \\mathbf{A}^2 \\mathbf{r}^{(0)}, \\ldots, \\mathbf{A}^{t} \\mathbf{r}^{(0)}\\}.\n", " \\end{eqnarray*}\n", " $$\n", " 0. $\\{\\mathbf{p}^{(0)},\\ldots,\\mathbf{p}^{(t)}\\}$ is contained in ${\\cal K}(\\mathbf{r}^{(0)}; t)$. \n", " 0. $\\mathbf{p}^{(0)}, \\ldots, \\mathbf{p}^{(t)}$ are conjugate with respect to $\\mathbf{A}$. \n", "The iterates $\\mathbf{x}^{(t)}$ converge to the solution in at most $n$ steps.\n", "\n", "* **CG algorithm (economical version)**: saves one matrix-vector multiplication.\n", "\n", " 0. Given $\\mathbf{x}^{(0)}$\n", " 0. Initialize: $\\mathbf{r}^{(0)} \\gets \\mathbf{A} \\mathbf{x}^{(0)} - \\mathbf{b}$, $\\mathbf{p}^{(0)} \\gets - \\mathbf{r}^{(0)}$, $t=0$\n", " 0. While $\\mathbf{r}^{(0)} \\ne \\mathbf{0}$\n", " 1. $\\alpha^{(t)} \\gets \\frac{\\mathbf{r}^{(t)T} \\mathbf{r}^{(t)}}{\\mathbf{p}^{(t)T} \\mathbf{A} \\mathbf{p}^{(t)}}$\n", " 2. $\\mathbf{x}^{(t+1)} \\gets \\mathbf{x}^{(t)} + \\alpha^{(t)} \\mathbf{p}^{(t)}$\n", " 3. $\\mathbf{r}^{(t+1)} \\gets \\mathbf{r}^{(t)} + \\alpha^{(t)} \\mathbf{A} \\mathbf{p}^{(t)}$\n", " 4. $\\beta^{(t+1)} \\gets \\frac{\\mathbf{r}^{(t+1)T} \\mathbf{r}^{(t+1)}}{\\mathbf{r}^{(t)T} \\mathbf{r}^{(t)}}$\n", " 5. $\\mathbf{p}^{(t+1)} \\gets - \\mathbf{r}^{(t+1)} + \\beta^{(t+1)} \\mathbf{p}^{(t)}$\n", " 6. $t \\gets t+1$\n", "\n", "* Computation cost per iteration is **one** matrix vector multiplication: $\\mathbf{A} \\mathbf{p}^{(t)}$. \n", "Consider PageRank problem, $\\mathbf{A}$ has dimension $n \\approx 10^{10}$ but is highly structured (sparse + low rank). Each matrix vector multiplication takes $O(n)$.\n", " \n", "* Theorem: If $\\mathbf{A}$ has $r$ distinct eigenvalues, $\\mathbf{x}^{(t)}$ converges to solution $\\mathbf{x}^*$ in at most $r$ steps." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pre-conditioned conjugate gradient (PCG)\n", "\n", "* Summary of conjugate gradient method for solving $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$ or equivalently minimizing $\\frac 12 \\mathbf{x}^T \\mathbf{A} \\mathbf{x} - \\mathbf{b}^T \\mathbf{x}$:\n", " * Each iteration needs one matrix vector multiplication: $\\mathbf{A} \\mathbf{p}^{(t+1)}$. For structured $\\mathbf{A}$, often $O(n)$ cost per iteration.\n", " * Guaranteed to converge in $n$ steps.\n", " \n", "* Two important bounds for conjugate gradient algorithm:\n", "\n", " Let $\\lambda_1 \\le \\cdots \\le \\lambda_n$ be the ordered eigenvalues of a pd $\\mathbf{A}$. \n", "$$\n", "\\begin{eqnarray*}\n", " \\|\\mathbf{x}^{(t+1)} - \\mathbf{x}^*\\|_{\\mathbf{A}}^2 &\\le& \\left( \\frac{\\lambda_{n-t} - \\lambda_1}{\\lambda_{n-t} + \\lambda_1} \\right)^2 \\|\\mathbf{x}^{(0)} - \\mathbf{x}^*\\|_{\\mathbf{A}}^2 \\\\\n", " \\|\\mathbf{x}^{(t+1)} - \\mathbf{x}^*\\|_{\\mathbf{A}}^2 &\\le& 2 \\left( \\frac{\\sqrt{\\kappa(\\mathbf{A})}-1}{\\sqrt{\\kappa(\\mathbf{A})}+1} \\right)^{t} \\|\\mathbf{x}^{(0)} - \\mathbf{x}^*\\|_{\\mathbf{A}}^2,\n", "\\end{eqnarray*}\n", "$$\n", "where $\\kappa(\\mathbf{A}) = \\lambda_n/\\lambda_1$ is the condition number of $\\mathbf{A}$.\n", "\n", "\n", "\n", "\n", "\n", "* Messages:\n", " * Roughly speaking, if the eigenvalues of $\\mathbf{A}$ occur in $r$ distinct clusters, the CG iterates will _approximately_ solve the problem after $O(r)$ steps. \n", " * $\\mathbf{A}$ with a small condition number ($\\lambda_1 \\approx \\lambda_n$) converges fast.\n", " \n", "* **Pre-conditioning**: Change of variables $\\widehat{\\mathbf{x}} = \\mathbf{C} \\mathbf{x}$ via a nonsingular $\\mathbf{C}$ and solve\n", "$$\n", "\t(\\mathbf{C}^{-T} \\mathbf{A} \\mathbf{C}^{-1}) \\widehat{\\mathbf{x}} = \\mathbf{C}^{-T} \\mathbf{b}.\n", "$$\n", "Choose $\\mathbf{C}$ such that \n", " * $\\mathbf{C}^{-T} \\mathbf{A} \\mathbf{C}^{-1}$ has small condition number, or \n", " * $\\mathbf{C}^{-T} \\mathbf{A} \\mathbf{C}^{-1}$ has clustered eigenvalues\n", " * Inexpensive solution of $\\mathbf{C}^T \\mathbf{C} \\mathbf{y} = \\mathbf{r}$\n", " \n", "* Preconditioned CG does not make use of $\\mathbf{C}$ explicitly, but rather the matrix $\\mathbf{M} = \\mathbf{C}^T \\mathbf{C}$.\n", "\n", "\n", "* **Preconditioned CG (PCG)** algorithm: \n", "\n", " 0. Given $\\mathbf{x}^{(0)}$, pre-conditioner $\\mathbf{M}$\n", " 0. $\\mathbf{r}^{(0)} \\gets \\mathbf{A} \\mathbf{x}^{(0)} - \\mathbf{b}$\n", " 0. solve $\\mathbf{M} \\mathbf{y}^{(0)} = \\mathbf{r}^{(0)}$ for $\\mathbf{y}^{(0)}$\n", " 0. $\\mathbf{p}^{(0)} \\gets - \\mathbf{r}^{(0)}$, $t=0$\n", " 0. While $\\mathbf{r}^{(0)} \\ne \\mathbf{0}$\n", " 1. $\\alpha^{(t)} \\gets \\frac{\\mathbf{r}^{(t)T} \\mathbf{y}^{(t)}}{\\mathbf{p}^{(t)T} \\mathbf{A} \\mathbf{p}^{(t)}}$\n", " 2. $\\mathbf{x}^{(t+1)} \\gets \\mathbf{x}^{(t)} + \\alpha^{(t)} \\mathbf{p}^{(t)}$\n", " 3. $\\mathbf{r}^{(t+1)} \\gets \\mathbf{r}^{(t)} + \\alpha^{(t)} \\mathbf{A} \\mathbf{p}^{(t)}$\n", " 4. Solve $\\mathbf{M} \\mathbf{y}^{(t+1)} \\gets \\mathbf{r}^{(t+1)}$ for $\\mathbf{y}^{(t+1)}$\n", " 5. $\\beta^{(t+1)} \\gets \\frac{\\mathbf{r}^{(t+1)T} \\mathbf{y}^{(t+1)}}{\\mathbf{r}^{(t)T} \\mathbf{r}^{(t)}}$\n", " 6. $\\mathbf{p}^{(t+1)} \\gets - \\mathbf{y}^{(t+1)} + \\beta^{(t+1)} \\mathbf{p}^{(t)}$\n", " 7. $t \\gets t+1$\n", "\n", " Remark: Only extra cost in the pre-conditioned CG algorithm is the need to solve the linear system $\\mathbf{M} \\mathbf{y} = \\mathbf{r}$.\n", " \n", "* Pre-conditioning is more like an art than science. Some choices include \n", " * Incomplete Cholesky. $\\mathbf{A} \\approx \\tilde{\\mathbf{L}} \\tilde{\\mathbf{L}}^T$, where $\\tilde{\\mathbf{L}}$ is a sparse approximate Cholesky factor. Then $\\tilde{\\mathbf{L}}^{-1} \\mathbf{A} \\tilde{\\mathbf{L}}^{-T} \\approx \\mathbf{I}$ (perfectly conditioned) and $\\mathbf{M} \\mathbf{y} = \\tilde{\\mathbf{L}} \\tilde {\\mathbf{L}}^T \\mathbf{y} = \\mathbf{r}$ is easy to solve. \n", " * Banded pre-conditioners. \n", " * Choose $\\mathbf{M}$ as a coarsened version of $\\mathbf{A}$.\n", " * Subject knowledge. Knowledge about the structure and origin of a problem is often the key to devising efficient pre-conditioner. For example, see recent work of Stein, Chen, Anitescu (2012) for pre-conditioning large covariance matrices. http://epubs.siam.org/doi/abs/10.1137/110834469" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example of PCG\n", "\n", "[Preconditioners.jl](https://github.com/mohamed82008/Preconditioners.jl) wraps a bunch of preconditioners.\n", "\n", "We use the Wathen matrix (sparse and positive definite) as a test matrix." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "include group.jl for user defined matrix generators\n", "verify download of index files...\n", "used remote site is https://sparse.tamu.edu/?per_page=All\n", "populating internal database...\n" ] }, { "data": { "text/plain": [ "30401×30401 SparseMatrixCSC{Float64,Int64} with 471601 stored entries:\n", " [1 , 1] = 2.18448\n", " [2 , 1] = -2.18448\n", " [3 , 1] = 0.728159\n", " [202 , 1] = -2.18448\n", " [203 , 1] = -2.91263\n", " [303 , 1] = 0.728159\n", " [304 , 1] = -2.91263\n", " [305 , 1] = 1.09224\n", " [1 , 2] = -2.18448\n", " [2 , 2] = 11.6505\n", " [3 , 2] = -2.18448\n", " [202 , 2] = 7.28159\n", " ⋮\n", " [30200, 30400] = 15.0471\n", " [30399, 30400] = -4.51414\n", " [30400, 30400] = 24.0754\n", " [30401, 30400] = -4.51414\n", " [30097, 30401] = 2.25707\n", " [30098, 30401] = -6.01885\n", " [30099, 30401] = 1.50471\n", " [30199, 30401] = -6.01885\n", " [30200, 30401] = -4.51414\n", " [30399, 30401] = 1.50471\n", " [30400, 30401] = -4.51414\n", " [30401, 30401] = 4.51414" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using BenchmarkTools, MatrixDepot, IterativeSolvers, LinearAlgebra, SparseArrays\n", "\n", "# Wathen matrix of dimension 30401 x 30401\n", "A = matrixdepot(\"wathen\", 100)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[1m Sparsity Pattern\u001b[22m\n", "\u001b[90m ┌──────────────────────────────────────────┐\u001b[39m \n", " \u001b[90m1\u001b[39m\u001b[90m │\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[31m> 0\u001b[39m\n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[34m< 0\u001b[39m\n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[35m⡀\u001b[39m\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m30401\u001b[39m\u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[35m⠈\u001b[39m\u001b[35m⠻\u001b[39m\u001b[35m⣦\u001b[39m\u001b[90m│\u001b[39m \n", "\u001b[90m └──────────────────────────────────────────┘\u001b[39m \n", "\u001b[90m 1\u001b[39m\u001b[90m \u001b[39m\u001b[90m 30401\u001b[39m\n", "\u001b[0m nz = 471601" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using UnicodePlots\n", "spy(A)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0005102687577359558" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# sparsity level\n", "count(!iszero, A) / length(A)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 1.16 MiB\n", " allocs estimate: 23\n", " --------------\n", " minimum time: 195.939 ms (0.00% GC)\n", " median time: 199.963 ms (0.00% GC)\n", " mean time: 200.412 ms (0.21% GC)\n", " maximum time: 210.798 ms (5.05% GC)\n", " --------------\n", " samples: 25\n", " evals/sample: 1" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# rhs\n", "b = ones(size(A, 1))\n", "# solve Ax=b by CG\n", "xcg = cg(A, b);\n", "@benchmark cg($A, $b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the incomplete cholesky preconditioner:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 22.336321 seconds (1.95 M allocations: 7.306 GiB, 1.01% gc time)\n" ] }, { "data": { "text/plain": [ "CholeskyPreconditioner{Float64,SparseMatrixCSC{Float64,Int64}}([1.478 0.0 … 0.0 0.0; -1.478 3.0767 … 0.0 0.0; … ; 0.0 0.0 … 3.48688 0.0; 0.0 0.0 … -0.209846 1.73063], 2)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Preconditioners\n", "@time p = CholeskyPreconditioner(A, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pre-conditioned conjugate gradient:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.2404710199864254e-7" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solver Ax=b by PCG\n", "xpcg = cg(A, b, Pl=p)\n", "# same answer?\n", "norm(xcg - xpcg)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 1.16 MiB\n", " allocs estimate: 32\n", " --------------\n", " minimum time: 13.705 ms (0.00% GC)\n", " median time: 14.927 ms (0.00% GC)\n", " mean time: 15.060 ms (0.71% GC)\n", " maximum time: 18.903 ms (10.66% GC)\n", " --------------\n", " samples: 332\n", " evals/sample: 1" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# PCG is >10 fold faster than CG\n", "@benchmark cg($A, $b, Pl=$p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other Krylov subspace methods\n", "\n", "* We leant about CG/PCG, which is for solving $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, $\\mathbf{A}$ pd.\n", "\n", "* **MINRES (minimum residual method)**: symmetric indefinite $\\mathbf{A}$.\n", "\n", "* **Bi-CG (bi-conjugate gradient)**: unsymmetric $\\mathbf{A}$.\n", "\n", "* **Bi-CGSTAB (Bi-CG stabilized)**: improved version of Bi-CG.\n", "\n", "* **GMRES (generalized minimum residual method)**: current _de facto_ method for unsymmetric $\\mathbf{A}$. E.g., PageRank problem.\n", "\n", "* **Lanczos method**: top eigen-pairs of a large symmetric matrix.\n", "\n", "* **Arnoldi method**: top eigen-pairs of a large unsymmetric matrix.\n", "\n", "* **Lanczos bidiagonalization** algorithm: top singular triplets of large matrix.\n", "\n", "* **LSQR**: least square problem $\\min \\|\\mathbf{y} - \\mathbf{X} \\beta\\|_2^2$. Algebraically equivalent to applying CG to the normal equation $(\\mathbf{X}^T \\mathbf{X} + \\lambda^2 I) \\beta = \\mathbf{X}^T \\mathbf{y}$.\n", "\n", "* **LSMR**: least square problem $\\min \\|\\mathbf{y} - \\mathbf{X} \\beta\\|_2^2$. Algebraically equivalent to applying MINRES to the normal equation $(\\mathbf{X}^T \\mathbf{X} + \\lambda^2 I) \\beta = \\mathbf{X}^T \\mathbf{y}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Software\n", "\n", "### Matlab \n", "\n", "* Iterative methods for solving linear equations: \n", " `pcg`, `bicg`, `bicgstab`, `gmres`, ...\n", "* Iterative methods for top eigen-pairs and singular pairs: \n", " `eigs`, `svds`, ...\n", "* Pre-conditioner: \n", " `cholinc`, `luinc`, ...\n", " \n", "* Get familiar with the **reverse communication interface (RCI)** for utilizing iterative solvers:\n", "```matlab\n", "x = gmres(A, b)\n", "x = gmres(@Afun, b)\n", "eigs(A)\n", "eigs(@Afun)\n", "```\n", " \n", "### Julia\n", "\n", "* `eigs` and `svds` in the [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl) package. [Numerical examples](http://hua-zhou.github.io/teaching/biostatm280-2019spring/slides/17-eigsvd/eigsvd.html#Lanczos/Arnoldi-iterative-method-for-top-eigen-pairs) later.\n", "\n", "* [`IterativeSolvers.jl`](https://github.com/JuliaMath/IterativeSolvers.jl) package. [CG numerical examples](http://hua-zhou.github.io/teaching/biostatm280-2019spring/slides/15-iterative/iterative.html#Numerical-examples)\n", "\n", "* Least squares examples:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Least squares example" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 1.42 GiB\n", " allocs estimate: 37941\n", " --------------\n", " minimum time: 5.171 s (1.41% GC)\n", " median time: 5.171 s (1.41% GC)\n", " mean time: 5.171 s (1.41% GC)\n", " maximum time: 5.171 s (1.41% GC)\n", " --------------\n", " samples: 1\n", " evals/sample: 1" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using BenchmarkTools, IterativeSolvers, LinearAlgebra, Random, SparseArrays\n", "\n", "Random.seed!(280) # seed\n", "n, p = 10000, 5000\n", "X = sprandn(n, p, 0.001) # iid standard normals with sparsity 0.01\n", "β = ones(p)\n", "y = X * β + randn(n)\n", "\n", "β̂_qr = X \\ y\n", "# least squares by QR\n", "@benchmark $X \\ $y" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(β̂_qr - β̂_lsqr) = 7.186597401024407e-5\n" ] }, { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 11.59 MiB\n", " allocs estimate: 1346\n", " --------------\n", " minimum time: 27.416 ms (0.00% GC)\n", " median time: 34.065 ms (11.22% GC)\n", " mean time: 34.335 ms (6.47% GC)\n", " maximum time: 45.718 ms (13.40% GC)\n", " --------------\n", " samples: 146\n", " evals/sample: 1" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "β̂_lsqr = lsqr(X, y)\n", "@show norm(β̂_qr - β̂_lsqr)\n", "# least squares by lsqr\n", "@benchmark lsqr($X, $y)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(β̂_qr - β̂_lsmr) = 0.02228791821986182\n" ] }, { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 4.06 MiB\n", " allocs estimate: 611\n", " --------------\n", " minimum time: 17.392 ms (0.00% GC)\n", " median time: 18.636 ms (0.00% GC)\n", " mean time: 19.384 ms (3.91% GC)\n", " maximum time: 26.789 ms (17.91% GC)\n", " --------------\n", " samples: 258\n", " evals/sample: 1" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "β̂_lsmr = lsmr(X, y)\n", "@show norm(β̂_qr - β̂_lsmr)\n", "# least squares by lsmr\n", "@benchmark lsmr($X, $y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Use LinearMaps in iterative solvers\n", "\n", "In many applications, it is advantageous to define linear maps indead of forming the actual (sparse) matrix. For a linear map, we need to specify how it acts on right- and left-multiplication on a vector. The [`LinearMaps.jl`](https://github.com/Jutho/LinearMaps.jl) package is exactly for this purpose and interfaces nicely with `IterativeSolvers.jl`, `Arnoldi.jl` and other iterative solver packages.\n", "\n", "Applications: \n", "1. The matrix is not sparse but admits special structure, e.g., easy + low rank (PageRank), Kronecker proudcts, etc. \n", "2. Less memory usage. \n", "3. Linear algebra on a standardized (centered and scaled) sparse matrix.\n", "\n", "Consider the differencing operator that takes differences between neighboring pixels." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearMaps.FunctionMap{Float64}(leftdiff!, mrightdiff!, 100, 100; ismutating=true, issymmetric=false, ishermitian=false, isposdef=false)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearMaps, IterativeSolvers\n", "\n", "# Overwrite y with A * x\n", "# left difference assuming periodic boundary conditions\n", "function leftdiff!(y::AbstractVector, x::AbstractVector) \n", " N = length(x)\n", " length(y) == N || throw(DimensionMismatch())\n", " @inbounds for i in 1:N\n", " y[i] = x[i] - x[mod1(i - 1, N)]\n", " end\n", " return y\n", "end\n", "\n", "# Overwrite y with A' * x\n", "# minus right difference\n", "function mrightdiff!(y::AbstractVector, x::AbstractVector) \n", " N = length(x)\n", " length(y) == N || throw(DimensionMismatch())\n", " @inbounds for i in 1:N\n", " y[i] = x[i] - x[mod1(i + 1, N)]\n", " end\n", " return y\n", "end\n", "\n", "# define linear map\n", "D = LinearMap{Float64}(leftdiff!, mrightdiff!, 100; ismutating=true) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Linear maps can be used like a regular matrix." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "size(D) = (100, 100)\n", "D * v = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n", "D' * v = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n" ] } ], "source": [ "@show size(D)\n", "v = ones(size(D, 2)) # vector of all 1s\n", "@show D * v\n", "@show D' * v;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we form the corresponding dense matrix, it will look like" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100×100 Array{Float64,2}:\n", " 1.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 -1.0\n", " -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 -1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 -1.0 1.0 … 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " ⋮ ⋮ ⋱ ⋮ \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 … -1.0 1.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 1.0" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Matrix(D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we form the corresponding sparse matrix, it will look like" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100×100 SparseMatrixCSC{Float64,Int64} with 200 stored entries:\n", " [1 , 1] = 1.0\n", " [2 , 1] = -1.0\n", " [2 , 2] = 1.0\n", " [3 , 2] = -1.0\n", " [3 , 3] = 1.0\n", " [4 , 3] = -1.0\n", " [4 , 4] = 1.0\n", " [5 , 4] = -1.0\n", " [5 , 5] = 1.0\n", " [6 , 5] = -1.0\n", " [6 , 6] = 1.0\n", " [7 , 6] = -1.0\n", " ⋮\n", " [95 , 95] = 1.0\n", " [96 , 95] = -1.0\n", " [96 , 96] = 1.0\n", " [97 , 96] = -1.0\n", " [97 , 97] = 1.0\n", " [98 , 97] = -1.0\n", " [98 , 98] = 1.0\n", " [99 , 98] = -1.0\n", " [99 , 99] = 1.0\n", " [100, 99] = -1.0\n", " [1 , 100] = -1.0\n", " [100, 100] = 1.0" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using SparseArrays\n", "sparse(D)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[1m Sparsity Pattern\u001b[22m\n", "\u001b[90m ┌──────────────────────────────────────────┐\u001b[39m \n", " \u001b[90m1\u001b[39m\u001b[90m │\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[90m│\u001b[39m \u001b[31m> 0\u001b[39m\n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \u001b[34m< 0\u001b[39m\n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[90m│\u001b[39m \n", " \u001b[90m100\u001b[39m\u001b[90m │\u001b[39m\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[0m⠀\u001b[34m⠈\u001b[39m\u001b[35m⠳\u001b[39m\u001b[35m⣄\u001b[39m\u001b[90m│\u001b[39m \n", "\u001b[90m └──────────────────────────────────────────┘\u001b[39m \n", "\u001b[90m 1\u001b[39m\u001b[90m \u001b[39m\u001b[90m 100\u001b[39m\n", "\u001b[0m nz = 200" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using UnicodePlots\n", "spy(sparse(D))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute top singular values using iterative method (Arnoldi)." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(SVD{Float64,Float64,Array{Float64,2}}([-0.1 0.0640197 -0.126101; 0.1 -0.0559754 0.129872; … ; -0.1 0.0793195 -0.117083; 0.1 -0.0718113 0.121832], [2.0, 1.99901, 1.99901], [0.1 -0.1 … 0.1 -0.1; -0.0600272 0.0518684 … -0.0756027 0.067949; 0.12805 -0.131566 … 0.119517 -0.124028]), 3, 35, 573, [0.0400615, 0.00610398, 0.0278884, 0.0776812, 0.08746, 0.00602873, 0.212726, 0.0333012, -0.0523244, 0.144026 … 0.083757, 0.0931424, -0.0237672, 0.00677622, 0.042755, 0.0557526, 0.0887591, 0.102905, 0.0669487, 0.0562545])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Arpack\n", "Arpack.svds(D, nsv=3)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100-element Array{Float64,1}:\n", " 2.0 \n", " 1.999013120731463 \n", " 1.9990131207314623 \n", " 1.9960534568565431 \n", " 1.9960534568565427 \n", " 1.99112392920616 \n", " 1.9911239292061598 \n", " 1.9842294026289555 \n", " 1.9842294026289549 \n", " 1.9753766811902758 \n", " 1.9753766811902744 \n", " 1.9645745014573777 \n", " 1.9645745014573768 \n", " ⋮ \n", " 0.37476262917144926 \n", " 0.3128689300804618 \n", " 0.31286893008046174 \n", " 0.25066646712860846 \n", " 0.25066646712860846 \n", " 0.18821662663702868 \n", " 0.18821662663702865 \n", " 0.1255810390586268 \n", " 0.12558103905862675 \n", " 0.06282151815625672 \n", " 0.06282151815625658 \n", " 2.5308506128625597e-18" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra\n", "# check solution against the direct method for SVD\n", "svdvals(Matrix(D))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute top eigenvalues of the Gram matrix `D'D` using iterative method (Arnoldi)." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([4.0, 3.99605, 3.99605], [0.1 -0.00939969 -0.141109; -0.1 0.000520858 0.14142; … ; 0.1 -0.0270112 -0.138818; -0.1 0.0182414 0.14024], 3, 30, 493, [0.103736, 0.0100948, 0.0228213, -0.0164503, 0.10995, 0.0865899, 0.124362, -0.161938, -0.0933529, 0.125877 … -0.0628246, -0.0637099, 0.0916849, 0.179155, 0.14872, 0.214858, 0.0147322, 0.146649, -0.206126, 0.0174464])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Arpack.eigs(D'D, nev=3, which=:LM)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further reading\n", "\n", "* Chapter 5 of [Numerical Optimization](https://ucla.worldcat.org/title/numerical-optimization/oclc/209918411&referer=brief_results) by Jorge Nocedal and Stephen Wright (1999).\n", "\n", "* Sections 11.3-11.5 of [Matrix Computations](https://ucla.worldcat.org/title/matrix-computations/oclc/824733531&referer=brief_results) by Gene Golub and Charles Van Loan (2013)." ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "135px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "skip_h1_title": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }