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.

  • $\begingroup$ 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. $\endgroup$ – J. M.'s ennui Oct 17 '15 at 1:56
  • $\begingroup$ Another note: as is well-known, Mathematica uses LAPACK behind the scenes. In relation to this, there is the function LinearAlgebra`LAPACK`LATRS[], which, though undocumented in Mathematica, is familiar to LAPACK users. The usage is a bit cumbersome, tho. $\endgroup$ – J. M.'s ennui Oct 17 '15 at 2:22
  • $\begingroup$ There is also the BLAS function LinearAlgebra`BLAS`TRSV[], which might be used internally by LinearSolveFunction[]. $\endgroup$ – J. M.'s ennui Oct 17 '15 at 23:14
  • $\begingroup$ 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. $\endgroup$ – J. M.'s ennui Oct 18 '15 at 6:23
  • $\begingroup$ @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? $\endgroup$ – 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]

lower triangular versus full solves

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]]}, 
                       LinearAlgebra`BLAS`TRSV["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]

LinearSolve[] versus BLAS

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

  • 1
    $\begingroup$ 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. $\endgroup$ – xyz Oct 20 '15 at 9:22
  • $\begingroup$ 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. $\endgroup$ – J. M.'s ennui Oct 20 '15 at 9:34
  • $\begingroup$ OK, @J.M. thanks for your explonation I never learned Fortran. I only learned the Visual Basic, C and Python. $\endgroup$ – xyz Oct 20 '15 at 9:41
  • $\begingroup$ @thank_you Have you seen this? $\endgroup$ – 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]

enter image description here

  • $\begingroup$ 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? $\endgroup$ – J. M.'s ennui Oct 17 '15 at 16:42
  • 2
    $\begingroup$ 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. $\endgroup$ – Jens Oct 17 '15 at 17:54
  • $\begingroup$ What about if the matrices are sparse? Also, where do I get RepeatedTiming? $\endgroup$ – fred Oct 19 '15 at 2:20
  • $\begingroup$ @fred, RepeatedTiming is a new function in version 10.1. $\endgroup$ – J. M.'s ennui Oct 20 '15 at 3:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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