28
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Dec 4, 2017 at 12:41
  • $\begingroup$ 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. $\endgroup$ Commented Dec 4, 2017 at 18:17
  • 4
    $\begingroup$ 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.... $\endgroup$ Commented Dec 12, 2017 at 12:42
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Mar 21, 2018 at 5:45
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Mar 22, 2018 at 2:42

3 Answers 3

Reset to default
22
$\begingroup$

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...

Improve this answer
$\endgroup$
2
  • $\begingroup$ Undocumented? Things are becoming more interesting! BTW, a natural question arises as where to find the documentation of undocumented usages/features? $\endgroup$ Commented Dec 12, 2017 at 13:44
  • $\begingroup$ @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... $\endgroup$ Commented Dec 12, 2017 at 15:41
21
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Dec 4, 2017 at 23:01
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Dec 5, 2017 at 5:07
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Dec 6, 2017 at 14:34
  • 1
    $\begingroup$ @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). $\endgroup$ Commented Dec 12, 2017 at 13:45
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Mar 21, 2018 at 7:45
13
$\begingroup$

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}}
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.