# 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. Oct 17, 2015 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. Oct 17, 2015 at 2:22
• There is also the BLAS function LinearAlgebraBLASTRSV[], which might be used internally by LinearSolveFunction[]. Oct 17, 2015 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. Oct 18, 2015 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, 2015 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][]}, {n, 10}];
timeTri = Table[{100 n, RepeatedTiming[tritab[[n]] @ vecs[[n]], .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][]},
{n, 10}];

(* timing LinearSolve[] calls *)
timeTri = Table[{100 n, RepeatedTiming[tritab[[n]]@vecs[[n]], .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, 2015 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. Oct 20, 2015 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, 2015 at 9:41
• @thank_you Have you seen this? Jan 29, 2017 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][]}, {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][]}, {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? Oct 17, 2015 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, 2015 at 17:54
• What about if the matrices are sparse? Also, where do I get RepeatedTiming?
– fred
Oct 19, 2015 at 2:20
• @fred, RepeatedTiming` is a new function in version 10.1. Oct 20, 2015 at 3:19