12
$\begingroup$

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.

Improve this question
$\endgroup$
9
  • $\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 eventual burnout ♦
    Oct 17, 2015 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 eventual burnout ♦
    Oct 17, 2015 at 2:22
  • $\begingroup$ There is also the BLAS function LinearAlgebra`BLAS`TRSV[], which might be used internally by LinearSolveFunction[]. $\endgroup$
    – J. M.'s eventual burnout ♦
    Oct 17, 2015 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 eventual burnout ♦
    Oct 18, 2015 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, 2015 at 12:56

2 Answers 2

Reset to default
7
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • 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, 2015 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 eventual burnout ♦
    Oct 20, 2015 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, 2015 at 9:41
  • $\begingroup$ @thank_you Have you seen this? $\endgroup$
    – Michael E2
    Jan 29, 2017 at 17:52
3
$\begingroup$

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

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

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.