{ "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
1.5  Julia types for triangular matrices
" ] }, { "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": [ "# 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/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trsv!), [`trsv`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trsv), [`trsm!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trsm!), [`trsm`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trsm)" ] }, { "cell_type": "code", "execution_count": 2, "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": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra, Random\n", "\n", "Random.seed!(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": 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": [ "(:data,)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fieldnames(typeof(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": [ "(Ptr{Float64} @0x000000010f12c170, Ptr{Float64} @0x000000010f12c170)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# same data\n", "pointer(Al.data), pointer(A)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 1.28031359532452 \n", " -4.485407033333146 \n", " 5.326125412123139 \n", " 0.4468195081389211\n", " -3.091688163812573 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Al \\ b # dispatched to BLAS function for triangular solve" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 1.28031359532452 \n", " -4.485407033333146\n", " 5.326125412123139\n", " 0.446819508138921\n", " -3.091688163812573" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# or use BLAS wrapper directly\n", "BLAS.trsv('L', 'N', 'N', A, b)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "\\begin{verbatim}\n", "trsv(ul, tA, dA, A, b)\n", "\\end{verbatim}\n", "Return the solution to \\texttt{A*x = b} or one of the other two variants determined by \\href{@ref stdlib-blas-trans}{\\texttt{tA}} and \\href{@ref stdlib-blas-uplo}{\\texttt{ul}}. \\href{@ref stdlib-blas-diag}{\\texttt{dA}} determines if the diagonal values are read or are assumed to be all ones.\n", "\n" ], "text/markdown": [ "```\n", "trsv(ul, tA, dA, A, b)\n", "```\n", "\n", "Return 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": [ "\u001b[36m trsv(ul, tA, dA, A, b)\u001b[39m\n", "\n", " Return the solution to \u001b[36mA*x = b\u001b[39m or one of the other two variants determined\n", " by \u001b[36mtA\u001b[39m and \u001b[36mul\u001b[39m. \u001b[36mdA\u001b[39m determines if the diagonal values are read or are assumed\n", " to be all ones." ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "?BLAS.trsv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Some other BLAS functions for triangular systems: [`trmv`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trmv), [`trmv!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trmv!), [`trmm`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.BLAS.trmm), [`trmm!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.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": "markdown", "metadata": {}, "source": [ "## Julia types for triangular matrices\n", "\n", "[LowerTriangular](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.LowerTriangular), UnitLowerTriangular, \n", "[UpperTriangular](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.UpperTriangular), UnitUpperTriangular. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "using BenchmarkTools, LinearAlgebra, Random\n", "\n", "Random.seed!(123) # seed\n", "A = randn(1000, 1000);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# if we don't tell Julia it's triangular: O(n^3) complexity\n", "# tril(A) returns a full triangular matrix, same as Matlab\n", "@benchmark eigvals(tril($A))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.94 KiB\n", " allocs estimate: 1\n", " --------------\n", " minimum time: 1.676 μs (0.00% GC)\n", " median time: 2.282 μs (0.00% GC)\n", " mean time: 3.285 μs (28.33% GC)\n", " maximum time: 4.328 ms (99.91% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 10" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# if we tell Julia it's triangular: O(n) complexity\n", "@benchmark eigvals(LowerTriangular($A))" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.64 MiB\n", " allocs estimate: 4\n", " --------------\n", " minimum time: 2.024 ms (0.00% GC)\n", " median time: 2.134 ms (0.00% GC)\n", " mean time: 2.428 ms (19.43% GC)\n", " maximum time: 45.379 ms (96.47% GC)\n", " --------------\n", " samples: 2051\n", " evals/sample: 1" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# if we don't tell Julia it's triangular: O(n^3) complexity\n", "# tril(A) returns a full triangular matrix, same as Matlab\n", "@benchmark det(tril($A))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.95 KiB\n", " allocs estimate: 2\n", " --------------\n", " minimum time: 2.025 μs (0.00% GC)\n", " median time: 2.595 μs (0.00% GC)\n", " mean time: 4.003 μs (29.86% GC)\n", " maximum time: 4.951 ms (99.90% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 9" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# if we tell Julia it's triangular: O(n) complexity\n", "@benchmark det(LowerTriangular($A))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 7.97 KiB\n", " allocs estimate: 3\n", " --------------\n", " minimum time: 2.189 μs (0.00% GC)\n", " median time: 2.858 μs (0.00% GC)\n", " mean time: 4.186 μs (29.45% GC)\n", " maximum time: 5.243 ms (99.89% GC)\n", " --------------\n", " samples: 10000\n", " evals/sample: 9" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark det(LowerTriangular(A))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Array{Float64,2}:\n", " 0.505563 0.490008 0.364867 0.064762 0.537664 \n", " 0.372991 0.1636 0.301699 0.614327 0.131543 \n", " 0.00596027 0.43842 0.607055 0.342615 0.0390959\n", " 0.216016 0.891323 0.820768 0.449411 0.31552 \n", " 0.944432 0.446089 0.967155 0.00461037 0.275575 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = rand(5, 5)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 LowerTriangular{Float64,Array{Float64,2}}:\n", " 0.505563 ⋅ ⋅ ⋅ ⋅ \n", " 0.372991 0.1636 ⋅ ⋅ ⋅ \n", " 0.00596027 0.43842 0.607055 ⋅ ⋅ \n", " 0.216016 0.891323 0.820768 0.449411 ⋅ \n", " 0.944432 0.446089 0.967155 0.00461037 0.275575" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "LowerTriangular(A)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 UnitLowerTriangular{Float64,Array{Float64,2}}:\n", " 1.0 ⋅ ⋅ ⋅ ⋅ \n", " 0.372991 1.0 ⋅ ⋅ ⋅ \n", " 0.00596027 0.43842 1.0 ⋅ ⋅ \n", " 0.216016 0.891323 0.820768 1.0 ⋅ \n", " 0.944432 0.446089 0.967155 0.00461037 1.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "LinearAlgebra.UnitLowerTriangular(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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": "103px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "skip_h1_title": true, "threshold": 4, "toc_cell": true, "toc_position": { "height": "352.28570556640625px", "left": "0px", "right": "717.1428833007813px", "top": "138.7142791748047px", "width": "160.57142639160156px" }, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }