# Does LinearSolve know that I gave it a Sparse lower triangular matrix?

Let $\mathbf L$ be a nonsingular numeric (non-symbolic) large, lower triangular sparse matrix. Let $\mathbf b$ be an appropriately sized vector. I want to solve $\mathbf L\mathbf x=\mathbf b$.

Does LinearSolve know that it can perform backsubstitution? There is an undocumented function called LUBackSubstitution. Wolfram's help website says its functionality has been superseded by LinearSolve. Is this true for sparse arrays as well?

Edit: All the responses so far have only addressed dense matrices, I was only interested in sparse matrices.

• Re: LUBackSubstitution[]: it used to be documented, but it was quietly shuffled off to the background after LinearSolveFunction[] (which does store the decomposed matrix in the same format internally) became available. I imagine LinearSolve[] is smart enough to detect triangular matrices. – J. M.'s technical difficulties Oct 17 '15 at 1:56
• Another note: as is well-known, Mathematica uses LAPACK behind the scenes. In relation to this, there is the function LinearAlgebraLAPACKLATRS[], which, though undocumented in Mathematica, is familiar to LAPACK users. The usage is a bit cumbersome, tho. – J. M.'s technical difficulties Oct 17 '15 at 2:22
• There is also the BLAS function LinearAlgebraBLASTRSV[], which might be used internally by LinearSolveFunction[]. – J. M.'s technical difficulties Oct 17 '15 at 23:14
• I tried using Oleksandr's trick here to see if I could intercept any calls from either of the two suspects I gave above. Unfortunately, LinearSolve[] is apparently opaque to this. – J. M.'s technical difficulties Oct 18 '15 at 6:23
• @rhermans While I am grateful for all the work people did in attempting to answer my question, neither of the responses uses a SparseArray. This is fundamental to my question. I should have put that in the title. Also this problem come from the discretization of a PDE so $1000 \times 1000$ would actually be on the "small" end of things. The dimension of the problem will be more like $100,000\times 100,000$. How should I address this? – fred Oct 23 '15 at 12:56

I decided to do the test I proposed in a comment myself. It is known that LinearSolve[], when applied to just a matrix, generates a LinearSolveFunction[] that internally stores the LU decomposition of a matrix. I feel that it would be more accurate for timing purposes to decouple the decomposition and backsubstitution phases, and take the timings only on the backsubstitution.

With that,

(* LU-decomposed matrices *)
lstab = Table[LinearSolve[RandomVariate[NormalDistribution[], {k, k}]],
{k, 100, 1000, 100}];

(* random vectors with matching dimensions *)
vecs = Table[RandomVariate[NormalDistribution[], k], {k, 100, 1000, 100}];

(* extract and transpose upper triangular factors for lower triangular example *)
tritab = Composition[LinearSolve, Transpose] /@ Through[lstab["getU"]];

(* ...and, timing! *)
timeFull = Table[{100 n, RepeatedTiming[lstab[[n]] @ vecs[[n]], .1][[1]]}, {n, 10}];
timeTri = Table[{100 n, RepeatedTiming[tritab[[n]] @ vecs[[n]], .1][[1]]}, {n, 10}];

ListLinePlot[{Legended[timeTri, "lower triangular"],
Legended[timeFull, "full matrix"]},
PlotMarkers -> Automatic, PlotRange -> All]


At the very least, we can guess that the backsubstitution behind the scenes is faster in the lower triangular case than in the full case, since the latter has to do two of them.

As a further comparison, I tried comparing the use of LinearSolve[] on a lower triangular matrix against a direct call to the internal BLAS function *TRSV:

(* lower triangular matrices *)
trimat = Transpose /@ Through[lstab["getU"]];

(* timing BLAS calls *)
timeBLAS = Table[{100 n, RepeatedTiming[
Block[{mat = trimat[[n]], x = vecs[[n]]},
LinearAlgebraBLASTRSV["L", "N", "N", mat, x]], .1][[1]]},
{n, 10}];

(* timing LinearSolve[] calls *)
timeTri = Table[{100 n, RepeatedTiming[tritab[[n]]@vecs[[n]], .1][[1]]}, {n, 10}];

ListLinePlot[{Legended[timeTri, "LinearSolve"],
Legended[timeBLAS, "BLAS"]},
PlotMarkers -> Automatic, PlotRange -> All]


If the overhead of LinearSolve[] is a concern, I suppose a direct BLAS call can be used instead.

• Although I can achieve many functions by ?LinearAlgebra* in the context LinearAlgebra, I cannot use these functions owning to that they are undocumented. So J.M. , could you tell me how did you know the usage of these function? thanks. – xyz Oct 20 '15 at 9:22
• I'm familiar with BLAS and LAPACK from using them in Fortran. It was merely a matter of matching what I knew about Fortran usage to how it might be called in Mathematica. – J. M.'s technical difficulties Oct 20 '15 at 9:34
• OK, @J.M. thanks for your explonation I never learned Fortran. I only learned the Visual Basic, C and Python. – xyz Oct 20 '15 at 9:41
• @thank_you Have you seen this? – Michael E2 Jan 29 '17 at 17:52

Weird. It seems that LinearSolve takes longer for a triangular matrix.

(* create a well-conditioned random lower triangular matrix *)
rndLowerTriangular[n_] :=
Normal[SparseArray[{{i_, j_} /; i > j -> .1, {i_, i_} -> 10}, {n,
n}]*RandomReal[{1, 2}, {n, n}]]

randomProblem[n_] :=
LinearSolve[rndLowerTriangular[n], RandomReal[{-1, 1}, n]]

Monitor[time =
Table[{n, RepeatedTiming[randomProblem[n], .1][[1]]}, {n, 100,
1000, 100}], n];

randomProblemFull[n_] :=
LinearSolve[RandomReal[{-1, 1}, {n, n}], RandomReal[{-1, 1}, n]]

Monitor[timeFull =
Table[{n, RepeatedTiming[randomProblemFull[n], .1][[1]]}, {n, 100,
1000, 100}], n];

ListPlot[{Legended[time, "Lower triangular"],
Legended[timeFull, "Full matrix"]}, PlotRange -> All, Mesh -> All,
Joined -> True]


• May I suggest another test? How about forming the LinearSolveFunction[] at once from the test matrix (to separate out the LU decomposition step), and just time the application of the function to a random vector? – J. M.'s technical difficulties Oct 17 '15 at 16:42
• I get pretty much identical timings if I make the matrices more comparable by defining rndFull[n_]:=(#+Transpose[#])&@rndLowerTriangular[n] and using that in the full matrix test. – Jens Oct 17 '15 at 17:54
• What about if the matrices are sparse? Also, where do I get RepeatedTiming? – fred Oct 19 '15 at 2:20
• @fred, RepeatedTiming` is a new function in version 10.1. – J. M.'s technical difficulties Oct 20 '15 at 3:19