Solving an upper triangular system by backsubstitution

# Solves Ax = b by backsubstitution, where the matrix A is upper triangular (with zeros below the diagonal) and b is a column vector with the same number of rows as A. The idea is to use the matrix structure for more efficient computation. See "Numerical methods in matrix computations", Åke Björck, 2014. def solve_triangular(A, b): n = len(b) x = [0] * n x[n-1] = b[n-1] / float(A[n-1][n-1]) for i in range(n-2, -1, -1): # (start, stop, step) s = sum([A[i][k] * x[k] for k in range(i + 1, n, 1)]) x[i] = (b[i] - s) / A[i][i] return x A = [ [1,5,4,3,6], [0,3,2,4,1], [0,0,5,3,1], [0,0,0,2,4], [0,0,0,0,5] ] b = [4,3,5,6,3] print(solve_triangular(A, b)) # [3.133333333333333, -1.4666666666666666, -0.2, 1.8, 0.6] # Check answer from scipy.linalg import solve_triangular print(solve_triangular(A, b)) # [3.13333333, -1.46666667, -0.2, 1.8, 0.6]