Numerical Stability of Forward Substitution in Python

Question

I am implementing some basic linear equation solvers in Python.

I have currently implemented forward and backward substitution for triangular systems of equations (so very straightforward to solve!), but the precision of the solutions becomes very poor even with systems of about 50 equations (50x50 coefficient matrix).

The following code performs the forward/backward substitution:

FORWARD_SUBSTITUTION = 1
BACKWARD_SUBSTITUTION = 2
def solve_triang_subst(A: np.ndarray, b: np.ndarray,
                       substitution=FORWARD_SUBSTITUTION) -> np.ndarray:
    """Solves a triangular system via
    forward or backward substitution.

    A must be triangular. FORWARD_SUBSTITUTION means A should be
    lower-triangular, BACKWARD_SUBSTITUTION means A should be upper-triangular.
    """
    rows = len(A)
    x = np.zeros(rows, dtype=A.dtype)
    row_sequence = reversed(range(rows)) if substitution == BACKWARD_SUBSTITUTION else range(rows)

    for row in row_sequence:
        delta = b[row] - np.dot(A[row], x)
        cur_x = delta / A[row][row]
        x[row] = cur_x

    return x

I am using numpy and 64-bit floats.

Simple Testing Tool

I have set up a simple test suite which generates coefficient matrices and x vectors, computes the b, and then uses forward or backward substitution to recover the x, comparing it to the its known value for validity.

The following code performs these checks:

import numpy as np
import scipy.linalg as sp_la

RANDOM_SEED = 1984
np.random.seed(RANDOM_SEED)


def check(sol: np.ndarray, x_gt: np.ndarray, description: str) -> None:
    if not np.allclose(sol, x_gt, rtol=0.1):
        print("Found inaccurate solution:")
        print(sol)
        print("Ground truth (not achieved...):")
        print(x_gt)
        raise ValueError("{} did not work!".format(description))

def fuzz_test_solving():
    N_ITERATIONS = 100
    refine_result = True
    for mode in [FORWARD_SUBSTITUTION, BACKWARD_SUBSTITUTION]:
        print("Starting mode {}".format(mode))
        for iteration in range(N_ITERATIONS):
            N = np.random.randint(3, 50)
            A = np.random.uniform(0.0, 1.0, [N, N]).astype(np.float64)

            if mode == BACKWARD_SUBSTITUTION:
                A = np.triu(A)
            elif mode == FORWARD_SUBSTITUTION:
                A = np.tril(A)
            else:
                raise ValueError()

            x_gt = np.random.uniform(0.0, 1.0, N).astype(np.float64)
            b = np.dot(A, x_gt)

            x_est = solve_triang_subst(A, b, substitution=mode,
                                       refine_result=refine_result)
            # TODO report error and count, don't throw!
            # Keep track of error norm!!
            check(x_est, x_gt,
                  "Mode {} custom triang iteration {}".format(mode, iteration))

if __name__ == '__main__':
    fuzz_test_solving()

Note that the maximum size of a test matrix is 49x49. Even in this case, the system cannot always compute decent solutions, and fails by more than a margin of 0.1. Here's an example of such a failure (this is doing backward substitution, so the biggest error is in the 0th coefficient; all the test data are sampled uniformly from [0, 1[):

Solution found with Mode 2 custom triang iteration 24:
[ 0.27876067  0.55200497  0.49499509  0.3259397   0.62420183  0.47041149
  0.63557676  0.41155446  0.47191956  0.74385864  0.03002819  0.4700286
  0.37989592  0.56527691  0.15072607  0.05659282  0.52587574  0.82252197
  0.65662833  0.50250729  0.74139748  0.10852731  0.27864265  0.42981232
  0.16327331  0.74097937  0.24411709  0.96934199  0.890266    0.9183985
  0.14842446  0.51806495  0.36966843  0.18227989  0.85399593  0.89615663
  0.39819336  0.90445931  0.21430972  0.61212349  0.85205597  0.66758689
  0.1793689   0.38067267  0.39104614  0.6765885   0.4118123 ]
Ground truth (not achieved...)
[ 0.20881608  0.71009766  0.44735271  0.31169033  0.63982328  0.49075813
  0.59669585  0.43844108  0.47764942  0.72222069  0.03497499  0.4707452
  0.37679884  0.56439738  0.15120397  0.05635977  0.52616387  0.82230625
  0.65670245  0.50251426  0.74139956  0.10845974  0.27864289  0.42981226
  0.1632732   0.74097939  0.24411707  0.96934199  0.89026601  0.91839849
  0.14842446  0.51806495  0.36966843  0.18227989  0.85399593  0.89615663
  0.39819336  0.90445931  0.21430972  0.61212349  0.85205597  0.66758689
  0.1793689   0.38067267  0.39104614  0.6765885   0.4118123 ]

I have also implemented the iterative refinement method described in Section 2.5 of [0], and while it did help a little, the results are still poor for larger matrices.

MATLAB Sanity Check

I also did this experiment in MATLAB, and even there, once there are more than 100 equations, the estimation error shoots up exponentially.

Here is the MATLAB code I used for this experiment:

err_norms = [];
range = 1:3:120;
for size=range

  A = rand(size, size);
  A = tril(A);
  x_gt = rand(size, 1);

  b = A * x_gt;
  x_sol = A\b;

  err_norms = [err_norms, norm(x_gt - x_sol)];
end

plot(range, err_norms);
set(gca, 'YScale', 'log')

And here is the resulting plot:

Main Question

My question is: Is this normal behavior, seeing as there is essentially no structure in the problem, given that I randomly generate the A matrix and x?

What about solving linear systems of 100s of equations for various practical applications? Are these limitations simply an accepted fact, and e.g., optimization algorithms are just naturally robust to these issues? Or am I missing some important facets of this problem?

[0]: Press, William H. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.

The generated matrix is dense. I am aware that if systems start exhibiting some specific structure (e.g., tridiagonal, block-diagonal, whatever), then much larger systems can be solved, but in my case, I'm asking about the most general, fully dense case. — Andrei Bârsan, Nov 18 '17 at 01:12

percusse · Accepted Answer · 2017-11-18T03:41:15.640

There are no limitations. This is a very fruitful exercise that we all came to realize; writing linear solvers are not that easy and that's why almost always LAPACK or its cousins in other languages are used with full confidence.

You are hit by almost singular matrices and because you are using matlab's backslash you don't see that matlab is switching to least squares solutions behind the scenes when near singularity is hit. If you just change A\b to linsolve(A,b) hence you restrict the solver to solve square systems you'll probably see lots of warnings on your console.

I didn't test it because I don't have a license anymore but if I write blindly this should show you the condition numbers of the matrices at each step.

err_norms = [];
range = 1:3:120;
for i=1:40
  size = range(i);
  A = rand(size, size);
  A = tril(A);
  x_gt = rand(size, 1);

  b = A * x_gt;
  x_sol = linsolve(A,b);

  err_norms = [err_norms, norm(x_gt - x_sol)];
  zzz(i) = rcond(A);
end

semilogy(range, err_norms);
figure,semilogy(range,zzz);

Note that because you are picking up numbers from a uniform distribution it becomes more and more likely to hit ill-conditioned matrices (wrt to inversion) as the rows have more probability to have rank deficiency. That's why the error becomes bigger and bigger. Sprinkle some identity matrix times a scalar and all errors should come back to eps*n levels.

But best, leave this to expert algorithms which have been tested through decades. It is really not that trivial to write any of these. You can read the Fortran codes, for example, dtrsm solves the triangular system.

On the Python side, you can use scipy.linalg.solve_triangular which uses ?trtrs routines from LAPACK.

Thanks for the detailed answer! This makes sense. I actually did check `?trtrs` in LAPACK (starting exactly with the method you pointed out) and, `CTRTRS` in BLAS (caled by `?trtrs`), and I didn't notice whether they were doing anything specific to address these issues (e.g., it doesn't seem to attempt to refine the solution, but at the same time, maybe I missed it since I'm not too experienced in reading Fortran). — Andrei Bârsan, Nov 18 '17 at 17:48
Yes, everything works as expected when "nudging" the diagonal away from zero, i.e., making the matrix less singular. Oh, and of course, in practical settings I will still use existing implementations of numerical algorithms. I was implementing these things for educational purposes. Thank you once again! — Andrei Bârsan, Nov 18 '17 at 17:58

Numerical Stability of Forward Substitution in Python

Simple Testing Tool

MATLAB Sanity Check

Main Question

1 Answers1