7

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,

BBSysDyn
  • 239
  • 2
  • 6
  • I flagged this for migration to SO –  Mar 05 '14 at 10:08
  • I'm not an expert, but I think this question will be better off on stackoverflow... –  Mar 05 '14 at 10:08
  • Yeah I did flag it. –  Mar 05 '14 at 10:10
  • 1
    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. –  Mar 05 '14 at 10:11
  • 2
    @Sabyasachi I am migrating to scicomp instead. Questions about SciPy and NumPy are explicitly on topic there. – Willie Wong Mar 05 '14 at 11:19
  • Is your matrix A sparse? It is the whay to know if the matrix is lower triangular without check all the upper triangular entries of the matrix. Otherwise, do you know if the matrix is triangular before to perform the solver? Then you can parse the solver with a flag to perform triangular solver in that case. – sebas Mar 05 '14 at 11:37
  • Let's say A is not sparse, and lower triangular so above the diagonal is all zero. My question is would numpy.linalg.solve detect this and use backsubstitution to solve the equation directly? – BBSysDyn Mar 05 '14 at 11:54

3 Answers3

10

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

sebas
  • 429
  • 2
  • 8
7

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.

Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
2

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.

BBSysDyn
  • 239
  • 2
  • 6
  • 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 05:53