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

1  Triangular systems
1.1  Lower triangular system
1.2  Upper triangular system
1.3  Implementation
1.4  Some algebraic facts of triangular system
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Triangular systems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider computer algorithms for solving linear equations $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, a ubiquitous task in statistics. \n", "\n", "Idea: turning original problem into an **easy** one, e.g., triangular system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lower triangular system\n", "\n", "To solve $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ is **lower triangular**\n", "\n", "$$\n", "\\begin{pmatrix}\n", " a_{11} & 0 & \\cdots & 0 \\\\\n", " a_{21} & a_{22} & \\cdots & 0 \\\\\n", " \\vdots & \\vdots & \\ddots & \\vdots \\\\\n", " a_{n1} & a_{n2} & \\cdots & a_{nn}\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n\n", "\\end{pmatrix} = \\begin{pmatrix}\n", "b_1 \\\\ b_2 \\\\ \\vdots \\\\ b_n\n", "\\end{pmatrix}.\n", "$$\n", "\n", "* **Forward substitution**: \n", "$$\n", "\\begin{eqnarray*}\n", " x_1 &=& b_1 / a_{11} \\\\\n", " x_2 &=& (b_2 - a_{21} x_1) / a_{22} \\\\\n", " x_3 &=& (b_3 - a_{31} x_1 - a_{32} x_2) / a_{33} \\\\\n", " &\\vdots& \\\\\n", " x_n &=& (b_n - a_{n1} x_1 - a_{n2} x_2 - \\cdots - a_{n,n-1} x_{n-1}) / a_{nn}.\n", "\\end{eqnarray*}\n", "$$\n", "\n", "* $1 + 3 + 5 + \\cdots + (2n-1) = n^2$ flops. \n", "\n", "* $\\mathbf{A}$ can be accessed by row ($ij$ loop) or column ($ji$ loop)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Upper triangular system\n", "\n", "To solve $\\mathbf{A} \\mathbf{x} = \\mathbf{b}$, where $\\mathbf{A} \\in \\mathbb{R}^{n \\times n}$ is upper triangular \n", "$$\n", "\\begin{pmatrix}\n", " a_{11} & \\cdots & a_{1,n-1} & a_{1n} \\\\\n", " \\vdots & \\ddots & \\vdots & \\vdots \\\\\n", " 0 & \\cdots & a_{n-1,n-1} & a_{n-1,n} \\\\\n", " 0 & 0 & 0 & a_{nn}\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "x_1 \\\\ \\vdots \\\\ x_{n-1} \\\\ x_n\n", "\\end{pmatrix} = \\begin{pmatrix}\n", "b_1 \\\\ \\vdots \\\\ b_{n-1} \\\\ b_n\n", "\\end{pmatrix}.\n", "$$\n", "\n", "* **Back substitution** \n", "$$\n", "\\begin{eqnarray*}\n", " x_n &=& b_n / a_{nn} \\\\\n", " x_{n-1} &=& (b_{n-1} - a_{n-1,n} x_n) / a_{n-1,n-1} \\\\\n", " x_{n-2} &=& (b_{n-2} - a_{n-2,n-1} x_{n-1} - a_{n-2,n} x_n) / a_{n-2,n-2} \\\\\n", " &\\vdots& \\\\\n", " x_1 &=& (b_1 - a_{12} x_2 - a_{13} x_3 - \\cdots - a_{1,n} x_{n}) / a_{11}.\n", "\\end{eqnarray*}\n", "$$\n", "\n", "* $n^2$ flops.\n", "\n", "* $\\mathbf{A}$ can be accessed by row ($ij$ loop) or column ($ji$ loop)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "\n", "* BLAS level 2 function: [`?trsv`](http://www.netlib.org/lapack/explore-html/d6/d96/dtrsv_8f.html) (triangular solve with one right hand side).\n", "\n", "* BLAS level 3 function: [`?trsm`](http://www.netlib.org/lapack/explore-html/de/da7/dtrsm_8f.html) (matrix triangular solve, i.e., multiple right hand sides).\n", "\n", "* Julia \n", " - The left divide `\\` operator in Julia is used for solving linear equations or least squares problem. \n", " - If `A` is a triangular matrix, the command `A \\ b` uses forward or backward substitution \n", " - Or we can call the BLAS wrapper functions directly: [`trsv!`](https://docs.julialang.org/en/stable/stdlib/linalg/?highlight=blas#Base.LinAlg.BLAS.trsv!), [`trsv`](https://docs.julialang.org/en/stable/stdlib/linalg/?highlight=blas#Base.LinAlg.BLAS.trsv), [`trsm`](https://docs.julialang.org/en/stable/stdlib/linalg/?highlight=blas#Base.LinAlg.BLAS.trsm), [`trsm!`](https://docs.julialang.org/en/stable/stdlib/linalg/?highlight=blas#Base.LinAlg.BLAS.trsm!)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Float64,2}:\n", " 1.19027 -0.664713 -0.339366 0.368002 -0.979539\n", " 2.04818 0.980968 -0.843878 -0.281133 0.260402\n", " 1.14265 -0.0754831 -0.888936 -0.734886 -0.468489\n", " 0.459416 0.273815 0.327215 -0.71741 -0.880897\n", " -0.396679 -0.194229 0.592403 -0.77507 0.277726" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "srand(123) # seed\n", "n = 5\n", "A = randn(n, n)\n", "b = randn(n)\n", "# a random matrix\n", "A" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Float64,2}:\n", " 1.19027 0.0 0.0 0.0 0.0 \n", " 2.04818 0.980968 0.0 0.0 0.0 \n", " 1.14265 -0.0754831 -0.888936 0.0 0.0 \n", " 0.459416 0.273815 0.327215 -0.71741 0.0 \n", " -0.396679 -0.194229 0.592403 -0.77507 0.277726" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tril(A) # create another triangular matrix" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 LowerTriangular{Float64,Array{Float64,2}}:\n", " 1.19027 ⋅ ⋅ ⋅ ⋅ \n", " 2.04818 0.980968 ⋅ ⋅ ⋅ \n", " 1.14265 -0.0754831 -0.888936 ⋅ ⋅ \n", " 0.459416 0.273815 0.327215 -0.71741 ⋅ \n", " -0.396679 -0.194229 0.592403 -0.77507 0.277726" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Al = LowerTriangular(A) # does not create extra matrix" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1-element Array{Symbol,1}:\n", " :data" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fieldnames(Al)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Float64,2}:\n", " 1.19027 -0.664713 -0.339366 0.368002 -0.979539\n", " 2.04818 0.980968 -0.843878 -0.281133 0.260402\n", " 1.14265 -0.0754831 -0.888936 -0.734886 -0.468489\n", " 0.459416 0.273815 0.327215 -0.71741 -0.880897\n", " -0.396679 -0.194229 0.592403 -0.77507 0.277726" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Al.data" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 1.28031\n", " -4.48541\n", " 5.32613\n", " 0.44682\n", " -3.09169" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tril(A) \\ b # dispatched to LowerTriangular(A) \\ b" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) at linalg/generic.jl:820" ], "text/plain": [ "\\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in Base.LinAlg at linalg/generic.jl:820" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check source code\n", "@which tril(A) \\ b" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 1.28031\n", " -4.48541\n", " 5.32613\n", " 0.44682\n", " -3.09169" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "LowerTriangular(A) \\ b # dispatched to A_ldiv_B" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\\(A::Union{LowerTriangular, UpperTriangular}, B::AbstractArray{T,1} where T) at linalg/triangular.jl:1624" ], "text/plain": [ "\\(A::Union{LowerTriangular, UpperTriangular}, B::AbstractArray{T,1} where T) in Base.LinAlg at linalg/triangular.jl:1624" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check source code\n", "@which LowerTriangular(A) \\ b" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 1.28031\n", " -4.48541\n", " 5.32613\n", " 0.44682\n", " -3.09169" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# or use BLAS wrapper directly\n", "LinAlg.BLAS.trsv('L', 'N', 'N', A, b)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/markdown": [ "```\n", "trsv(ul, tA, dA, A, b)\n", "```\n", "\n", "Returns the solution to `A*x = b` or one of the other two variants determined by [`tA`](@ref stdlib-blas-trans) and [`ul`](@ref stdlib-blas-uplo). [`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or are assumed to be all ones.\n" ], "text/plain": [ "```\n", "trsv(ul, tA, dA, A, b)\n", "```\n", "\n", "Returns the solution to `A*x = b` or one of the other two variants determined by [`tA`](@ref stdlib-blas-trans) and [`ul`](@ref stdlib-blas-uplo). [`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or are assumed to be all ones.\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "?LinAlg.BLAS.trsv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Some other BLAS functions for triangular systems: [`trmv`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.BLAS.trmv), [`trmv!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.BLAS.trmv!), [`trmm`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.BLAS.trmm), [`trmm!`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.BLAS.trmm!)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some algebraic facts of triangular system\n", "\n", "* Eigenvalues of a triangular matrix $\\mathbf{A}$ are diagonal entries $\\lambda_i = a_{ii}$. \n", "\n", "* Determinant $\\det(\\mathbf{A}) = \\prod_i a_{ii}$.\n", "\n", "* The product of two upper (lower) triangular matrices is upper (lower) triangular.\n", "\n", "* The inverse of an upper (lower) triangular matrix is upper (lower) triangular.\n", "\n", "* The product of two unit upper (lower) triangular matrices is unit upper (lower) triangular.\n", "\n", "* The inverse of a unit upper (lower) triangular matrix is unit upper (lower) triangular." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 15.55 MiB\n", " allocs estimate: 22\n", " --------------\n", " minimum time: 45.093 ms (0.00% GC)\n", " median time: 50.140 ms (3.11% GC)\n", " mean time: 50.127 ms (2.49% GC)\n", " maximum time: 55.333 ms (0.00% GC)\n", " --------------\n", " samples: 100\n", " evals/sample: 1" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using BenchmarkTools\n", "\n", "srand(123) # seed\n", "A = randn(1000, 1000)\n", "\n", "# if we don't tell Julia it's triangular\n", "@benchmark eigvals(tril(A))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 8.02 KiB\n", " allocs estimate: 4\n", " --------------\n", " minimum time: 2.232 μs (0.00% GC)\n", " median time: 2.728 μs (0.00% GC)\n", " mean time: 3.355 μs (10.53% GC)\n", " maximum time: 206.033 μs (93.15% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 9" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# if we tell Julia it's triangular\n", "@benchmark eigvals(LowerTriangular(A))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.64 MiB\n", " allocs estimate: 7\n", " --------------\n", " minimum time: 2.564 ms (0.00% GC)\n", " median time: 3.212 ms (0.00% GC)\n", " mean time: 4.040 ms (21.96% GC)\n", " maximum time: 10.679 ms (31.87% GC)\n", " --------------\n", " samples: 1233\n", " evals/sample: 1" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# if we don't tell Julia it's triangular\n", "@benchmark det(tril(A))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 8.03 KiB\n", " allocs estimate: 5\n", " --------------\n", " minimum time: 2.589 μs (0.00% GC)\n", " median time: 3.088 μs (0.00% GC)\n", " mean time: 3.684 μs (12.47% GC)\n", " maximum time: 250.764 μs (96.42% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 9" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# if we tell Julia it's triangular\n", "@benchmark det(LowerTriangular(A))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 0.6.2\n", "Commit d386e40c17 (2017-12-13 18:08 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", " BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n", " LAPACK: libopenblas64_\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-3.9.1 (ORCJIT, skylake)\n" ] } ], "source": [ "versioninfo()" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "103px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }