6
$\begingroup$

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?

Thanks in advance,

Improve this question
$\endgroup$
7
  • $\begingroup$ I flagged this for migration to SO $\endgroup$
    – Sabyasachi
    Mar 5, 2014 at 10:08
  • $\begingroup$ I'm not an expert, but I think this question will be better off on stackoverflow... $\endgroup$
    – 5xum
    Mar 5, 2014 at 10:08
  • $\begingroup$ Yeah I did flag it. $\endgroup$
    – Sabyasachi
    Mar 5, 2014 at 10:10
  • 1
    $\begingroup$ 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. $\endgroup$
    – Sabyasachi
    Mar 5, 2014 at 10:11
  • 2
    $\begingroup$ @Sabyasachi I am migrating to scicomp instead. Questions about SciPy and NumPy are explicitly on topic there. $\endgroup$
    – Willie Wong
    Mar 5, 2014 at 11:19

3 Answers 3

Reset to default
9
$\begingroup$

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).

Improve this answer
$\endgroup$
1
  • $\begingroup$ Excellent. solve_triangular(A,b,lower=True) did the trick. Thanks. $\endgroup$
    – BBSysDyn
    Mar 7, 2014 at 22:06
7
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – sebas
    Mar 5, 2014 at 20:31
  • $\begingroup$ Don't worry about it. It happens. $\endgroup$
    – Geoff Oxberry
    Mar 5, 2014 at 20:37
  • $\begingroup$ 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)$. $\endgroup$
    – Federico Poloni
    Mar 7, 2014 at 13:16
  • $\begingroup$ It does; the link discusses that in more detail. $\endgroup$
    – Geoff Oxberry
    Mar 7, 2014 at 18:07
2
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – beldaz
    Mar 14, 2018 at 5:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.