---
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)))
```
