versioninfo()
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.
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).
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).
using LinearAlgebra, Random
Random.seed!(123) # seed
n = 5
A = randn(n, n)
b = randn(n)
# a random matrix
A
Al = LowerTriangular(A) # does not create extra matrix
dump(Al)
Al.data
# same data
pointer(Al.data), pointer(A)
Al \ b # dispatched to BLAS function for triangular solve
# or use BLAS wrapper directly
BLAS.trsv('L', 'N', 'N', A, b)
?BLAS.trsv
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.
LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular.
A
LowerTriangular(A)
LinearAlgebra.UnitLowerTriangular(A)
using BenchmarkTools, LinearAlgebra, Random
Random.seed!(123) # seed
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)))
# if we tell Julia it's triangular: O(n) complexity
@benchmark eigvals($(LowerTriangular(A)))
# 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)))
# if we tell Julia it's triangular: O(n) complexity
@benchmark det($(LowerTriangular(A)))
@benchmark det($(LowerTriangular(A)))