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

1  A list of \"easy\" linear systems
1.1  Diagonal matrix
1.2  Bidiagonal, tridiagonal, and banded matrices
1.3  Triangular matrix
1.4  Block diagonal matrix
1.5  Kronecker product
1.6  Sparse matrix
1.6.1  Matrix-vector multiplication
1.6.2  Solve linear equation
1.7  Easy plus low rank
1.7.1  Solve linear equation
1.7.2  Matrix-vector multiplication
1.7.3  Determinant
1.8  Easy plus border
1.9  Orthogonal matrix
1.10  Toeplitz matrix
1.11  Circulant matrix
1.12  Vandermonde matrix
1.13  Cauchy-like matrix
1.14  Structured-rank matrix
" ] }, { "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": [ "# A list of \"easy\" linear systems\n", "\n", "Consider $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$. Or, consider matrix inverse (if you want). $\\mathbf{A}$ can be huge. Keep massive data in mind: 1000 Genome Project, NetFlix, Google PageRank, finance, spatial statistics, ... We should be alert to many easy linear systems. \n", "\n", "Don't blindly use `A \\ b` and `inv` in Julia or `solve` function in R. **Don't waste computing resources by bad choices of algorithms!**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Diagonal matrix\n", "\n", "Diagonal $\\mathbf{A}$: $n$ flops. Use `Diagonal` type of Julia." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "using BenchmarkTools, LinearAlgebra, Random\n", "\n", "# generate random data\n", "Random.seed!(280)\n", "n = 1000\n", "A = diagm(0 => randn(n)) # a diagonal matrix stored as Matrix{Float64}\n", "b = randn(n);" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in LinearAlgebra at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:892" ], "text/plain": [ "\\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in LinearAlgebra at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:892" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# should give link: https://github.com/JuliaLang/julia/blob/5b637df34396034b0dd353e603ab3d61322369fb/stdlib/LinearAlgebra/src/generic.jl#L956\n", "@which A \\ b" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 15.89 KiB\n", " allocs estimate: 3\n", " --------------\n", " minimum time: 761.706 μs (0.00% GC)\n", " median time: 833.682 μs (0.00% GC)\n", " mean time: 872.122 μs (0.22% GC)\n", " maximum time: 9.602 ms (91.45% GC)\n", " --------------\n", " samples: 5665\n", " evals/sample: 1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check `istril(A)` and `istriu(A)` (O(n^2)), then call `Diagonal(A) \\ b` (O(n))\n", "@benchmark $A \\ $b" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 15.89 KiB\n", " allocs estimate: 3\n", " --------------\n", " minimum time: 2.752 μs (0.00% GC)\n", " median time: 3.731 μs (0.00% GC)\n", " mean time: 5.059 μs (28.25% GC)\n", " maximum time: 4.938 ms (99.86% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 9" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# O(n) computation, no extra array allocation\n", "@benchmark Diagonal($A) \\ $b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bidiagonal, tridiagonal, and banded matrices\n", "\n", "Bidiagonal, tridiagonal, or banded $\\mathbf{A}$: Band LU, band Cholesky, ... roughly $O(n)$ flops. \n", "* Use [`Bidiagonal`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.Bidiagonal), [`Tridiagonal`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.Tridiagonal), [`SymTridiagonal`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.SymTridiagonal) types of Julia." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1000×1000 SymTridiagonal{Float64,Array{Float64,1}}:\n", " 0.126238 0.618244 ⋅ … ⋅ ⋅ ⋅ \n", " 0.618244 -2.34688 1.10206 ⋅ ⋅ ⋅ \n", " ⋅ 1.10206 1.91661 ⋅ ⋅ ⋅ \n", " ⋅ ⋅ -0.447244 ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋮ ⋱ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ \n", " ⋅ ⋅ ⋅ 0.709662 ⋅ ⋅ \n", " ⋅ ⋅ ⋅ 0.542815 -0.363838 ⋅ \n", " ⋅ ⋅ ⋅ -0.363838 -1.04034 0.321148\n", " ⋅ ⋅ ⋅ ⋅ 0.321148 -0.276537" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Random.seed!(280) \n", "\n", "n = 1000\n", "dv = randn(n)\n", "ev = randn(n - 1)\n", "b = randn(n) # rhs\n", "# symmetric tridiagonal matrix\n", "A = SymTridiagonal(dv, ev)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.65 MiB\n", " allocs estimate: 5\n", " --------------\n", " minimum time: 9.502 ms (0.00% GC)\n", " median time: 10.496 ms (0.00% GC)\n", " mean time: 10.635 ms (4.82% GC)\n", " maximum time: 47.393 ms (80.36% GC)\n", " --------------\n", " samples: 470\n", " evals/sample: 1" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# convert to a full matrix\n", "Afull = Matrix(A)\n", "\n", "# LU decomposition (2/3) n^3 flops!\n", "@benchmark $Afull \\ $b" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 23.86 KiB\n", " allocs estimate: 5\n", " --------------\n", " minimum time: 13.397 μs (0.00% GC)\n", " median time: 15.377 μs (0.00% GC)\n", " mean time: 22.577 μs (27.96% GC)\n", " maximum time: 42.611 ms (99.94% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# specialized algorithm for tridiagonal matrix, O(n) flops\n", "@benchmark $A \\ $b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Triangular matrix\n", "\n", "Triangular $\\mathbf{A}$: $n^2$ flops to solve linear system." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.94 KiB\n", " allocs estimate: 1\n", " --------------\n", " minimum time: 707.663 μs (0.00% GC)\n", " median time: 768.855 μs (0.00% GC)\n", " mean time: 775.331 μs (0.00% GC)\n", " maximum time: 1.453 ms (0.00% GC)\n", " --------------\n", " samples: 6362\n", " evals/sample: 1" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Random.seed!(280)\n", "\n", "n = 1000\n", "A = tril(randn(n, n)) # a lower-triangular matrix stored as Matrix{Float64}\n", "b = randn(n)\n", "\n", "# check istril() then triangular solve\n", "@benchmark $A \\ $b" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.94 KiB\n", " allocs estimate: 1\n", " --------------\n", " minimum time: 313.525 μs (0.00% GC)\n", " median time: 346.465 μs (0.00% GC)\n", " mean time: 375.613 μs (0.00% GC)\n", " maximum time: 1.111 ms (0.00% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# triangular solve directly; save the cost of istril()\n", "@benchmark LowerTriangular($A) \\ $b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Block diagonal matrix\n", "\n", "Block diagonal: Suppose $n = \\sum_b n_b$. For linear equations, $(\\sum_b n_b)^3$ (without using block diagonal structure) vs $\\sum_b n_b^3$ (using block diagonal structure). \n", "\n", "Julia has a [`blockdiag`](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#SparseArrays.blockdiag) function that generates a **sparse** matrix. **Anyone interested writing a `BlockDiagonal.jl` package?**" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1000×1000 SparseMatrixCSC{Float64,Int64} with 969 stored entries:\n", " [31 , 1] = 2.07834\n", " [53 , 1] = -1.11883\n", " [58 , 1] = -0.66448\n", " [14 , 4] = 1.11793\n", " [96 , 5] = 1.22813\n", " [81 , 8] = -0.919643\n", " [48 , 9] = 1.0185\n", " [49 , 9] = -0.996332\n", " [15 , 10] = 1.30841\n", " [28 , 10] = -0.818757\n", " [39 , 11] = 1.08248\n", " [82 , 11] = -0.0102294\n", " ⋮\n", " [935 , 987] = 0.677319\n", " [956 , 987] = -0.900804\n", " [967 , 987] = -0.438788\n", " [971 , 991] = 0.176756\n", " [929 , 992] = -1.17384\n", " [974 , 993] = 1.59235\n", " [967 , 994] = 0.542169\n", " [994 , 995] = 0.627832\n", " [998 , 997] = 0.60382\n", " [935 , 998] = 0.342675\n", " [947 , 998] = 0.482228\n", " [975 , 1000] = 0.991598" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using SparseArrays\n", "\n", "Random.seed!(280)\n", "\n", "B = 10 # number of blocks\n", "ni = 100\n", "A = blockdiag([sprandn(ni, ni, 0.01) for b in 1:B]...)" ] }, { "cell_type": "code", "execution_count": 12, "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[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[31m> 0\u001b[39m\n", " \u001b[90m │\u001b[39m\u001b[35m⢽\u001b[39m\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[31m⠁\u001b[39m\u001b[0m⠀\u001b[31m⠉\u001b[39m\u001b[31m⠈\u001b[39m\u001b[35m⢦\u001b[39m\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[35m⣐\u001b[39m\u001b[31m⡂\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[31m⠈\u001b[39m\u001b[35m⠈\u001b[39m\u001b[35m⠉\u001b[39m\u001b[31m⠁\u001b[39m\u001b[35m⢷\u001b[39m\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[34m⢪\u001b[39m\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[34m⠈\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⠁\u001b[39m\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[35m⢘\u001b[39m\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[34m⠈\u001b[39m\u001b[35m⠑\u001b[39m\u001b[35m⠙\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⢊\u001b[39m\u001b[35m⡆\u001b[39m\u001b[35m⢦\u001b[39m\u001b[35m⣄\u001b[39m\u001b[34m⣲\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[35m⢐\u001b[39m\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[31m⠈\u001b[39m\u001b[35m⠉\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⠛\u001b[39m\u001b[35m⡤\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⣤\u001b[39m\u001b[35m⣤\u001b[39m\u001b[34m⡄\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[35m⣿\u001b[39m\u001b[35m⣿\u001b[39m\u001b[35m⡺\u001b[39m\u001b[35m⠊\u001b[39m\u001b[34m⠆\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[35m⠝\u001b[39m\u001b[35m⠪\u001b[39m\u001b[35m⠷\u001b[39m\u001b[35m⠚\u001b[39m\u001b[35m⣂\u001b[39m\u001b[35m⡤\u001b[39m\u001b[35m⣠\u001b[39m\u001b[35m⣠\u001b[39m\u001b[31m⡄\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[35m⣻\u001b[39m\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[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[34m⢃\u001b[39m\u001b[35m⣠\u001b[39m\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[35m⠸\u001b[39m\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[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[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[35m⢸\u001b[39m\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[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[35m⡧\u001b[39m\u001b[0m⠀\u001b[34m⡀\u001b[39m\u001b[34m⣀\u001b[39m\u001b[34m⡀\u001b[39m\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[31m⢠\u001b[39m\u001b[35m⣾\u001b[39m\u001b[35m⢽\u001b[39m\u001b[35m⡿\u001b[39m\u001b[35m⢷\u001b[39m\u001b[90m│\u001b[39m \n", " \u001b[90m1000\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[35m⠘\u001b[39m\u001b[35m⢻\u001b[39m\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 1000\u001b[39m\n", "\u001b[0m nz = 969" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using UnicodePlots\n", "spy(A)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Kronecker product\n", "\n", "Use\n", "$$\n", "\\begin{eqnarray*}\n", " (\\mathbf{A} \\otimes \\mathbf{B})^{-1} &=& \\mathbf{A}^{-1} \\otimes \\mathbf{B}^{-1} \\\\\n", " (\\mathbf{C}^T \\otimes \\mathbf{A}) \\text{vec}(\\mathbf{B}) &=& \\text{vec}(\\mathbf{A} \\mathbf{B} \\mathbf{C}) \\\\\n", " \\text{det}(\\mathbf{A} \\otimes \\mathbf{B}) &=& [\\text{det}(\\mathbf{A})]^p [\\text{det}(\\mathbf{B})]^m, \\quad \\mathbf{A} \\in \\mathbb{R}^{m \\times m}, \\mathbf{B} \\in \\mathbb{R}^{p \\times p}\n", "\\end{eqnarray*} \n", "$$\n", "to avoid forming and doing costly computation on the potentially huge Kronecker $\\mathbf{A} \\otimes \\mathbf{B}$.\n", "\n", "**Anyone interested writing a package?**" ] }, { "cell_type": "code", "execution_count": 13, "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": [ "50×50 Array{Float64,2}:\n", " 1.0 0.5 0.333333 0.25 … 0.0208333 0.0204082 0.02\n", " 0.5 1.0 0.666667 0.5 0.0416667 0.0408163 0.04\n", " 0.333333 0.666667 1.0 0.75 0.0625 0.0612245 0.06\n", " 0.25 0.5 0.75 1.0 0.0833333 0.0816327 0.08\n", " 0.2 0.4 0.6 0.8 0.104167 0.102041 0.1 \n", " 0.166667 0.333333 0.5 0.666667 … 0.125 0.122449 0.12\n", " 0.142857 0.285714 0.428571 0.571429 0.145833 0.142857 0.14\n", " 0.125 0.25 0.375 0.5 0.166667 0.163265 0.16\n", " 0.111111 0.222222 0.333333 0.444444 0.1875 0.183673 0.18\n", " 0.1 0.2 0.3 0.4 0.208333 0.204082 0.2 \n", " 0.0909091 0.181818 0.272727 0.363636 … 0.229167 0.22449 0.22\n", " 0.0833333 0.166667 0.25 0.333333 0.25 0.244898 0.24\n", " 0.0769231 0.153846 0.230769 0.307692 0.270833 0.265306 0.26\n", " ⋮ ⋱ \n", " 0.025641 0.0512821 0.0769231 0.102564 0.8125 0.795918 0.78\n", " 0.025 0.05 0.075 0.1 0.833333 0.816327 0.8 \n", " 0.0243902 0.0487805 0.0731707 0.097561 … 0.854167 0.836735 0.82\n", " 0.0238095 0.047619 0.0714286 0.0952381 0.875 0.857143 0.84\n", " 0.0232558 0.0465116 0.0697674 0.0930233 0.895833 0.877551 0.86\n", " 0.0227273 0.0454545 0.0681818 0.0909091 0.916667 0.897959 0.88\n", " 0.0222222 0.0444444 0.0666667 0.0888889 0.9375 0.918367 0.9 \n", " 0.0217391 0.0434783 0.0652174 0.0869565 … 0.958333 0.938776 0.92\n", " 0.0212766 0.0425532 0.0638298 0.0851064 0.979167 0.959184 0.94\n", " 0.0208333 0.0416667 0.0625 0.0833333 1.0 0.979592 0.96\n", " 0.0204082 0.0408163 0.0612245 0.0816327 0.979592 1.0 0.98\n", " 0.02 0.04 0.06 0.08 0.96 0.98 1.0 " ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using MatrixDepot, LinearAlgebra\n", "\n", "A = matrixdepot(\"lehmer\", 50) # a pd matrix" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100×100 Array{Float64,2}:\n", " 0.88678 0.117729 -0.0552891 … 1.52453e-18 -3.11225e-18\n", " 0.117729 0.732177 0.227559 -4.47253e-18 9.13044e-18\n", " -0.0552891 0.227559 0.401882 9.04011e-18 -1.84549e-17\n", " -0.000257948 -0.00903794 0.0464366 -1.07207e-17 2.18857e-17\n", " -0.00361503 0.0096735 -0.0214526 3.98827e-17 -8.14185e-17\n", " 0.00042184 -0.00124191 0.0024795 … -1.12045e-16 2.28734e-16\n", " -8.70147e-5 0.000242067 -0.00049947 2.82936e-16 -5.77599e-16\n", " 6.41788e-6 -9.33946e-5 0.000342935 -4.14338e-16 8.45851e-16\n", " -4.65439e-5 0.000107093 -0.000195109 4.39702e-16 -8.9763e-16 \n", " -2.08953e-6 7.84307e-7 1.28448e-5 -8.06945e-17 1.64734e-16\n", " 4.92703e-7 -1.17028e-6 1.39914e-6 … 1.90854e-16 -3.89625e-16\n", " -6.70515e-7 1.93176e-6 -3.73725e-6 -1.01747e-15 2.07716e-15\n", " 2.43687e-7 -7.52096e-7 1.61857e-6 7.99215e-16 -1.63158e-15\n", " ⋮ ⋱ \n", " 4.43802e-18 -1.30198e-17 2.63164e-17 -0.00372067 0.00622875 \n", " -2.63663e-18 7.73509e-18 -1.56346e-17 0.00322562 -0.00443629 \n", " 3.01875e-18 -8.85614e-18 1.79005e-17 … -0.00241141 0.00639837 \n", " -2.7251e-18 7.99463e-18 -1.61592e-17 0.00302995 -0.00631032 \n", " 3.23187e-18 -9.48136e-18 1.91642e-17 -0.0042048 0.00769745 \n", " -2.97621e-18 8.73132e-18 -1.76482e-17 0.0045267 -0.00739848 \n", " 2.79334e-18 -8.19483e-18 1.65638e-17 -0.00460457 0.00736033 \n", " -7.34735e-19 2.1555e-18 -4.3568e-18 … 0.00382355 -0.00836409 \n", " 6.96272e-19 -2.04266e-18 4.12873e-18 -0.0394007 0.0187525 \n", " -1.59219e-18 4.67103e-18 -9.44134e-18 0.304724 -0.0667661 \n", " 1.52453e-18 -4.47253e-18 9.04011e-18 0.645854 0.106467 \n", " -3.11225e-18 9.13044e-18 -1.84549e-17 0.106467 0.175967 " ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = matrixdepot(\"oscillate\", 100) # pd matrix" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "M = kron(A, B)\n", "c = ones(size(M, 2)) # rhs\n", "# Method 1: form Kronecker product and Cholesky solve\n", "x1 = cholesky(Symmetric(M)) \\ c;" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Method 2: use (A ⊗ B)^{-1} = A^{-1} ⊗ B^{-1}\n", "m, p = size(A, 1), size(B, 1)\n", "x2 = vec(transpose(cholesky(Symmetric(A)) \\ \n", " transpose(cholesky(Symmetric(B)) \\ reshape(c, p, m))));" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.358949337519373e-7" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# relative error\n", "norm(x1 - x2) / norm(x1)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 381.51 MiB\n", " allocs estimate: 11\n", " --------------\n", " minimum time: 630.446 ms (3.03% GC)\n", " median time: 711.411 ms (3.10% GC)\n", " mean time: 702.294 ms (5.36% GC)\n", " maximum time: 741.308 ms (2.74% GC)\n", " --------------\n", " samples: 8\n", " evals/sample: 1" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using BenchmarkTools\n", "\n", "# Method 1: form Kronecker and Cholesky solve\n", "@benchmark cholesky(Symmetric(kron(A, B))) \\ c" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 176.59 KiB\n", " allocs estimate: 24\n", " --------------\n", " minimum time: 131.506 μs (0.00% GC)\n", " median time: 236.375 μs (0.00% GC)\n", " mean time: 273.470 μs (9.59% GC)\n", " maximum time: 6.839 ms (69.43% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Method 2: use (A ⊗ B)^{-1} = A^{-1} ⊗ B^{-1}\n", "@benchmark vec(transpose(cholesky(Symmetric(A)) \\ \n", " transpose(cholesky(Symmetric(B)) \\ reshape(c, p, m))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sparse matrix\n", "\n", "Sparsity: sparse matrix decomposition or iterative method. \n", "* The easiest recognizable structure. Familiarize yourself with the sparse matrix computation tools in Julia, Matlab, R (`Matrix` package), MKL (sparse BLAS), ... as much as possible." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.001994776158751544" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using MatrixDepot\n", "\n", "Random.seed!(280)\n", "\n", "# a 7701-by-7701 sparse pd matrix\n", "A = matrixdepot(\"wathen\", 50)\n", "# random generated rhs\n", "b = randn(size(A, 1))\n", "Afull = Matrix(A)\n", "count(!iszero, A) / length(A) # sparsity" ] }, { "cell_type": "code", "execution_count": 21, "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[90m7701\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 7701\u001b[39m\n", "\u001b[0m nz = 118301" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using UnicodePlots\n", "spy(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrix-vector multiplication" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 60.27 KiB\n", " allocs estimate: 2\n", " --------------\n", " minimum time: 22.048 ms (0.00% GC)\n", " median time: 25.870 ms (0.00% GC)\n", " mean time: 25.786 ms (0.00% GC)\n", " maximum time: 29.732 ms (0.00% GC)\n", " --------------\n", " samples: 194\n", " evals/sample: 1" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# dense matrix-vector multiplication\n", "@benchmark $Afull * $b" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 60.27 KiB\n", " allocs estimate: 2\n", " --------------\n", " minimum time: 99.190 μs (0.00% GC)\n", " median time: 135.390 μs (0.00% GC)\n", " mean time: 150.205 μs (5.97% GC)\n", " maximum time: 4.153 ms (96.20% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# sparse matrix-vector multiplication\n", "@benchmark $A * $b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solve linear equation" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 452.52 MiB\n", " allocs estimate: 8\n", " --------------\n", " minimum time: 1.576 s (0.67% GC)\n", " median time: 1.631 s (1.28% GC)\n", " mean time: 1.624 s (1.85% GC)\n", " maximum time: 1.658 s (1.27% GC)\n", " --------------\n", " samples: 4\n", " evals/sample: 1" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve via dense Cholesky\n", "xchol = cholesky(Afull) \\ b\n", "@benchmark cholesky($Afull) \\ $b" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xchol - xcholsp) = 3.7385578057004605e-15\n" ] }, { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 13.45 MiB\n", " allocs estimate: 55\n", " --------------\n", " minimum time: 15.046 ms (0.00% GC)\n", " median time: 17.306 ms (4.78% GC)\n", " mean time: 17.089 ms (2.95% GC)\n", " maximum time: 20.649 ms (4.04% GC)\n", " --------------\n", " samples: 293\n", " evals/sample: 1" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve via sparse Cholesky\n", "xcholsp = cholesky(A) \\ b\n", "@show norm(xchol - xcholsp)\n", "@benchmark cholesky($A) \\ $b" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(xcg - xchol) = 2.1854385431265016e-7\n" ] }, { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 302.20 KiB\n", " allocs estimate: 23\n", " --------------\n", " minimum time: 29.459 ms (0.00% GC)\n", " median time: 30.978 ms (0.00% GC)\n", " mean time: 31.294 ms (0.12% GC)\n", " maximum time: 41.731 ms (0.00% GC)\n", " --------------\n", " samples: 160\n", " evals/sample: 1" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# sparse solve via conjugate gradient\n", "using IterativeSolvers\n", "\n", "xcg = cg(A, b)\n", "@show norm(xcg - xchol)\n", "@benchmark cg($A, $b)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Easy plus low rank\n", "\n", "Easy plus low rank: $\\mathbf{U} \\in \\mathbb{R}^{n \\times r}$, $\\mathbf{V} \\in \\mathbb{R}^{r \\times n}$, $r \\ll n$. Woodbury formula\n", "\\begin{eqnarray*}\n", "\t(\\mathbf{A} + \\mathbf{U} \\mathbf{V}^T)^{-1} &=& \\mathbf{A}^{-1} - \\mathbf{A}^{-1} \\mathbf{U} (\\mathbf{I}_r + \\mathbf{V}^T \\mathbf{A}^{-1} \\mathbf{U})^{-1} \\mathbf{V}^T \\mathbf{A}^{-1} \\\\\n", " \\text{det}(\\mathbf{A} + \\mathbf{U} \\mathbf{V}^T) &=& \\text{det}(\\mathbf{A}) \\text{det}(\\mathbf{I}_r + \\mathbf{V} \\mathbf{A}^{-1} \\mathbf{U}^T).\n", "\\end{eqnarray*}\n", "\n", "* Keep HW2 Q2 (multivariate density) and PageRank problem in mind. \n", "* [`WoodburyMatrices.jl`](https://github.com/timholy/WoodburyMatrices.jl) package can be useful." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1000×1000 Symmetric{Float64,Array{Float64,2}}:\n", " 1.8571 0.513107 0.872146 … 0.764278 -0.241331 0.54921 \n", " 0.513107 4.57505 -0.636972 -1.86465 -1.92237 -1.72569 \n", " 0.872146 -0.636972 4.81387 1.99357 1.99337 3.66327 \n", " -0.516414 -0.996711 -0.0919924 0.262832 0.612402 0.621834\n", " 0.193686 1.68244 -0.770028 -0.723437 -1.4868 -1.32247 \n", " 1.6567 0.0634435 -0.901968 … -0.241872 -0.0356772 -0.39826 \n", " 0.553996 -0.274515 2.21265 0.219437 2.20382 2.60902 \n", " 0.402356 1.89288 -1.13032 -0.771441 -1.96862 -1.93483 \n", " -1.07744 -1.63881 1.78016 0.96551 1.7292 1.91326 \n", " -2.21617 -2.90695 -2.55971 -0.47867 0.855389 -0.933916\n", " 1.29975 0.779828 4.12459 … 1.87358 0.737112 2.84136 \n", " -0.80833 1.44882 1.67581 -0.139063 -0.107873 0.818132\n", " -2.32469 -4.83109 -2.31796 -0.0346402 2.65564 0.591075\n", " ⋮ ⋱ \n", " 0.334995 0.0914829 1.60026 0.0937123 1.40804 1.92919 \n", " 0.390722 2.36157 -1.58383 -1.92435 -1.3291 -1.88114 \n", " 1.19218 0.845478 1.9362 … -0.0890571 1.36046 2.01013 \n", " 0.915526 0.889395 -0.606443 -0.428052 -0.817555 -1.2017 \n", " -0.778472 2.1444 1.50696 0.00689644 -1.28104 -0.141234\n", " 0.275366 -0.25866 0.416593 0.481534 0.0874531 0.316543\n", " -0.289749 -1.16146 0.550825 0.698152 0.789054 0.949917\n", " 0.439213 -0.379216 1.0951 … 0.626399 0.624574 0.8568 \n", " -0.592534 -0.235847 -1.11586 -0.601497 -0.0651787 -0.573737\n", " 0.764278 -1.86465 1.99357 3.21548 0.932773 1.97505 \n", " -0.241331 -1.92237 1.99337 0.932773 3.43991 2.8539 \n", " 0.54921 -1.72569 3.66327 1.97505 2.8539 4.85234 " ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using BenchmarkTools, Random, WoodburyMatrices\n", "\n", "Random.seed!(280)\n", "n = 1000\n", "r = 5\n", "\n", "A = Diagonal(rand(n))\n", "B = randn(n, r)\n", "D = Diagonal(rand(r))\n", "b = randn(n)\n", "# Woodbury structure: W = A + B * D * B'\n", "W = SymWoodbury(A, B, D)\n", "Wfull = Symmetric(Matrix(W)) # stored as a Matrix{Float64}" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(48200, 8000056)" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compares storage\n", "Base.summarysize(W), Base.summarysize(Wfull)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solve linear equation" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.64 MiB\n", " allocs estimate: 7\n", " --------------\n", " minimum time: 6.344 ms (0.00% GC)\n", " median time: 6.739 ms (0.00% GC)\n", " mean time: 7.086 ms (8.22% GC)\n", " maximum time: 9.931 ms (17.51% GC)\n", " --------------\n", " samples: 705\n", " evals/sample: 1" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve via Cholesky\n", "@benchmark cholesky($Wfull) \\ $b" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 75.45 KiB\n", " allocs estimate: 24\n", " --------------\n", " minimum time: 21.200 μs (0.00% GC)\n", " median time: 43.413 μs (0.00% GC)\n", " mean time: 52.796 μs (16.47% GC)\n", " maximum time: 3.121 ms (97.47% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# solve using Woodbury formula\n", "@benchmark $W \\ $b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrix-vector multiplication" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.94 KiB\n", " allocs estimate: 1\n", " --------------\n", " minimum time: 33.573 μs (0.00% GC)\n", " median time: 37.166 μs (0.00% GC)\n", " mean time: 45.578 μs (0.00% GC)\n", " maximum time: 385.600 μs (0.00% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# multiplication without using Woodbury structure\n", "@benchmark $Wfull * $b" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 24.06 KiB\n", " allocs estimate: 5\n", " --------------\n", " minimum time: 4.087 μs (0.00% GC)\n", " median time: 4.960 μs (0.00% GC)\n", " mean time: 7.489 μs (28.11% GC)\n", " maximum time: 391.177 μs (96.95% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 7" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# multiplication using Woodbury structure\n", "@benchmark $W * $b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Determinant" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 8.13 MiB\n", " allocs estimate: 7\n", " --------------\n", " minimum time: 12.483 ms (0.00% GC)\n", " median time: 13.279 ms (0.00% GC)\n", " mean time: 13.585 ms (3.69% GC)\n", " maximum time: 18.074 ms (0.00% GC)\n", " --------------\n", " samples: 368\n", " evals/sample: 1" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# determinant without using Woodbury structure\n", "@benchmark det($Wfull)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 138.05 KiB\n", " allocs estimate: 28\n", " --------------\n", " minimum time: 28.345 μs (0.00% GC)\n", " median time: 88.347 μs (0.00% GC)\n", " mean time: 109.560 μs (19.29% GC)\n", " maximum time: 4.632 ms (96.85% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 1" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# determinant using Woodbury structure\n", "@benchmark det($W)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Easy plus border\n", "\n", "Easy plus border: For $\\mathbf{A}$ pd and $\\mathbf{V}$ full row rank,\n", "$$\n", "\t\\begin{pmatrix}\n", "\t\\mathbf{A} & \\mathbf{V}^T \\\\\n", "\t\\mathbf{V} & \\mathbf{0}\n", "\t\\end{pmatrix}^{-1} = \\begin{pmatrix}\n", "\t\\mathbf{A}^{-1} - \\mathbf{A}^{-1} \\mathbf{V}^T (\\mathbf{V} \\mathbf{A}^{-1} \\mathbf{V}^T)^{-1} \\mathbf{V} \\mathbf{A}^{-1} & \\mathbf{A}^{-1} \\mathbf{V}^T (\\mathbf{V} \\mathbf{A}^{-1} \\mathbf{V}^T)^{-1} \\\\\n", "\t(\\mathbf{V} \\mathbf{A}^{-1} \\mathbf{V}^T)^{-1} \\mathbf{V} \\mathbf{A}^{-1} & - (\\mathbf{V} \\mathbf{A}^{-1} \\mathbf{V}^T)^{-1}\n", "\t\\end{pmatrix}.\n", "$$\n", "**Anyone interested writing a package?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Orthogonal matrix\n", "\n", "Orthogonal $\\mathbf{A}$: $n^2$ flops **at most**. Why? Permutation matrix, Householder matrix, Jacobi matrix, ... take less." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Toeplitz matrix\n", "\n", "Toeplitz systems (constant diagonals):\n", "$$\n", "\t\\mathbf{T} = \\begin{pmatrix}\n", "\tr_0 & r_1 & r_2 & r_3 \\\\\n", "\tr_{-1} & r_0 & r_1 & r_2 \\\\\n", "\tr_{-2} & r_{-1} & r_0 & r_1 \\\\\n", "\tr_{-3} & r_{-2} & r_{-1} & r_0\n", "\t\\end{pmatrix}.\n", "$$\n", "$\\mathbf{T} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{T}$ is pd and Toeplitz, can be solved in $O(n^2)$ flops. Durbin algorithm (Yule-Walker equation), Levinson algorithm (general $\\mathbf{b}$), Trench algorithm (inverse). These matrices occur in auto-regressive models and econometrics.\n", "\n", "* [`ToeplitzMatrices.jl`](https://github.com/JuliaMatrices/ToeplitzMatrices.jl) package can be useful." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Circulant matrix\n", "\n", "Circulant systems: Toeplitz matrix with wraparound\n", "$$\n", "\tC(\\mathbf{z}) = \\begin{pmatrix}\n", "\tz_0 & z_4 & z_3 & z_2 & z_1 \\\\\n", "\tz_1 & z_0 & z_4 & z_3 & z_2 \\\\\n", "\tz_2 & z_1 & z_0 & z_4 & z_3 \\\\\n", "\tz_3 & z_2 & z_1 & z_0 & z_4 \\\\\n", "\tz_4 & z_3 & z_2 & z_1 & z_0\n", "\t\\end{pmatrix},\n", "$$\n", "FFT type algorithms: DCT (discrete cosine transform) and DST (discrete sine transform)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vandermonde matrix\n", "\n", "Vandermonde matrix: such as in interpolation and approximation problems\n", "$$\n", "\t\\mathbf{V}(x_0,\\ldots,x_n) = \\begin{pmatrix}\n", "\t1 & 1 & \\cdots & 1 \\\\\n", "\tx_0 & x_1 & \\cdots & x_n \\\\\n", "\t\\vdots & \\vdots & & \\vdots \\\\\n", "\tx_0^n & x_1^n & \\cdots & x_n^n\n", "\t\\end{pmatrix}.\n", "$$\n", "$\\mathbf{V} \\mathbf{x} = \\mathbf{b}$ or $\\mathbf{V}^T \\mathbf{x} = \\mathbf{b}$ can be solved in $O(n^2)$ flops." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cauchy-like matrix\n", "\n", "Cauchy-like matrices:\n", "$$\n", "\t\\Omega \\mathbf{A} - \\mathbf{A} \\Lambda = \\mathbf{R} \\mathbf{S}^T,\n", "$$\n", "where $\\Omega = \\text{diag}(\\omega_1,\\ldots,\\omega_n)$ and $\\Lambda = \\text{diag}(\\lambda_1,\\ldots, \\lambda_n)$. $O(n)$ flops for LU and QR." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Structured-rank matrix\n", "\n", "Structured-rank problems: semiseparable matrices (LU and QR takes $O(n)$ flops), quasiseparable matrices, ..." ] } ], "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": "30px", "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 }