Numerical Stability of Forward Substitution in Python

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);endsemilogy(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.