{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  Conjugate Gradient and Krylov Subspace Methods
1.1  Introduction
1.2  Conjugate gradient (CG) method
1.3  Pre-conditioned conjugate gradient (PCG)
1.3.1  Example of PCG
1.4  Other Krylov subspace methods
1.5  Software
1.5.1  Matlab
1.5.2  Julia
1.5.2.1  Least squares example
1.5.2.2  Use LinearMaps in iterative solvers
1.6  Further reading
" ] }, { "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 }