{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## 1D finite differences\n", "Jeff Fessler, University of Michigan\n", "2019-01-31 \n", "2020-01-18 Julia 1.3.1, LinearMapsAA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first-order finite difference operation is important for regularization. \n", "It is especially important for it to be memory and compute efficient\n", "for 2D and higher dimensional problems,\n", "but for simplicity this notebook first explores it for 1D problems.\n", "\n", "If `x` is a 1D array of length $N$ and `y` is its first-order finite difference\n", "then \n", "`y[n] = x[n+1] - x[n]` for $n=1,\\ldots,N-1$.\n", "\n", "There are surprisingly many ways to compute `y` from `x`\n", "in a high-level language like Julia.\n", "After loading the key packages,\n", "the next cell gives a table of some possible methods." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using BenchmarkTools\n", "# using LinearAlgebra: diagm - far too slow to even illustrate\n", "using Printf\n", "# using DSP # for conv() - also too slow to bother\n", "# using FastConv - currently not loading properly so excluded. Will be slower anyway.\n", "# using DirectConvolution - ditto.\n", "\n", "do_time = ENV[\"HOST\"] == \"jf28.local\" # set to true to run timing tests\n", "do_time = false;" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7×2 Array{Any,2}:\n", " #3 \"myloop()\" \n", " #4 \"diff()\" \n", " #5 \"comprehension\" \n", " #7 \"compreh. inbounds\"\n", " #9 \"circshift\" \n", " #10 \"index\" \n", " #11 \"index views\" " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# here it is written as a loop as some Fortran or C programmers might do\n", "function myloop(x::AbstractVector{<:Number})\n", " N = length(x)\n", " y = similar(x, N-1)\n", " for n=1:(N-1)\n", " @inbounds y[n] = x[n+1] - x[n]\n", " end\n", " return y\n", "end\n", "\n", "# the remaining methods all use built-in Julia functions\n", "tryfun = [\n", " x -> myloop(x) \"myloop()\";\n", " x -> diff(x) \"diff()\";\n", " x -> [x[n+1]-x[n] for n=1:(length(x)-1)] \"comprehension\";\n", " x -> (@inbounds [x[n+1]-x[n] for n=1:(length(x)-1)]) \"compreh. inbounds\";\n", " x -> (circshift(x, -1) - x)[1:(length(x)-1)] \"circshift\";\n", " x -> x[2:end] - x[1:(end-1)] \"index\";\n", "# x -> (@inbounds x[2:end] - x[1:(end-1)]) \"index inbound\";\n", "# x -> (x[2:N] - x[1:(N-1)]) \"index N\";\n", "# x -> (@inbounds x[2:N] - x[1:(N-1)]) \"index inbound N\";\n", " x -> (@views x[2:end] - x[1:(end-1)]) \"index views\";\n", "# x -> conv(x, eltype(x).([1, -1]))[2:end-1] \"conv\"; # slow!\n", "]\n", "nfun = size(tryfun,1)\n", "tryfun" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next cell creates some test data\n", "and then verifies that all of the `nfun` methods above\n", "produce the same output\n", "to within numerical precision.\n", "We use `Float32` rather than the default 64-bit float\n", "because single precision suffices\n", "for most imaging applications\n", "and runs faster because of less memory." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "N = 2^15 # adjust this based on your computer speed\n", "x = randn(Float32, N) # faster than Float64\n", "\n", "y = diff(x)\n", "for it = 1:size(tryfun,1) # warm-up\n", " fun = tryfun[it,1]\n", " @assert isapprox(y, fun(x))\n", "end\n", "label = tryfun[:,2];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use `@btime` macro in `BenchmarkTools` package to time these `nfun` methods." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "if do_time\n", " for it = 1:size(tryfun,1)\n", " @printf(\"%20s \", tryfun[it,2])\n", " @btime y = $tryfun[$it,1]($x)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results of course depend on the computer and Julia version. \n", "\n", "Here are the results for a 2017 iMac with 4.2GHz Intel Core i7 (Kaby Lake I7-7700K) with Julia 1.3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", " myloop() 7.637 μs (2 allocations: 128.08 KiB)\n", " diff() 7.500 μs (4 allocations: 128.17 KiB)\n", " comprehension 22.183 μs (4 allocations: 128.13 KiB)\n", " compreh. inbounds 29.439 μs (4 allocations: 128.13 KiB)\n", " circshift 21.929 μs (6 allocations: 384.23 KiB)\n", " index 20.490 μs (6 allocations: 384.23 KiB)\n", " index views 7.771 μs (4 allocations: 128.17 KiB)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Which approach is the most like well-written single-thread C-code?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A. `myloop(x)` \n", "B. `x[2:end] - x[1:(end-1)]` (index) \n", "C. `circshift(x ...` \n", "D. `(@views x[2:end] - x[1:(end-1)])` (index views) \n", "E. `[x[n+1]-x[n] for n=1:(length(x)-1)]` (comprehension)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sparse matrices\n", "These function calls are all linear, but they do not \"look\" like matrix operations.\n", "In contrast, if we create a sparse matrix `C` using `spdiagm`\n", "then we can make the code look a lot like the math\n", "by typing `y = C * x`" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "using SparseArrays: spdiagm\n", "# note: the extra [1:(N-1),:] indexing at the end eliminates a useless (?) row of zeros\n", "Spn = (N,T) -> spdiagm(0 => -ones(T, N-1), 1 => ones(T, N-1))[1:(N-1),:] # sparse\n", "C = Spn(length(x),eltype(x))\n", "@assert y == C*x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But using sparse matrix multiplication is much slower than using `diff`:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "if do_time\n", " @btime $C * $x\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On same computer:\n", "```\n", "63.513 μs (2 allocations: 128.08 KiB)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### LinearMaps\n", "\n", "The `LinearMapsAA` package lets us do the same thing with any (linear) function!" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearMapAA: 32767×32768\n", "(none = nothing,)\n", "LinearMaps.FunctionMap{Float32}(#23, 32767, 32768; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearMapsAA\n", "Cdiff = LinearMapAA(x -> diff(x), (N-1, N))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# it works as you would hope:\n", "@assert y == Cdiff*x" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearMapAA: 4×5\n", "(none = nothing,)\n", "LinearMaps.FunctionMap{Float32}(#25, 4, 5; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)\n" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# we can examine the elements for small sizes\n", "Ctiny = LinearMapAA(x -> diff(x), (4, 5))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×5 Array{Int64,2}:\n", " -1 1 0 0 0\n", " 0 -1 1 0 0\n", " 0 0 -1 1 0\n", " 0 0 0 -1 1" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Int.(Matrix(Ctiny))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element Array{Int64,1}:\n", " 1\n", " -1\n", " 0\n", " 0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# How does Matrix() work here? Think about this:\n", "Ctiny * [0, 1, 0, 0, 0]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element Array{Int64,1}:\n", " 2\n", " 1\n", " 3\n", " -7" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# It acts like a matrix for multiplication:\n", "Ctiny * [1, 3, 4, 7, 0]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element Array{Int64,1}:\n", " 1\n", " 1\n", " 1\n", " 1" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Ctiny * (1:5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Timing the LinearMap (and sparse) versions" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# make LinearMap from each of the functions\n", "\n", "DD = similar(Array{Any}, nfun+1)\n", "\n", "for it = 1:size(tryfun,1) # warm-up\n", " fun = tryfun[it,1]\n", " DD[it] = LinearMapAA(fun, (N-1, N))\n", " @assert isapprox(y, DD[it] * x)\n", "end\n", "\n", "# include sparse array in the collection\n", "if true # spdiagm\n", " label = [tryfun[:,2]..., \"spdiagm\"]\n", " DD[length(label)] = C\n", " @assert isapprox(y, DD[length(label)] * x)\n", "end\n", "\n", "if do_time\n", " for it = 1:length(label)\n", " @printf(\"%20s \", label[it])\n", " @btime y = $DD[$it] * $x\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Results on the same computer: \n", "```\n", " myloop() 7.757 μs (2 allocations: 128.08 KiB)\n", " diff() 7.647 μs (4 allocations: 128.17 KiB)\n", " comprehension 22.295 μs (4 allocations: 128.13 KiB)\n", " compreh. inbounds 29.454 μs (4 allocations: 128.13 KiB)\n", " circshift 21.520 μs (6 allocations: 384.23 KiB)\n", " index 20.608 μs (6 allocations: 384.23 KiB)\n", " index views 7.730 μs (4 allocations: 128.17 KiB)\n", " spdiagm 64.011 μs (2 allocations: 128.08 KiB)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Take away:\n", "Using `spdiagm` to perform `C*x` is $\\approx 10 \\times$ slower\n", "than using the `LinearMapAA` with the fastest function! \n", "Also, `LinearMapAA` requires *very little* overhead." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Transpose\n", "But we are not done yet.\n", "Our regularizer gradient needs multiply operations\n", "with both `C` and `C'`\n", "and we have `C'` for the sparse array,\n", "but not yet for the `LinearMapAA` version.\n", "First look at the transpose to see what it should be:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×4 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:\n", " -1 0 0 0\n", " 1 -1 0 0\n", " 0 1 -1 0\n", " 0 0 1 -1\n", " 0 0 0 1" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Int.(Matrix(Ctiny))'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the variable type above is `Adjoint` that is a generalization of transpose: \n", "https://en.wikipedia.org/wiki/Adjoint\n", "\n", "We can define many Julia functions to perform `C' * y`.\n", "Here is a C/Fortran style loop:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "myadj (generic function with 1 method)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function myadj(y::AbstractVector{<:Number})\n", " N = length(y)\n", " x = similar(y, N+1)\n", " x[1] = -y[1]\n", " @views x[2:end-1] .= y[1:(end-1)] - y[2:end]\n", " x[N+1] = y[N]\n", " return x\n", "end" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# the remaining methods all use built-in Julia functions\n", "adjfun = [\n", " y -> myadj(y) \"myadj()\";\n", " y -> [-y[1]; -diff(y); y[end]] \"diff()\";\n", " y -> [-y[1]; [y[n]-y[n+1] for n=1:(length(y)-1)]; y[end]] \"comprehension\";\n", " y -> [-y[1]; (@views y[1:end-1] - y[2:end]); y[end]] \"index views\";\n", "# y -> conv(y, eltype(y).([1, -1]))[???] \"conv\"; # slow, not done!\n", "]\n", "nadj = size(adjfun,1)\n", "\n", "z = C' * y\n", "for it = 1:size(adjfun,1) # warm-up and check\n", " fun = adjfun[it,1]\n", " @assert isapprox(z, fun(y))\n", "end\n", "adjlabel = tryfun[:,2];" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# time these methods\n", "if do_time\n", " for it = 1:size(adjfun,1)\n", " @printf(\"%20s \", adjfun[it,2])\n", " @btime z = $adjfun[$it,1]($y)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On the same computer the adjoint timing results were:\n", "```\n", " myadj() 15.089 μs (7 allocations: 256.30 KiB)\n", " diff() 25.894 μs (39 allocations: 385.31 KiB)\n", " comprehension 31.864 μs (37 allocations: 257.19 KiB)\n", " index views 17.372 μs (37 allocations: 257.23 KiB)\n", "```\n", "\n", "Here, index views and `myadj()` are the two fastest." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now make a new `LinearMapAA` that has both forward and adjoint multiplication. \n", "Use the `@views` functions for both for simplicity." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearMapAA: 4×5\n", "(none = nothing,)\n", "LinearMaps.FunctionMap{Float32}(#38, #39, 4, 5; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)\n" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ClmN = N -> LinearMapAA(\n", " x -> (@views x[2:end] - x[1:(end-1)]), # forward function\n", " y -> [-y[1]; (@views y[1:end-1] - y[2:end]); y[end]], (N-1, N)) # adjoint function\n", "Clm5 = ClmN(5)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearMapAA: 5×4\n", "(none = nothing,)\n", "LinearMaps.TransposeMap{Float32,LinearMaps.FunctionMap{Float32,var\"#38#41\",var\"#39#42\"}}(LinearMaps.FunctionMap{Float32}(#38, #39, 4, 5; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false))\n" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Clm5'" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4×5 Array{Int64,2}:\n", " -1 1 0 0 0\n", " 0 -1 1 0 0\n", " 0 0 -1 1 0\n", " 0 0 0 -1 1" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "5×4 Array{Int64,2}:\n", " -1 0 0 0\n", " 1 -1 0 0\n", " 0 1 -1 0\n", " 0 0 1 -1\n", " 0 0 0 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# check for adjoint consistency\n", "display(Int.(Matrix(Clm5)))\n", "display(Int.(Matrix(Clm5')))\n", "@assert Matrix(Clm5)' == Matrix(Clm5')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Timing test to compare `LinearMapAA` versus `SparseArray`" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "Clm = ClmN(N)\n", "@assert Clm * x == C * x\n", "@assert Clm' * y == C' * y\n", "if do_time\n", " @printf(\"%22s \", \"linearmap C\")\n", " @btime $Clm * $x;\n", " @printf(\"%22s \", \"linearmap C'\")\n", " @btime $Clm' * $y\n", " @printf(\"%22s \", \"sparse C\")\n", " @btime $C * $x\n", " @printf(\"%22s \", \"sparse C'\")\n", " @btime $C' * $y\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Timing results on same computer:\n", "```\n", " linearmap C 7.848 μs (4 allocations: 128.17 KiB)\n", " linearmap C' 18.617 μs (42 allocations: 257.39 KiB)\n", " sparse C 63.818 μs (2 allocations: 128.08 KiB)\n", " sparse C' 66.526 μs (3 allocations: 128.09 KiB)\n", "```\n", "Using a `LinearMapAA` is much faster than using a `SparseArray` in both directions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check adjoint consistency for large LinearMaps" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "u = randn(size(C,1))\n", "v = randn(size(C,2))\n", "@assert u'*(C*v) ≈ (C' * u)' * v" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Illustrate Array{Any}" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 2.0\n", " 2.0\n", " 2.0\n", " 2.0\n", " 2.0" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = ones(5); b = 2*ones(5)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×2 Array{Float64,2}:\n", " 1.0 2.0\n", " 1.0 2.0\n", " 1.0 2.0\n", " 1.0 2.0\n", " 1.0 2.0" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[a b]" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element Array{Any,1}:\n", " [1.0, 1.0, 1.0, 1.0, 1.0]\n", " [2.0, 2.0, 2.0, 2.0, 2.0]\n", " \"hi\" \n", " var\"#43#44\"() " ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr = [a, b, \"hi\", x -> x^2]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×2 Array{Any,2}:\n", " 2 \"hi\"" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[2 \"hi\"]" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "data": { "text/plain": [ "1×2 Array{Any,2}:\n", " [1.0, 1.0, 1.0, 1.0, 1.0] \"hi\"" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[[a] \"hi\"]" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "# code from xijia quan in w19 - code for prof. fessler to study later.\n", "# N-N first and second order\n", "#=\n", "using Distributed\n", "a = randn(100,100)\n", "@time begin\n", " tmp_C1 = diff(a,dims=1)\n", " tmp_C0 = a[1,:] - a[end,:]\n", " C1 = vcat(tmp_C1, tmp_C0')\n", "end # 1st order\n", "\n", "@time begin\n", " tmp_C2 = -diff(C1,dims=1)\n", " tmp_C20 = C1[end,:] - C1[1,:]\n", " C2 = vcat(tmp_C2, tmp_C20')\n", "end # 2rd order\n", "\n", "# slow\n", "C3 = Matrix{Float64}(undef,size(a))\n", "@time begin\n", " @sync @distributed for i=1:size(a,2)\n", " tmp_C3 = diff(a[:,i])\n", " tmp_C30 = a[1,i] - a[end,i]\n", " C3[:,i] = vcat(tmp_C3, tmp_C30)\n", " end # 1st order\n", "end\n", "\n", "C4 = Matrix{Float64}(undef,size(a))\n", "@time begin\n", " @sync @distributed for i=1:size(a,2)\n", " tmp_C4 = -diff(C3[:,i])\n", " tmp_C40 = C3[end,i] - C3[1,i]\n", " C4[:,i] = vcat(tmp_C4, tmp_C40)\n", " end # 2rd order\n", "end\n", "=#" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.3.0", "language": "julia", "name": "julia-1.3" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.3.0" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false } }, "nbformat": 4, "nbformat_minor": 2 }