# Can numpy.linalg.solve use back substitution when possible?

The question is if Python Numpy library can use back subsitution to solve Ax=b if possible, that is, if A is lower triangular? Do numerical linear algebra packages do this? I would think Numpy would detect the triangular state and use the proper approach, but a Google search returns things like scipy.linalg.lu_solve or scipy.linalg.cho_solve, which I assume are to be used in case when we know we have a triangular matrix?

• I flagged this for migration to SO
– Sabyasachi
Mar 5, 2014 at 10:08
• I'm not an expert, but I think this question will be better off on stackoverflow...
– 5xum
Mar 5, 2014 at 10:08
• Yeah I did flag it.
– Sabyasachi
Mar 5, 2014 at 10:10
• Also I am not that well versed in Python, but I would assume that Numpy would call the proper scipy function if it detects that such an optimization would be beneficial. Pure speculation though.
– Sabyasachi
Mar 5, 2014 at 10:11
• @Sabyasachi I am migrating to scicomp instead. Questions about SciPy and NumPy are explicitly on topic there. Mar 5, 2014 at 11:19

Looking at the information of nympy.linalg.solve for dense matrices, it seems that they are calling LAPACK subroutine gesv, which perform the LU factorization of your matrix (without checking if the matrix is already lower triangular) and then solves the system. So the answer is NO.

Otherwise, it makes sense. If you do not have an easy (cheap) way to verify that your matrix is triangular (lower or upper), it can be very expensive to check this things.

If you know that your matrix is lower triangular, it is better to solve it in scipy with solve_triangular, while the matrix is still dense square matrix (so you are consumming a lot of effort).

• Excellent. solve_triangular(A,b,lower=True) did the trick. Thanks. Mar 7, 2014 at 22:06

No. The numpy.linalg.solve method uses LAPACK's DGESV, which is a general linear equation solver driver. If you know that your matrix is triangular, you should use a driver specialized for that matrix structure.

scipy.linalg.solve does something similar.

MATLAB detects triangularity in a solve if you use the backslash operator; see this page for pseudocode.

• Sorry @Geoff to duplicate the answer, I was writting my answer at the same time as you answered without realizing that. I agree completely with your answer. Mar 5, 2014 at 20:31
• Don't worry about it. It happens. Mar 5, 2014 at 20:37
• I think that Matlab also detects permuted triangularity, so that in [L,U]=lu(A);b=U\(L\b) the second instruction runs in $O(n^2)$. Mar 7, 2014 at 13:16
• It does; the link discusses that in more detail. Mar 7, 2014 at 18:07

I needed this for $y^T\Sigma^{-1}y$ calculation which can be solved by

$$y^T\Sigma^{-1}y= y^TL^{-T}L^{-1}y$$

where

$$\Sigma = LL^T$$

$$= (y^TL^{-T})L^{-1}y$$

$$= (L^{-1}y)^TL^{-1}y$$

$$= |L^{-1}y|^2$$

So we need $L^{-1}y$ which is the linear solution of $Lx=y$. Based on the two previous answers, I used (first vanilla solution)

import numpy.linalg as lin
Sigma = np.array([[10., 2.],[2., 5.]])
y = np.array([[1.],[2.]])
print np.dot(np.dot(y.T,lin.inv(Sigma)),y)


Now with solve_triangular

import scipy.linalg as slin
L = lin.cholesky(Sigma)
x = slin.solve_triangular(L,y,lower=True)
print np.dot(x.T,x)


Both give 0.804.

• This would be a valuable answer if you were to give a comparison for a large matrix. For a 2x2 example the costs are likely dominated by basic setup/housekeeping operations, which explains why there's no difference. Mar 14, 2018 at 5:53