Issue with scipy.linalg.solve when using Random Lower Triangular matrix

[akihis2]

This is question about odd behaviour I have encountered while using Numpy.

I am observing very odd behaviour with 'scipy.linalg.solve' when I am using it to solve randomly generated lower triangular matrix.

I solve the randomly generated lower triangular Numpy 2D-array with randomly generated right-hand-side using 'scipy.linalg.solve' and 'scipy.linalg.solve_triangular' and compare their solution by calculating the norm of the difference between them. For small matrix, the error is negligible, but as I make the matrix larger (100x100), the error becomes very significant (I have seen it up to 10^23). I am setting the 'lower' flag to 'True' in both case to tell the solver that the input matrix is lower triangular matrix.

I have also exported the randomly generated matrix and the right-hand-side to Matlab to see which solver seems to output "correct" solution according to Matlab. The Matlab solution seems to agree with the solution from 'scipy.linalg.solve_triangular'. I have also hand calculated first few solution using back-substitution, and again, the solution seems to agree with result from 'scipy.linalg.solve_triangular'.

This is very odd behaviour which only seems to be encountered when you use lower triangular matrix (which should be easy to solve by the solver). If instead, I use upper triangular matrix, the error calculated between the two solver is always equal to zero (for the cases I have tested).

Anyone know what might be causing this odd behaviour with 'scipy.linalg.solve'? I have always thought that the 'scipy.linalg.solve' is direct substitute for Matlab 'mldivide ()', but it seems to have more limitation than 'mldivide' that I was not aware of, until I tested this.

I have attached the toy code I was using to test this behaviour, in case I just had stupid programming mistakes.

import numpy as np
import scipy.linalg

arr_size = 100;

a = np.array(np.random.randint(1,arr_size,(arr_size,arr_size)), dtype=np.double)
a = np.tril(a)

b = np.array(np.random.randint(1,arr_size,arr_size),  dtype=np.double)

x_solve = scipy.linalg.solve(a,b,lower=True)
x_trian = scipy.linalg.solve_triangular(a,b,lower=True)

error = np.linalg.norm(x_solve - x_trian)
print error

[quantheory]

Large randomly generated triangular matrices are often ill-conditioned, so that a small change in b (or the algorithm used) leads to a very large change in the calculated x. Without knowing too much about the specific details of solve, I'm betting that it uses a method which is inferior for ill-conditioned lower-triangular matrices? In fact, I bet that if your starting matrices are random floats (not integers), the solve_triangular method will also do poorly.