Triangular systems

Contents

\(\newcommand{\re}{\mathbb{R}}\)

Triangular systems#

Given \(L \in \re^{m \times m}\) lower-triangular, non-singular matrix, we want to solve

\[ L x = b \]

The \(i\)’th equation is

\[ \sum_{j=0}^{m-1} L_{ij} x_j = b_i, \qquad 0 \le i \le m-1 \]

and since \(L\) is lower-triangular, this simplifies to

\[ \sum_{j=0}^{i} L_{ij} x_j = b_i \]

which can be solved for \(x_i\)

\[ x_i = \frac{1}{L_{ii}} \left[b_i - \sum_{j=0}^{i-1} L_{ij} x_j\right] \]

We can solve using forward-substitution, i.e., solve for \(x_0\), then \(x_1\), and so on.

Forward-substitution: For \(i=0,1,2,\ldots,m-1\)

\[ x_i = \frac{1}{L_{ii}} \left[b_i - \sum_{j=0}^{i-1} L_{ij} x_j\right] \]

Since \(L\) is non-singular, the diagonal entries are non-zero, \(L_{ii} \ne 0\).

import numpy as np

We can compute the sum

\[ \sum_{j=0}^{i-1} L_{ij} x_j \]

using a for loop, but this is slow in Python. Instead, we can compute it as a dot product

\[ \sum_{j=0}^{i-1} L_{ij} x_j = \texttt{np.dot(L[i,0:i], x[0:i])} \]
def LSolve(L,b):
    m = L.shape[0]
    # solve Lx = b
    x = np.empty_like(b)
    for i in range(m):
        x[i] = (b[i] - np.dot(L[i,0:i], x[0:i])) / L[i,i]
    return x

Create a random matrix and right hand side vector \(b\) and solve \(L x = b\).

m = 10
A = np.random.rand(m,m) # random matrix
L = np.tril(A)          # lower-triangular part
b = np.random.rand(m)
x = LSolve(L,b)
print('||Lx-b|| = ',np.linalg.norm(L@x - b,np.inf))

||Lx-b|| =  8.881784197001252e-16

Exercise#

Write a code to solve an upper-triangular, non-singular system and test it.

Contents