# 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 '14 at 10:08
• I'm not an expert, but I think this question will be better off on stackoverflow... – 5xum Mar 5 '14 at 10:08
• Yeah I did flag it. – Sabyasachi Mar 5 '14 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 '14 at 10:11
• @Sabyasachi I am migrating to scicomp instead. Questions about SciPy and NumPy are explicitly on topic there. – Willie Wong Mar 5 '14 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. – BBSysDyn Mar 7 '14 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. – sebas Mar 5 '14 at 20:31
• Don't worry about it. It happens. – Geoff Oxberry Mar 5 '14 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)$. – Federico Poloni Mar 7 '14 at 13:16
• It does; the link discusses that in more detail. – Geoff Oxberry Mar 7 '14 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. – beldaz Mar 14 '18 at 5:53