# 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`

.

`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#

`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

\(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#

\(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:

\(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\).

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#

`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#

\(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#

\(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

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

\(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#

\(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#

\(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)\).

*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\)#

\(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\)#

\(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`