--- title: Triangular Systems subtitle: Biostat/Biomath M257 author: Dr. Hua Zhou @ UCLA date: today format: html: theme: cosmo embed-resources: true number-sections: true toc: true toc-depth: 4 toc-location: left code-fold: false jupyter: jupytext: formats: 'ipynb,qmd' text_representation: extension: .qmd format_name: quarto format_version: '1.0' jupytext_version: 1.14.5 kernelspec: display_name: Julia 1.8.5 language: julia name: julia-1.8 --- System information (for reproducibility): ```{julia} versioninfo() ``` Load packages: ```{julia} using Pkg Pkg.activate(pwd()) Pkg.instantiate() Pkg.status() ``` For the next couple of lectures, we consider computer algorithms for solving linear equations $\mathbf{A} \mathbf{x} = \mathbf{b}$, a ubiquitous task in statistics. Key idea: turning original problem into an **easy** one, e.g., triangular system. ## Lower triangular system To solve $\mathbf{A} \mathbf{x} = \mathbf{b}$, where $\mathbf{A} \in \mathbb{R}^{n \times n}$ is **lower triangular** $$ \begin{pmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}. $$ * **Forward substitution**: $$ \begin{eqnarray*} x_1 &=& b_1 / a_{11} \\ x_2 &=& (b_2 - a_{21} x_1) / a_{22} \\ x_3 &=& (b_3 - a_{31} x_1 - a_{32} x_2) / a_{33} \\ &\vdots& \\ x_n &=& (b_n - a_{n1} x_1 - a_{n2} x_2 - \cdots - a_{n,n-1} x_{n-1}) / a_{nn}. \end{eqnarray*} $$ * $1 + 3 + 5 + \cdots + (2n-1) = n^2$ flops. * $\mathbf{A}$ can be accessed by row ($ij$ loop) or column ($ji$ loop). ## Upper triangular system To solve $\mathbf{A} \mathbf{x} = \mathbf{b}$, where $\mathbf{A} \in \mathbb{R}^{n \times n}$ is upper triangular $$ \begin{pmatrix} a_{11} & \cdots & a_{1,n-1} & a_{1n} \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & a_{n-1,n-1} & a_{n-1,n} \\ 0 & 0 & 0 & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_{n-1} \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ \vdots \\ b_{n-1} \\ b_n \end{pmatrix}. $$ * **Back substitution** $$ \begin{eqnarray*} x_n &=& b_n / a_{nn} \\ x_{n-1} &=& (b_{n-1} - a_{n-1,n} x_n) / a_{n-1,n-1} \\ 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} \\ &\vdots& \\ x_1 &=& (b_1 - a_{12} x_2 - a_{13} x_3 - \cdots - a_{1,n} x_{n}) / a_{11}. \end{eqnarray*} $$ * $n^2$ flops. * $\mathbf{A}$ can be accessed by row ($ij$ loop) or column ($ji$ loop). ## Implementation * BLAS level 2 function: [trsv](http://www.netlib.org/lapack/explore-html/d6/d96/dtrsv_8f.html) (triangular solve with one right hand side). * 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). * Julia - The left divide `\` operator in Julia is used for solving linear equations or least squares problem. - If `A` is a triangular matrix, the command `A \ b` uses forward or backward substitution - 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) ```{julia} using LinearAlgebra, Random Random.seed!(123) # seed n = 5 A = randn(n, n) b = randn(n) # a random matrix A ``` ```{julia} Al = LowerTriangular(A) # does not create extra matrix ``` ```{julia} dump(Al) ``` ```{julia} Al.data ``` ```{julia} # same data pointer(Al.data), pointer(A) ``` ```{julia} Al \ b # dispatched to BLAS function for triangular solve ``` ```{julia} # or use BLAS wrapper directly BLAS.trsv('L', 'N', 'N', A, b) ``` ```{julia} ?BLAS.trsv ``` * 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!) ## Some algebraic facts of triangular system * Eigenvalues of a triangular matrix $\mathbf{A}$ are diagonal entries $\lambda_i = a_{ii}$. * Determinant $\det(\mathbf{A}) = \prod_i a_{ii}$. * The product of two upper (lower) triangular matrices is upper (lower) triangular. * The inverse of an upper (lower) triangular matrix is upper (lower) triangular. * The product of two unit upper (lower) triangular matrices is unit upper (lower) triangular. * The inverse of a unit upper (lower) triangular matrix is unit upper (lower) triangular. ## Julia types for triangular matrices [LowerTriangular](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.LowerTriangular), UnitLowerTriangular, [UpperTriangular](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.UpperTriangular), UnitUpperTriangular. ```{julia} A ``` ```{julia} LowerTriangular(A) ``` ```{julia} LinearAlgebra.UnitLowerTriangular(A) ``` ```{julia} using BenchmarkTools, LinearAlgebra, Random Random.seed!(123) # seed A = randn(1000, 1000); ``` ```{julia} # if we don't tell Julia it's triangular: O(n^3) complexity # tril(A) returns a full triangular matrix, same as Matlab @benchmark eigvals($(tril(A))) ``` ```{julia} # if we tell Julia it's triangular: O(n) complexity @benchmark eigvals($(LowerTriangular(A))) ``` ```{julia} # if we don't tell Julia it's triangular: O(n^3) complexity # tril(A) returns a full triangular matrix, same as Matlab @benchmark det($(tril(A))) ``` ```{julia} # if we tell Julia it's triangular: O(n) complexity @benchmark det($(LowerTriangular(A))) ``` ```{julia} @benchmark det($(LowerTriangular(A))) ```