Triangular Systems of Linear Equations

Advanced Statistical Computing

Joong-Ho Won

Seoul National University

October 2023

versioninfo()

Julia Version 1.9.3
Commit bed2cd540a1 (2023-08-24 14:43 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (x86_64-apple-darwin22.4.0)
  CPU: 8 × Intel(R) Core(TM) i5-8279U CPU @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 2 on 8 virtual cores

using Pkg
Pkg.activate(pwd())
Pkg.instantiate()
Pkg.status()

Status `~/Dropbox/class/M1399.000200/2023/M1300_000200-2023fall/lectures/06-trisys/Project.toml`
  [6e4b80f9] BenchmarkTools v1.3.2

Triangular systems

We consider computer algorithms for solving linear equations \(\mathbf{A} \mathbf{x} = \mathbf{b}\), a ubiquitous task in statistics.

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{bmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}. \]

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{bmatrix} 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{bmatrix} \begin{bmatrix} x_1 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix} = \begin{bmatrix} b_1 \\ \vdots \\ b_{n-1} \\ b_n \end{bmatrix}. \]

Back substitution (backsolve): \[ \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: BLAS.trsv (triangular solve with one right hand side: \(\mathbf{A}x=b\)).
BLAS level 3 function: BLAS.trsm (matrix triangular solve, i.e., multiple right hand sides: \(\mathbf{A}\mathbf{X}=\alpha\mathbf{C}\)).

The BLAS triangular system solve is done in place, i.e., \(\mathbf{b}\) is overwritten by \(\mathbf{x}\).

    # forward substitution
    for i=1:n
        for j=1:i-1
            b[i] -= A[i, j] * b[j]
        end
    end
    # backsolve
    for i=n:-1:1
        for j=i+1:n
            b[i] -= A[i, j] * b[j]
        end
        b[i] /= A[i, i]
    end

In 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
- Imagine \(\frac{b}{A}=A^{-1}b\) to memorize.
Or we can call the BLAS wrapper functions directly: trsv!, trsv, trsm!, trsm

using LinearAlgebra, Random

Random.seed!(123) # seed
n = 5
A = randn(n, n)
b = randn(n)
# a random matrix
A

5×5 Matrix{Float64}:
  0.808288   0.229819    1.21928    -0.890077   2.00811
 -1.12207   -0.421769    0.292914    0.854242   0.76503
 -1.10464   -1.35559    -0.0311481   0.341782   0.180254
 -0.416993   0.0694591   0.315833   -0.31887    2.02891
  0.287588  -0.117323   -2.16238    -0.337454  -1.08822

Al = LowerTriangular(A) # does not create an extra matrix

5×5 LowerTriangular{Float64, Matrix{Float64}}:
  0.808288    ⋅           ⋅           ⋅          ⋅ 
 -1.12207   -0.421769     ⋅           ⋅          ⋅ 
 -1.10464   -1.35559    -0.0311481    ⋅          ⋅ 
 -0.416993   0.0694591   0.315833   -0.31887     ⋅ 
  0.287588  -0.117323   -2.16238    -0.337454  -1.08822

fieldnames(typeof(Al))

(:data,)

Al.data

5×5 Matrix{Float64}:
  0.808288   0.229819    1.21928    -0.890077   2.00811
 -1.12207   -0.421769    0.292914    0.854242   0.76503
 -1.10464   -1.35559    -0.0311481   0.341782   0.180254
 -0.416993   0.0694591   0.315833   -0.31887    2.02891
  0.287588  -0.117323   -2.16238    -0.337454  -1.08822

# same data
pointer(Al.data), pointer(A)

(Ptr{Float64} @0x0000000103a08b40, Ptr{Float64} @0x0000000103a08b40)

Al \ b # dispatched to BLAS function for triangular solve

5-element Vector{Float64}:
    0.8706777634658246
   -2.6561704782961275
   79.95736189623058
   74.31241954954415
 -181.43591336155936

# or use BLAS wrapper directly
BLAS.trsv('L', 'N', 'N', A, b)

5-element Vector{Float64}:
    0.8706777634658246
   -2.6561704782961275
   79.95736189623058
   74.31241954954415
 -181.43591336155936

?BLAS.trsv

trsv(ul, tA, dA, A, b)

Return the solution to A*x = b or one of the other two variants determined by tA and ul. dA determines if the diagonal values are read or are assumed to be all ones.

Some other BLAS functions for triangular systems: trmv, trmv!, trmm, 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, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular.

using BenchmarkTools, LinearAlgebra, Random

Random.seed!(123) # seed

A = rand(5, 5)

5×5 Matrix{Float64}:
 0.700342    0.685739   0.586295   0.0330756   0.0579227
 0.228912    0.495799   0.166108   0.328311    0.693948
 0.0681775   0.1248     0.323899   0.00168736  0.929565
 0.575665    0.891646   0.690664   0.703312    0.858179
 0.00944764  0.0463539  0.0620355  0.732952    0.136627

LowerTriangular(A)

5×5 LowerTriangular{Float64, Matrix{Float64}}:
 0.700342     ⋅          ⋅          ⋅         ⋅ 
 0.228912    0.495799    ⋅          ⋅         ⋅ 
 0.0681775   0.1248     0.323899    ⋅         ⋅ 
 0.575665    0.891646   0.690664   0.703312   ⋅ 
 0.00944764  0.0463539  0.0620355  0.732952  0.136627

LinearAlgebra.UnitLowerTriangular(A)

5×5 UnitLowerTriangular{Float64, Matrix{Float64}}:
 1.0          ⋅          ⋅          ⋅         ⋅ 
 0.228912    1.0         ⋅          ⋅         ⋅ 
 0.0681775   0.1248     1.0         ⋅         ⋅ 
 0.575665    0.891646   0.690664   1.0        ⋅ 
 0.00944764  0.0463539  0.0620355  0.732952  1.0

A = randn(1000, 1000);
# 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))

BenchmarkTools.Trial: 85 samples with 1 evaluation.
 Range (min … max):  56.034 ms … 77.472 ms  ┊ GC (min … max): 0.00% … 3.46%
 Time  (median):     57.644 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   59.215 ms ±  4.118 ms  ┊ GC (mean ± σ):  1.53% ± 2.13%
  ▂▄█▅ ▇    ▄                                                  
  ███████▆█▆█▅▅▃▁▃▁▁▁▁▁▁▁▁▁▃▃▁▁▁▁▃▁▅▁▃▁▅▁▃▃▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▃▁▃ ▁
  56 ms           Histogram: frequency by time        72.3 ms <
 Memory estimate: 15.56 MiB, allocs estimate: 16.

# if we tell Julia it's triangular: O(n) complexity
@benchmark eigvals(LowerTriangular($A))

BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.603 μs … 880.045 μs  ┊ GC (min … max):  0.00% … 99.05%
 Time  (median):     2.251 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   3.569 μs ±  28.734 μs  ┊ GC (mean ± σ):  30.92% ±  3.84%
  ▃▄▅▅▄▅▇▇██▇▆▅▃▁▁                          ▁▁▂▂              ▂
  ████████████████▇▇▄▄▄▅▆▇▅▆▇▇█▇▆██▇▆▄▅▄▄▅▆▇██████▇██▇███▅▆▆▆ █
  1.6 μs       Histogram: log(frequency) by time      5.94 μs <
 Memory estimate: 7.94 KiB, allocs estimate: 1.

@benchmark det(tril($A))

BenchmarkTools.Trial: 2407 samples with 1 evaluation.
 Range (min … max):  1.105 ms …  13.096 ms  ┊ GC (min … max):  0.00% … 86.75%
 Time  (median):     1.807 ms               ┊ GC (median):     0.00%
 Time  (mean ± σ):   2.066 ms ± 949.353 μs  ┊ GC (mean ± σ):  19.92% ± 24.28%
    ▁     ▁  ▄█▇▃                                              
  ▄██▄▃▃▂▅█▇█████▇▅▄▃▃▂▂▂▂▁▂▁▁▁▁▁▁▁▂▁▂▁▁▁▁▂▂▃▄▆▅▅▃▃▃▃▂▂▂▂▂▂▂▂ ▃
  1.11 ms         Histogram: frequency by time        4.38 ms <
 Memory estimate: 7.64 MiB, allocs estimate: 3.

# if we tell Julia it's triangular: O(n) complexity
@benchmark det(LowerTriangular($A))

BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.757 μs …  1.524 ms  ┊ GC (min … max):  0.00% … 99.47%
 Time  (median):     2.406 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):   3.628 μs ± 31.077 μs  ┊ GC (mean ± σ):  31.64% ±  3.84%
         ▁▆█▆▂                                                
  ▃▄▃▃▄▅██████▆▃▃▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  1.76 μs        Histogram: frequency by time        5.77 μs <
 Memory estimate: 7.94 KiB, allocs estimate: 1.