Sometimes, to find the inverse of a matrix is a labor-consuming task, or even "disgusting", especially when the matrix is "ill". It is said, e.g. in Julia, that the left/right division operation is more stable than the bare inv. So, does there exist a similar functionality in Mathematica?

Of very close relation, there are at least three functions that deal with matrix inversion in Wolfram: Inverse, PseudoInverse and LinearSolve (LinearSolve[#, IdentityMatrix[Length[#]]] &). Then according to what principles can one make the decision choosing which to use?


Update

The answers make me somewhat "embarrassed" because they all are very informative and instructive so that it is difficult to choose one to make a green tick.

  • 2
    The Julia documentation has some info on what criteria they use to determine which algorithm to use. I think the closest to that is LinearSolve, which also chooses the best algorithm (Method) depending on the type of matrix you give to it. – Thies Heidecke Dec 4 '17 at 12:41
  • The pseudoinverse, aka the Moore-Penrose inverse, is not exactly the same as the inverse. If you need the inverse, use Inverse. If LinearSolve would be the way to go, clearly this is how Inverse would have been implemented. – C. E. Dec 4 '17 at 18:17
  • 3
    I've been wishing someone would emphasize that LinearSolve[A, b] is the closest Mathematica equivalent to A \ b, other than in a link. But it doesn't seem worth a new answer.... – Michael E2 Dec 12 '17 at 12:42
  • 2
    You forget LeastSquares[] in your list. Use that if you don't have square matrices, but you want a minimum norm solution to your linear equation. – J. M. is computer-less Mar 21 at 5:45
  • 1
    @J.M. It was not forgotten, but in fact I had not known it. But it is plausibly included in "at least". Anyway, thank you for letting me hear of it. – Αλέξανδρος Ζεγγ Mar 22 at 2:42

The documentation for MATLAB's mrdivide states that

The operators / and \ are related to each other by the equation B/A = (A'\B')'

In light of this we may write:

LinearSolve[Transpose[A], b]

This question was also answered by Daniel Lichtblau here.

The equivalency follows from the following two properties of the transpose: $$ \left(A^\mathrm{-1}\right)^\mathrm{T} = \left(A^\mathrm{T}\right)^\mathrm{-1},\quad \left(\mathbf{A}\mathbf{B}\right)^\mathrm{T}=\mathbf{B}^\mathrm{T}\mathbf{A}^\mathrm{T} $$ We can show it like this: $$ x\mathbf{A} = \mathbf{B}\Leftrightarrow x = \mathbf{B}\mathbf{A}^{-1}\Leftrightarrow x^\mathrm{T}=\left(\mathbf{B}\mathbf{A}^{-1}\right)^\mathrm{T}=\left(\mathbf{A}^{-1}\right)^\mathrm{T}\mathbf{B}^\mathrm{T}=\left(\mathbf{A}^{\mathrm{T}}\right)^{-1}\mathbf{B}^\mathrm{T}\\ \Leftrightarrow \mathbf{A}^\mathrm{T}x^\mathrm{T} = \mathbf{B}^\mathrm{T} $$

Due to the way Mathematica handles the concept of row and column vectors, we do not need to use Transpose on vectors. This leads to the solution above.

  • 2
    Still, this requires the transposition of the matrix A. Even more severe: If A has already been factorized (e.g., LU-factorization with LinearSolve[A]), the Transpose[A] had to be factorized, again. This can be a performance issue (in particular with UMFPACK that is used as the default for SparseArrays. But common direct solvers based on LU-decomposition (LAPACK, Pardiso...) also provide a solver for the transpose that reuses the first factorization of A. What's merely lacking (as far as I know) is a user working interface for that. – Henrik Schumacher Dec 4 '17 at 23:01
  • 1
    @HenrikSchumacher Yes, Pardiso's API is underutilized. (And the Pardiso method is not documented yet.) However, this is not very surprising from such a high-level language as WL. Other considerations than such performance considerations take precedence. – C. E. Dec 5 '17 at 5:07
  • 1
    @CElliott I think there is a misunderstanding of what this answer does. The point of this answer is indeed that it solves xA = b without explicitly inverting A, which is not desirable for various reason. It achieves this by using LinearSolve. How LinearSolve solves the problem depends on the method chosen, but a decomposition such as the QR decomposition or the Cholesky decomposition will probably be involved, yes. – C. E. Dec 6 '17 at 14:34
  • 1
    @CElliott Here is an apology for inv(A)*b. Note that for regression, one would probably use LinearModelFit, whose methods are unknown to me but seem to never lose to QR on random examples (the SS for LinearModelFit are sometimes a couple of machine epsilon less than QR). – Michael E2 Dec 12 '17 at 13:45
  • 1
    @Michael, LinearModelFit[] (and Fit[], and FindFit[] in the linear case) all use SVD, which is definitely much more stable (tho more expensive) than QR for least squares. – J. M. is computer-less Mar 21 at 7:45

I just stumbled upon it! At least when storing the factorization in a LinearSolveFunction object, we can use it for the transposed solve by supplying a further (not documented?) string variable to it: When calling sol[b, 1234] with a LinearSolveFunction object sol, Mathematica tells us that it only accepts the strings "N", "T", "C", and "J". As it is the LAPACK standard, "N" means solving with no transposition and "T" means solving with transposition. "C" means conjugate transpose, and "J" means conjugate (maybe because "J" is a reflection of "C" in the same way as the transpose of a matrix is its reflection along the main diagonal?). Here a usage example:

SeedRandom[123];
n = 5;
A = RandomComplex[{-1 - I, 1 + I}, {n, n}];
b = RandomComplex[{-1 - I, 1 + I}, {n}];
sol = LinearSolve[A];
x = sol[b];
y = sol[b, "T"];
z = sol[b, "C"];
w = sol[b, "J"];

Max[Abs[A.x - b]]
Max[Abs[y.A - b]]
Max[Abs[z.Conjugate[A] - b]]
Max[Abs[Conjugate[A].w - b]]

3.33067*10^-16

4.8473*10^-16

4.44089*10^-16

3.14018*10^-16

A word of warning:

This seems to work with LAPACK as backend for dense matrices and with UMFPACK (Method->"Multifrontal") and the nonsymmetric iterative solvers (Method->{"Krylov",Method->"GMRES" and Method->{"Krylov",Method->"BiCGSTAB") as backend for sparse matrices but not with the (also almost undocumented) Pardiso solver (Method->"Pardiso"), although Pardiso has these capabilities and although it would have been easy to implement it...

  • Undocumented? Things are becoming more interesting! BTW, a natural question arises as where to find the documentation of undocumented usages/features? – Αλέξανδρος Ζεγγ Dec 12 '17 at 13:44
  • @AlexanderZeng That's the sad truth: The API for the various linear solvers are rarely documented. If found this feature by pure luck due to a typo and recalled that there was this open question... – Henrik Schumacher Dec 12 '17 at 15:41

Henrik has already written on the undocumented second argument of a LinearSolveFunction[], so let me just put out a short demo wrapper showing how to use the built-in, but undocumented LAPACK routines to solve a linear equation, given the output of LUDecomposition[]:

Options[backsub] = {Mode -> Automatic};
backsub[lu_, perm_, opts : OptionsPattern[]][rhs_] := Module[{xx = rhs, piv, switch},
    piv = LinearAlgebra`LAPACK`PermutationToPivot[perm];
    switch = Switch[OptionValue[Mode],
                    Normal | Automatic, "N",
                    Transpose, "T",
                    ConjugateTranspose, "C",
                    _, "N"];
    LinearAlgebra`LAPACK`GETRS[switch, lu, piv, xx]; xx]

For example, after using LUDecomposition[] as usual,

mat = N[{{1, 2 - I, 3}, {1 + 4 I, 5, 6 + 3 I}, {7 - 5 I, 8 - 2 I, 9}}];
{lu, perm, cond} = LUDecomposition[mat]
   {{{7. - 5. I, 8. - 2. I, 9. + 0. I},
     {-0.175676 + 0.445946 I, 5.51351 - 3.91892 I, 7.58108 - 1.01351 I},
     {0.0945946 + 0.0675676 I, 0.249262 - 0.0679268 I, 0.32782 + 0.15948 I}},
    {3, 2, 1}, 115.139}

we can do this:

new = {{14 - 2 I, 24 - 7 I, 24 + 7 I},
       {29 + 13 I, 36 - 7 I, 36 + 7 I},
       {50 - 9 I, 42 + 6 I, 42 - 6 I}};

bs = backsub[lu, perm];
bs[new] // Chop
   {{1., -3.2 + 7.4 I, -8.8 + 4.6 I},
    {2., -25.4 + 3.4 I, -25. - 32.2 I},
    {3., 24.8667 - 15.5333 I, 38.3333 + 13.9333 I}}

bsh = backsub[lu, perm, Mode -> ConjugateTranspose];
bsh[new] // Chop
   {{40.1333 - 9.2 I, -21.8 + 74.4 I, 1.},
    {-6.86667 + 0.466667 I, 13.2 - 17.6 I, 2.},
    {-3.4 - 0.533333 I, 9. - 8. I, 3.}}

bst = backsub[lu, perm, Mode -> Transpose];
bst[new] // Chop
   {{23.9333 - 115.4 I, 1., -21.8 - 74.4 I},
    {2.93333 + 32.9333 I, 2., 13.2 + 17.6 I},
    {6.6 + 14.5333 I, 3., 9. + 8. I}}

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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