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{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 (triangular solve with one right hand side).

  • BLAS level 3 function: ?trsm (matrix triangular solve, i.e., multiple right hand sides).

  • Julia

    • 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!, trsv, trsm, trsm!
In [1]:
srand(123) # seed
n = 5
A = randn(n, n)
b = randn(n)

tril(A) # create another triangular matrix
Out[1]:
5×5 Array{Float64,2}:
  1.19027    0.0         0.0        0.0      0.0     
  2.04818    0.980968    0.0        0.0      0.0     
  1.14265   -0.0754831  -0.888936   0.0      0.0     
  0.459416   0.273815    0.327215  -0.71741  0.0     
 -0.396679  -0.194229    0.592403  -0.77507  0.277726
In [2]:
LowerTriangular(A) # does not create extra matrix
Out[2]:
5×5 LowerTriangular{Float64,Array{Float64,2}}:
  1.19027     ⋅           ⋅          ⋅        ⋅      
  2.04818    0.980968     ⋅          ⋅        ⋅      
  1.14265   -0.0754831  -0.888936    ⋅        ⋅      
  0.459416   0.273815    0.327215  -0.71741   ⋅      
 -0.396679  -0.194229    0.592403  -0.77507  0.277726
In [ ]:
tril(A) \ b # dispatched to LowerTriangular(A) \ b
In [24]:
LowerTriangular(A) \ b # dispatched to A_ldiv_B
Out[24]:
5-element Array{Float64,1}:
  1.28031
 -4.48541
  5.32613
  0.44682
 -3.09169
In [25]:
# or use BLAS wrapper directly
Base.LinAlg.BLAS.trsv('L', 'N', 'N', A, b)
Out[25]:
5-element Array{Float64,1}:
  1.28031
 -4.48541
  5.32613
  0.44682
 -3.09169

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.