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

2. 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 and b.

Algorithm 2.19

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
$\quad\quad b_i^{(k+1)} = b_i^{(k)} - l_{i,k} b_{k}^{(k)}$
$\quad$ end
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#

Algorithm 2.20

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

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
$\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
end

Forward substitution with a unit lower triangular matrix#

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

Algorithm 2.21

$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

Backward substitution with an upper triangular matrix#

Algorithm 2.22

$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

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:

Algorithm 2.23

$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

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$.

Algorithm 2.24

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
$\quad\quad b_{p_i} = b_{p_i} - l_{p_i,k} b_{p_k}$
$\quad$ end
end

The Doolittle LU factorization algorithm with maximal element partial pivoting#

Algorithm 2.25

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
$\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
end

Forward substitution with maximal element partial pivoting#

Algorithm 2.26

$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

Backward substitution with maximal element partial pivoting#

Algorithm 2.27

$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

Tridiagonal matrix algorithms#

Describe a tridiagonal matrix with three 1D arrays as $$ T = \left[\begin{array}{cccccc} 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

\[\begin{split}\left[\begin{array}{c} b_0 \\ b_2 \\ \vdots \\ b_{n-1} \end{array}\right]\end{split}\]

Doolittle factorization of a tridiagonal matrix#

The factorization has the form $ T = L U$ with $$ L = \left[\begin{array}{cccccc} 1\\ L_0 & 1 \\ & L_1 & 1 \\ && \ddots & \ddots \\ &&& L_{n-3} & 1 \\ &&&& L_{n-2} & 1 \end{array}\right], \; U = \left[\begin{array}{cccccc} 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.

Algorithm 2.28

$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

Forward substitution with a tridiagonal matrix#

Algorithm 2.29

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

Backward substitution with a tridiagonal matrix#

Algorithm 2.30

$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

Banded matrix algorithms#

Doolittle factorization of a matrix of bandwidth $p$#

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

In addition 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)$.

Algorithm 2.31

The top row is unchanged:
for j in [0, p+1)
$\quad$ $u_{0,j} = a_{0,j}$
end
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
$\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
end

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

Algorithm 2.32

$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

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

Algorithm 2.33

$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

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

Contents

2. Linear algebra algorithms using 0-based indexing and semi-open intervals#

The naive Gaussian elimination algorithm#

The LU factorization with \(L\) unit lower triangular, by Doolittle’s direct method#