RandomizedSVD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomized Singular Value Decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tsvd (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "# Halko et. al.\n",
    "function tsvd(A::AbstractMatrix, k::Integer, numiter::Integer, numoversample::Integer)\n",
    "    m, n = size(A)\n",
    "    Q = rand(n, k + numoversample)\n",
    "    Q = A * Q\n",
    "    for i = 1:numiter\n",
    "        Q = transpose(transpose(Q) * A)\n",
    "        Q = Matrix(qr(Q).Q)\n",
    "        Q = A * Q\n",
    "        Q = Matrix(qr(Q).Q)\n",
    "    end\n",
    "    QA = transpose(Q) * A\n",
    "    Va, s, R = svd(transpose(QA), full=false)\n",
    "    U = Q * transpose(R)\n",
    "    U[:, 1:k], s[1:k], transpose(Va)[1:k, :]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.662047768214454e-16"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate a rank 1 matrix\n",
    "A = rand(2 ^ 12) * rand(2 ^ 10)'\n",
    "\n",
    "# Perform truncated SVD\n",
    "U, S, Vt = tsvd(A, 20, 4, 16)\n",
    "\n",
    "# Calculate error\n",
    "norm(A .- U * Diagonal(S) * Vt) / norm(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.3717996541638573e-15"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Do the same using builtin svd routine\n",
    "U′, S′, Vt′ = svd(A)\n",
    "norm(A .- U′ * Diagonal(S′) * transpose(Vt′)) / norm(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the extreme case of a rank-1 matrix, the truncated svd has great accuracy.\n",
    "Due the relatively high default tolerance settings on the bulitin svd routine, it even outperforms the builtin here.\n",
    "\n",
    "In terms of runtime performance, the builtin svd is no match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  26.95 MiB\n",
       "  allocs estimate:  133\n",
       "  --------------\n",
       "  minimum time:     30.781 ms (0.00% GC)\n",
       "  median time:      33.515 ms (5.03% GC)\n",
       "  mean time:        34.796 ms (3.13% GC)\n",
       "  maximum time:     52.626 ms (4.78% GC)\n",
       "  --------------\n",
       "  samples:          144\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using BenchmarkTools\n",
    "\n",
    "@benchmark tsvd(A, 20, 4, 16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  104.13 MiB\n",
       "  allocs estimate:  13\n",
       "  --------------\n",
       "  minimum time:     238.864 ms (0.13% GC)\n",
       "  median time:      245.052 ms (1.94% GC)\n",
       "  mean time:        246.937 ms (1.83% GC)\n",
       "  maximum time:     268.918 ms (1.73% GC)\n",
       "  --------------\n",
       "  samples:          21\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark svd(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Block-Krylov SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "bksvd (generic function with 3 methods)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Musco et. al.\n",
    "function bksvd(A, bsize=20, iter=4)\n",
    "    l = ones(size(A, 1), 1)\n",
    "    u = zeros(1, size(A, 2))\n",
    "    \n",
    "    K = zeros(size(A, 2), bsize * iter)\n",
    "    # Random block initialization\n",
    "    block = randn(size(A, 2), bsize)\n",
    "    block = Matrix(qr(block).Q)\n",
    "    # Preallocate space for temporary products\n",
    "    T = zeros(size(A, 2), bsize)\n",
    "\n",
    "    # Construct and orthonormalize Krylov subspace\n",
    "    # Orthogonalize at each step using economy size QR decomposition\n",
    "    for i = 1:iter\n",
    "        T = A * block - l * (u * block)\n",
    "        block = A' * T - u' * (l' * T)\n",
    "        block = Matrix(qr(block).Q)\n",
    "        K[:, ((i - 1) * bsize + 1):(i * bsize)] = block\n",
    "    end\n",
    "    Q = Matrix(qr(K).Q)\n",
    "    \n",
    "    # Rayleigh-Ritz postprocessing with economy size dense SVD\n",
    "    T = A * Q - l * (u * Q)\n",
    "\n",
    "    Ut, St, Vt = svd(T)\n",
    "    U = Ut[:, 1:bsize]\n",
    "    S = St[1:bsize]\n",
    "    V = Q * Vt'\n",
    "    U, S, V[:, 1:bsize]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.050727967517019e-16"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U′, S′, V′ = bksvd(A, 20)\n",
    "norm(A .- U′ * Diagonal(S′) * transpose(V′)) / norm(A)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.2",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}