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.
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).
srand(123) # seed
n = 5
A = randn(n, n)
b = randn(n)
# a random matrix
A
tril(A) # create another triangular matrix
Al = LowerTriangular(A) # does not create extra matrix
fieldnames(Al)
Al.data
tril(A) \ b # dispatched to LowerTriangular(A) \ b
# check source code
@which tril(A) \ b
LowerTriangular(A) \ b # dispatched to A_ldiv_B
# check source code
@which LowerTriangular(A) \ b
# or use BLAS wrapper directly
LinAlg.BLAS.trsv('L', 'N', 'N', A, b)
?LinAlg.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.
using BenchmarkTools
srand(123) # seed
A = randn(1000, 1000)
# if we don't tell Julia it's triangular
@benchmark eigvals(tril(A))
# if we tell Julia it's triangular
@benchmark eigvals(LowerTriangular(A))
# if we don't tell Julia it's triangular
@benchmark det(tril(A))
# if we tell Julia it's triangular
@benchmark det(LowerTriangular(A))
versioninfo()