Linear algebra algorithms using 0-based indexing and semi-open intervals¶

This section describes some of the core algorithms of linear algebra using the same indexing conventions as in most modern programming languages: Python, C, Java, C++, javascript, Objective-C, C#, Swift, etc. (In fact, almost everything except Matlab and Fortran.)

The key elements of this are:

Indices for vectors and other arrays start at 0.
Ranges of indices are described with semi-open intervals $[a,b)$.
This “index interval” notation has two virtues: it emphasizes the mathematical fact that the order in which things are done is irrelevant (such as within sums), and it more closely resembles the way that most programming languages specify index ranges. For example, the indices $i$ of a Python array with $n$ elements are $0 \leq i < n$, or $[0,n)$, and the Python notations range(n), range(0,n), :n and 0:n all describe this. Similarly, in Java, C C++ etc., one can loop over the indices $i \in [a,b)$ with for(i=a, i<b, i+=1)

The one place that the indexing is still a bit tricky is counting backwards!
For this, note that the index range $i = b, b-1, \dots a$ is $b \geq i > a-1$, which in Python is range(b, a-1, -1).

I include Python code for comparison for just the three most basic algorithms: “naive” LU factorization and forward and backward substitution, without pivoting. The rest are good exercises for learning how to program loops and sums.

The naive Gaussian elimination algorithm¶

In this careful version, the original matrix $A$ is called $A^{(0)}$, and the new versions at each stage are called $A^{(1)}$, $A^{(2)}$, and so on to $A^{(n-1)}$, which is the row-reduced form also called $U$; likewise with the right-hand sides $b^{(0)}=b$, $b^{(1)}$ up to $b^{(n-1)} = c$.

However, in software all those super-scripts can be ignored, just updating arrays A anf b.

for k in [0, n-1)
$\quad$ for i in [k+1, n)
$\quad\quad l_{i,k} = a_{i,k}^{(k)}/a_{k,k}^{(k)}$
$\quad\quad$ for j in [k+1, n)
$\quad\quad\quad a_{i,j}^{(k+1)} = a_{i,j}^{(k)} - l_{i,k} a_{k,j}^{(k)}$
$\quad\quad$ end for
$\quad\quad b_i^{(k+1)} = b_i^{(k)} - l_{i,k} b_{k}^{(k)}$
$\quad$ end for
end for

Actually this skips formulas for some elements of the new matrix $A^{(k+1)}$, because they are either zero or are unchanged from $A^{(k)}$:
the rows before $i=k$ are unchanged, and for columns before $j=k$, the new entries are zeros.

The LU factorization with $L$ unit lower triangular, by Doolittle’s direct method¶

Reference: the textbook’s section Matrix Factorizations.

for j in [0, n)
$\quad$ $u_{0,j} = a_{0,j}$
end for
for i in [1, n)
$\quad$ $l_{i,0} = a_{i,0}/u_{0,0}$
end for

for k in [1, n)
$\quad$ for j in [k, n)
$\quad\quad \displaystyle u_{k,j} = a_{k,j} - \sum_{s \in [0, k)} l_{k,s}u_{s,j}$
$\quad$ end for
$\quad$ for i in [k+1, n)
$\quad\quad \displaystyle l_{i,k} = \left( a_{i,k} - \sum_{s \in [0, k)} l_{i,s}u_{s,k} \right)/u_{k,k}$
$\quad$ end for
end for

Forward substitution with a unit lower triangular matrix¶

For $L$ unit lower triangular, solving $Lc = b$ by forward substitution is

$c_0 = b_0$
for i in [1, n)
$\quad \displaystyle c_i = b_i - \sum_{j \in [0,i)} l_{i,j} c_j$
end for

(Update: redundant division by $l_{0,0}$ removed.)

Backward substitution with an upper triangular matrix¶

$x_{n-1} = c_{n-1}/u_{n-1,n-1}$
for i from n-2 down to 0
$\displaystyle \quad x_i = \frac{c_i - \sum_{j \in [i+1, n)} u_{i,j} x_j}{u_{i,i}}$
end for

Counting backwards in Python¶

For the Python implementation, we need a range of indices that change by a step of -1 instead of 1. This can be done with an optional third argument to range: range(a, b, step) generates a succession values starting with a, each value differing from its predecessor by step, and stopping just before b. That last rule requires care whe the step is negative: for example, range(3, 0, -1) gives the sequence $\{3, 2, 1\}$. So to count down to zero, one has to use $b = -1$! That is, to count down from $n-1$ and end at 0, one uses range(n-1, -1, -1).

import numpy as np

def solveUpperTriangular(U, c):
    n = len(c)
    x = np.zeros(n)
    x[n-1] = c[n-1]/U[n-1,n-1]
    for i in range(n-2, -1, -1):
        x[i] = ( c[i] - sum(U[i,i+1:n] * x[i+1:n]) )/U[i,i]
    return x

Another way to count downwards¶

A range of integers $i$ from $b$ down to $a$ is also $i = b-d$ for $d$ from 0 up to $b-a$,
so one can use for d in [0, b-a)

This gives the alternative for the algorithm:

$x_{n-1} = c_{n-1}/u_{n-1,n-1}$
for d in [2, n+1):
$\quad i = n - d$
$\displaystyle \quad x_i = \frac{c_i - \sum_{j \in [i+1, n)} u_{i,j} x_j}{u_{i,i}}$
end for

This corresponds to the alternative Python code:

def solveUpperTriangular(U, c):
    n = len(c)
    x = np.zeros(n)
    x[n-1] = c[n-1]/U[n-1,n-1]
    for nmi in range(2, n+1):  # nmi is n-i
        i = n - nmi
        x[i] = ( c[i] - sum(U[i,i+1:n] * x[i+1:n]) )/U[i,i]
    return x

Versions with maximal element partial pivoting¶

Apart the choices of pivot rows and updating of the permutation vector $p$, the only change from the non-pivoting version is that all row indices change from $i$ to $p_i$ and so on, in both $U$ and $c$; column indices are unchanged.

Gaussian elimination with maximal element partial pivoting¶

In the following description, I will discard the above distinction between the successive matrices $A^{(k)}$ and vectors $b^{(k)}$, and instead refer to $A$ and $b$ like variable arrays in a programming language, with their elements being updated. Likewise, the permutation will be stored in a variable array $p$.

Initialize the permuation vector as $p = [0, 1, \dots, n-1]$
for k in [0, n-1)
$\quad$ Search elements $a_{p_i,k}$ for $i \in [k, n)$ and find the index r of the one with largest absolute value.
$\quad$ If $r\neq k$, swap $p_k$ with $p_r$
$\quad$ for i in [k+1, n)
$\quad\quad l_{p_i,k} = a_{p_i,k}/a_{p_k,k}$
$\quad\quad$ for j in [k+1, n)
$\quad\quad\quad a_{p_i,j} = a_{p_i,j} - l_{p_i,k} a_{p_k,j}$
$\quad\quad$ end for
$\quad\quad b_{p_i} = b_{p_i} - l_{p_i,k} b_{p_k}$
$\quad$ end for
end for

The Doolittle LU factorization algorithm with maximal element partial pivoting¶

for k in [0, n)
$\quad$ Search elements $a_{p_i,k}$ for $i \in [k, n)$ and find the index r of the one with largest absolute value.
$\quad$ If $r\neq k$, swap $p_k$ with $p_r$
$\quad$ for j in [k, n)
$\quad\quad$ Note that for $k=0$, the sum can [and should!] be omitted in the following line:
$\quad\quad \displaystyle u_{p_k,j} = a_{p_k,j} - \sum_{s \in [0, k)} l_{p_k,s}u_{p_s,j}$
$\quad$ end for
$\quad$ for i in [k+1, n)
$\quad\quad$ Note that for $k=0$, the sum can [and should!] be omitted in the following line:
$\quad\quad \displaystyle l_{p_i,k} = \left( a_{p_i,k} - \sum_{s \in [0, k)} l_{p_i,s}u_{p_s,k} \right)/u_{p_k,k}$
$\quad$ end for
end for

Forward substitution with maximal element partial pivoting¶

$c_{p_0} = b_{p_0}/l_{p_0,0}$
for i in [1, n)
$\quad \displaystyle c_{p_i} = b_{p_i} - \sum_{j \in [0, i)} l_{p_i,j} c_{p_j}$
end for

Backward substitution with maximal element partial pivoting¶

$x_{n-1} = c_{p_{n-1}}/u_{p_{n-1},n-1}$
for i from n-2 down to 0
$\displaystyle \quad x_i = \frac{c_{p_i} - \displaystyle \sum_{j \in [i+1, n)} u_{p_i,j} x_j}{u_{p_i,i}}$
end for

Tridiagonal matrix algorithms¶

Describe a tridiagonal matrix with three 1D arrays as $$ T = \left[\begin{array}{cc} d_0 & u_0\\ l_0 & d_1 & u_1 & \\ & l_1 & d_2 & u_2 \\ && \ddots & \ddots & \ddots \\ &&& l_{n-3} & d_{n-2} & u_{n-2} \\ &&&& l_{n-2} & d_{n-1} \end{array}\right] $$ with all "missing" entries being zeros, and the right had side of system $T x = b$ as $$\left[\begin{array}{c} b_0 \\ b_2 \\ \vdots \\ b_{n-1} \end{array}\right]$$

Doolittle factorization of a tridiagonal matrix¶

The factorization has the form $ T = L U$ with $$ L = \left[\begin{array}{cc} 1\\ L_0 & 1 \\ & L_1 & 1 \\ && \ddots & \ddots \\ &&& L_{n-3} & 1 \\ &&&& L_{n-2} & 1 \end{array}\right], \; U = \left[\begin{array}{cc} D_0 & u_0\\ & D_1 & u_1 \\ && D_2 & u_2\\ &&& \ddots & \ddots \\ &&&& D_{n-2} & u_{n-2} \\ &&&&& D_{n-1} \end{array}\right] $$ so just the arrays $L$ and $D$ are to be computed.

$D_0 = d_0$
for i in [1, n)
$\quad L_{i-1} = l_{i-1}/D_{i-1}$
$\quad D_{i} = d_{i} - L_{i-1} u_{i-1}$
end for

Forward substitution with a tridiagonal matrix¶

$c_0 = b_0$
for i in [1, n)
$\quad c_{i} = b_{i} - L_{i-1} c_{i-1}$
end for

Backward substitution with a tridiagonal matrix¶

$x_{n-1} = c_{n-1}/D_{n-1}$
for i from n-2 down to 0
$\quad x_{i} = (c_{i} - u_{i} x_{i+1})/D_i$
end for

Banded matrix algorithms¶

Doolittle factorization of a matrix of bandwidth $p$¶

That is, $A_{i,j} = 0$ when $|i-j| > p$.

In addito to loops stopping at the point beyond which valeus would be zero or unchanged the top row and top-right diagonal overlapping at the element of “1-based” indices $(1,p+1)$, which is now at 0-based indices $(0,p)$.

The top row is unchanged:
for j in [0, p+1)
$\quad$ $u_{0,j} = a_{0,j}$
end for
The top non-zero diagonal is also unchanged:
for k in [1, n - p)
$\quad$ $u_{k,k+p} = a_{k,k+p}$
end for
The left column requires no sums:
for i in [1, p+1)
$\quad$ $l_{i,0} = a_{i,0}/u_{0,0}$
end for

The main loop:
for k in [1, n)
$\quad$ for j in [k, min(n, k+p))
$\quad\quad$ $u_{k,j} = a_{k,j} - \displaystyle\sum_{s \in [\max(0, j-p), k)} l_{k,s} u_{s,j}$
$\quad$ end for
$\quad$ for i in [k+1, min(n,k+p+1))
$\quad\quad$ $l_{i,k} = \left( a_{i,k} - \displaystyle \sum_{s \in [\max(0,i-p), k)} l_{i,s} u_{s,k} \right)/u_{k,k}$
$\quad$ end for
end for

Forward substitution with a unit lower-triangular matrix of bandwidth $p$¶

$c_0 = b_0/l_{0,0}$
for i in [1, n)
$\quad \displaystyle c_i = b_i - \sum_{j \in [max(0,i-p),i)} l_{i,j} c_j$
end for

Backward substitution with an upper-triangular matrix of bandwidth $p$¶

$x_{n-1} = c_{n-1}/u_{n-1,n-1}$
for i from n-2 down to 0
$\displaystyle \quad x_i = \frac{c_i - \sum_{j \in [i+1, \min(n, i+p+1))} u_{i,j} x_j}{u_{i,i}}$
end for

This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International

MATH 375. Elementary Numerical Analysis (with Python)