Back-substitution

In [2]: 
import numpy as np
In [11]: 
n = 5

A = np.random.randn(n, n) * np.tri(n).T
print(A)

x = np.random.randn(n)
print(x)

b = np.dot(A, x)

[[-1.26236737 -0.8644523   1.55110419 -0.94165954 -0.71166821]
 [-0.         -1.89991829 -1.12215066  0.16162471 -0.5094088 ]
 [-0.         -0.         -0.52611369  1.03649351 -1.03046035]
 [-0.         -0.          0.          0.22869562 -0.45786146]
 [-0.         -0.          0.         -0.          0.19889282]]
[ 1.35615426 -0.7539793  -0.04295377  0.12033124 -1.9996183 ]
In [16]: 
xcomp = np.zeros(n)

for i in range(n-1, -1, -1):
    tmp = b[i]
    for j in range(n-1, i, -1):
        tmp -= xcomp[j]*A[i,j]
        
    xcomp[i] = tmp/A[i,i]

Now compare the computed $x$ against the reference solution.

In [19]: 
print(x)
print(xcomp)

[ 1.35615426 -0.7539793  -0.04295377  0.12033124 -1.9996183 ]
[ 1.35615426 -0.7539793  -0.04295377  0.12033124 -1.9996183 ]

Questions/comments:

  • Can this fail?
  • What's the operation count?