using pinverse vs gauss-jordan to solve a possibly singular system of linear equations?

Question

so i have a possibly singular coefficient matrix A . (it's sometimes singular and sometimes not !) (in the case of singularity it has many solutions we only need one)

i know using an inverse to find the answer to a system of linear equations is an awful idea .

so it seems even in this case we should use methods like Gauss-Jordan to solve the system even when it is singular . is there any method to solve singular systems like Gauss-Jordan.(to give a answer)

and i read the matlab document on solving this types of systems and it says that we should use the pseudo-inverse function pinv like : x = pinv(A)*b

this just looks as awful as taking the inverse (since for a non-singular matrix pseudo-inverse is same as inverse). (am I correct?)

and is it the case that we don't have any other more efficient choices of finding and answer in this case ??

Answer 1

It is not generally true that the pseudoinverse is the same as the inverse. In the case that a matrix is invertible, then they coincide. The pseudoinverse can be computed in different ways; pinv in MATLAB is generally good at selecting the right algorithm for the job (typically something SVD-like).

It's important to note that for a singular or otherwise ill-conditioned problem, the resulting value will not necessarily represent the solution, but rather the solution that minimizes the value $\| Ax-b\|$.

Answer 2

The answer is framed using the tools in Mathematica, which are the most familiar set. The translation to MATLAB or Octave, etc., should be direct.

Given a system matrix $\mathbf{A}$ and a data vector $b$, solve the linear system $$ \mathbf{A} x = b $$

---

The first attempt is to try a direct linear solution

x=LinearSolve[A,b]

You can trust the application, or you can specify a solution method. Current options are

1. banded matrix solver
2. Cholesky method for positive definite Hermitian matrices
3. iterative Krylov sparse solver
4. direct sparse LU decomposition

   ---

   Next, try for a least squares solutions

x = LeastSquares[A, b]

This opens the door for the experts who have written the application which chooses a good method for solution.

If there is no solution, we have learned that the data vector is in the null space.

---

Of course, there is the pseudoinverse solution

x = PseudoInverse[A].b

As noted by @Shiyue in his post What is the best way to compute the pseudoinverse of a matrix?, the application writers use fast methods to assemble the pseudoinverse.

---

There are slower options. One is to construct the pseudoinverse manually.

{u,s,v}=SingularValueDecomposition[A] x = v.PseudoInverse[s].ConjugateTranspose[u].b

The $\mathbf{QR}$ decomposition:

{q,r}=QRDecomposition[A] x = PseudoInverse[r].Conjugate[q].b