Existence and uniqueness

As discussed above, for any matrix ${\displaystyle A}$ there is one and only one pseudoinverse ${\displaystyle A^{+}}$.^[8]

A matrix satisfying only the first of the conditions given above, namely ${\textstyle AA^{+}A=A}$, it is known as a generalized inverse. If the matrix also satisfies the second condition, namely ${\textstyle A^{+}AA^{+}=A^{+}}$, it is called a generalized reflexive inverse. Generalized inverses always exist but are not in general unique. Uniqueness is a consequence of the last two conditions.

Basic properties

Proofs for the properties below can be found at b:Topics in Abstract Algebra/Linear algebra.

• If ${\displaystyle A}$ has real entries, then so does ${\displaystyle A^{+}}$.
• If ${\displaystyle A}$ is invertible, its pseudoinverse is its inverse. That is, ${\displaystyle A^{+}=A^{-1}}$.^[9]: 243
• The pseudoinverse of the pseudoinverse is the original matrix: ${\displaystyle \left(A^{+}\right)^{+}=A}$.^[9]: 245
• Pseudoinversion commutes with transposition, complex conjugation, and taking the conjugate transpose:^[9]: 245
  $${\displaystyle \left(A^{\operatorname {T} }\right)^{+}=\left(A^{+}\right)^{\operatorname {T} },\quad \left({\overline {A}}\right)^{+}={\overline {A^{+}}},\quad \left(A^{*}\right)^{+}=\left(A^{+}\right)^{*}.}$$

  
• The pseudoinverse of a scalar multiple of ${\displaystyle A}$ is the reciprocal multiple of ${\displaystyle A^{+}}$:
  $${\displaystyle \left(\alpha A\right)^{+}=\alpha ^{-1}A^{+}}$$

  
  for ${\displaystyle \alpha \neq 0}$.
• Ker and image of the pseudoinverse coincide with those of the conjugate transpose: ${\displaystyle \ker \left(A^{+}\right)=\ker \left(A^{*}\right)}$ and ${\displaystyle \operatorname {ran} \left(A^{+}\right)=\operatorname {ran} \left(A^{*}\right)}$.

Identities

The following identity formula can be used to cancel or expand certain subexpressions involving pseudoinverses:

$${\displaystyle A={}A{}A^{*}{}A^{+*}{}={}A^{+*}{}A^{*}{}A.}$$

Equivalently, substituting ${\displaystyle A^{+}}$ for ${\displaystyle A}$ gives

$${\displaystyle A^{+}={}A^{+}{}A^{+*}{}A^{*}{}={}A^{*}{}A^{+*}{}A^{+},}$$

while substituting ${\displaystyle A^{*}}$ for ${\displaystyle A}$ gives

$${\displaystyle A^{*}={}A^{*}{}A{}A^{+}{}={}A^{+}{}A{}A^{*}.}$$

Reduction to Hermitian case

The computation of the pseudoinverse is reducible to its construction in the Hermitian case. This is possible through the equivalences:

$${\displaystyle A^{+}=\left(A^{*}A\right)^{+}A^{*},}$$

$${\displaystyle A^{+}=A^{*}\left(AA^{*}\right)^{+},}$$

as ${\displaystyle A^{*}A}$ and ${\displaystyle AA^{*}}$ are Hermitian.

Pseudoinverse of products

The equality ${\displaystyle (AB)^{+}=B^{+}A^{+}}$ does not hold in general. Rather, suppose ${\displaystyle A\in \mathbb {K} ^{m\times n},\ B\in \mathbb {K} ^{n\times p}}$. Then the following are equivalent:^[10]

1. $${\displaystyle (AB)^{+}=B^{+}A^{+}}$$

   
2. $${\displaystyle {\begin{aligned}A^{+}ABB^{*}A^{*}&=BB^{*}A^{*},\\BB^{+}A^{*}AB&=A^{*}AB.\end{aligned}}}$$

   
3. $${\displaystyle {\begin{aligned}\left(A^{+}ABB^{*}\right)^{*}&=A^{+}ABB^{*},\\\left(A^{*}ABB^{+}\right)^{*}&=A^{*}ABB^{+}.\end{aligned}}}$$

   
4. $${\displaystyle A^{+}ABB^{*}A^{*}ABB^{+}=BB^{*}A^{*}A}$$

   
5. $${\displaystyle {\begin{aligned}A^{+}AB&=B(AB)^{+}AB,\\BB^{+}A^{*}&=A^{*}AB(AB)^{+}.\end{aligned}}}$$

   

The following are sufficient conditions for ${\displaystyle (AB)^{+}=B^{+}A^{+}}$:

1. ${\displaystyle A}$ has orthonormal columns (then ${\displaystyle A^{*}A=A^{+}A=I_{n}}$),   or
2. ${\displaystyle B}$ has orthonormal rows (then ${\displaystyle BB^{*}=BB^{+}=I_{n}}$),   or
3. ${\displaystyle A}$ has linearly independent columns (then ${\displaystyle A^{+}A=I}$ ) and ${\displaystyle B}$ has linearly independent rows (then ${\displaystyle BB^{+}=I}$),   or
4. ${\displaystyle B=A^{*}}$, or
5. ${\displaystyle B=A^{+}}$.

The following is a necessary condition for ${\displaystyle (AB)^{+}=B^{+}A^{+}}$:

1. ${\displaystyle (A^{+}A)(BB^{+})=(BB^{+})(A^{+}A)}$

The fourth sufficient condition yields the equalities

$${\displaystyle {\begin{aligned}\left(AA^{*}\right)^{+}&=A^{+*}A^{+},\\\left(A^{*}A\right)^{+}&=A^{+}A^{+*}.\end{aligned}}}$$

Here is a counterexample where ${\displaystyle (AB)^{+}\neq B^{+}A^{+}}$:

$${\displaystyle {\Biggl (}{\begin{pmatrix}1&1\\0&0\end{pmatrix}}{\begin{pmatrix}0&0\\1&1\end{pmatrix}}{\Biggr )}^{+}={\begin{pmatrix}1&1\\0&0\end{pmatrix}}^{+}={\begin{pmatrix}{\tfrac {1}{2}}&0\\{\tfrac {1}{2}}&0\end{pmatrix}}\quad \neq \quad {\begin{pmatrix}{\tfrac {1}{4}}&0\\{\tfrac {1}{4}}&0\end{pmatrix}}={\begin{pmatrix}0&{\tfrac {1}{2}}\\0&{\tfrac {1}{2}}\end{pmatrix}}{\begin{pmatrix}{\tfrac {1}{2}}&0\\{\tfrac {1}{2}}&0\end{pmatrix}}={\begin{pmatrix}0&0\\1&1\end{pmatrix}}^{+}{\begin{pmatrix}1&1\\0&0\end{pmatrix}}^{+}}$$

Projectors

${\displaystyle P=AA^{+}}$ and ${\displaystyle Q=A^{+}A}$ are orthogonal projection operators, that is, they are Hermitian (${\displaystyle P=P^{*}}$, ${\displaystyle Q=Q^{*}}$) and idempotent (${\displaystyle P^{2}=P}$ and ${\displaystyle Q^{2}=Q}$). The following hold:

• ${\displaystyle PA=AQ=A}$ and ${\displaystyle A^{+}P=QA^{+}=A^{+}}$
• ${\displaystyle P}$ is the orthogonal projector onto the range of ${\displaystyle A}$ (which equals the orthogonal complement of the kernel of ${\displaystyle A^{*}}$).
• ${\displaystyle Q}$ is the orthogonal projector onto the range of ${\displaystyle A^{*}}$ (which equals the orthogonal complement of the kernel of ${\displaystyle A}$).
• ${\displaystyle I-Q=I-A^{+}A}$ is the orthogonal projector onto the kernel of ${\displaystyle A}$.
• ${\displaystyle I-P=I-AA^{+}}$ is the orthogonal projector onto the kernel of ${\displaystyle A^{*}}$.^[8]

The last two properties imply the following identities:

• ${\displaystyle A\,\ \left(I-A^{+}A\right)=\left(I-AA^{+}\right)A\ \ =0}$
• ${\displaystyle A^{*}\left(I-AA^{+}\right)=\left(I-A^{+}A\right)A^{*}=0}$

Another property is the following: if ${\displaystyle A\in \mathbb {K} ^{n\times n}}$ is Hermitian and idempotent (true if and only if it represents an orthogonal projection), then, for any matrix ${\displaystyle B\in \mathbb {K} ^{m\times n}}$ the following equation holds:^[11]

$${\displaystyle A(BA)^{+}=(BA)^{+}}$$

This can be proven by defining matrices ${\displaystyle C=BA}$, ${\displaystyle D=A(BA)^{+}}$, and checking that ${\displaystyle D}$ is indeed a pseudoinverse for ${\displaystyle C}$ by verifying that the defining properties of the pseudoinverse hold, when ${\displaystyle A}$ is Hermitian and idempotent.

From the last property it follows that, if ${\displaystyle A\in \mathbb {K} ^{n\times n}}$ is Hermitian and idempotent, for any matrix ${\displaystyle B\in \mathbb {K} ^{n\times m}}$

$${\displaystyle (AB)^{+}A=(AB)^{+}}$$

Finally, if ${\displaystyle A}$ is an orthogonal projection matrix, then its pseudoinverse trivially coincides with the matrix itself, that is, ${\displaystyle A^{+}=A}$.

Geometric construction

If we view the matrix as a linear map ${\displaystyle A:\mathbb {K} ^{n}\to \mathbb {K} ^{m}}$ over the field ${\displaystyle \mathbb {K} }$ then ${\displaystyle A^{+}:\mathbb {K} ^{m}\to \mathbb {K} ^{n}}$ can be decomposed as follows. We write ${\displaystyle \oplus }$ for the direct sum, ${\displaystyle \perp }$ for the orthogonal complement, ${\displaystyle \ker }$ for the kernel of a map, and ${\displaystyle \operatorname {ran} }$ for the image of a map. Notice that ${\displaystyle \mathbb {K} ^{n}=\left(\ker A\right)^{\perp }\oplus \ker A}$ and ${\displaystyle \mathbb {K} ^{m}=\operatorname {ran} A\oplus \left(\operatorname {ran} A\right)^{\perp }}$. The restriction ${\displaystyle A:\left(\ker A\right)^{\perp }\to \operatorname {ran} A}$ is then an isomorphism. This implies that ${\displaystyle A^{+}}$ on ${\displaystyle \operatorname {ran} A}$ is the inverse of this isomorphism, and is zero on ${\displaystyle \left(\operatorname {ran} A\right)^{\perp }.}$

In other words: To find ${\displaystyle A^{+}b}$ for given ${\displaystyle b}$ in ${\displaystyle \mathbb {K} ^{m}}$, first project ${\displaystyle b}$ orthogonally onto the range of ${\displaystyle A}$, finding a point ${\displaystyle p(b)}$ in the range. Then form ${\displaystyle A^{-1}(\{p(b)\})}$, that is, find those vectors in ${\displaystyle \mathbb {K} ^{n}}$ that ${\displaystyle A}$ sends to ${\displaystyle p(b)}$. This will be an affine subspace of ${\displaystyle \mathbb {K} ^{n}}$ parallel to the kernel of ${\displaystyle A}$. The element of this subspace that has the smallest length (that is, is closest to the origin) is the answer ${\displaystyle A^{+}b}$ we are looking for. It can be found by taking an arbitrary member of ${\displaystyle A^{-1}(\{p(b)\})}$ and projecting it orthogonally onto the orthogonal complement of the kernel of ${\displaystyle A}$.

This description is closely related to the minimum-norm solution to a linear system.

Limit relations

The pseudoinverse are limits:

$${\displaystyle A^{+}=\lim _{\delta \searrow 0}\left(A^{*}A+\delta I\right)^{-1}A^{*}=\lim _{\delta \searrow 0}A^{*}\left(AA^{*}+\delta I\right)^{-1}}$$

(see Tikhonov regularization). These limits exist even if ${\displaystyle \left(AA^{*}\right)^{-1}}$ or ${\displaystyle \left(A^{*}A\right)^{-1}}$ do not exist.^[8]: 263

Continuity

In contrast to ordinary matrix inversion, the process of taking pseudoinverses is not continuous: if the sequence ${\displaystyle \left(A_{n}\right)}$ converges to the matrix ${\displaystyle A}$ (in the maximum norm or Frobenius norm, say), then ${\displaystyle (A_{n})^{+}}$ need not converge to ${\displaystyle A^{+}}$. However, if all the matrices ${\displaystyle A_{n}}$ have the same rank as ${\displaystyle A}$, ${\displaystyle (A_{n})^{+}}$ will converge to ${\displaystyle A^{+}}$.^[12]

Derivative

The derivative of a real-valued pseudoinverse matrix that has constant rank at a point ${\displaystyle x}$ may be calculated in terms of the derivative of the original matrix:^[13]

$${\displaystyle {\frac {\mathrm {d} }{\mathrm {d} x}}A^{+}(x)=-A^{+}\left({\frac {\mathrm {d} }{\mathrm {d} x}}A\right)A^{+}~+~A^{+}A^{+{\textsf {T}}}\left({\frac {\mathrm {d} }{\mathrm {d} x}}A^{\operatorname {T} }\right)\left(I-AA^{+}\right)~+~\left(I-A^{+}A\right)\left({\frac {\text{d}}{{\text{d}}x}}A^{\operatorname {T} }\right)A^{+{\textsf {T}}}A^{+}.}$$

For a complex matrix, the transpose is replaced with conjugate transpose.^[14] For a real-valued symmetric matrix, the Magnus-Neudecker derivative is established.^[15]

Scalars

It is also possible to define a pseudoinverse for scalars and vectors. This amounts to treating these as matrices. The pseudoinverse of a scalar ${\displaystyle x}$ is zero if ${\displaystyle x}$ is zero and the reciprocal of ${\displaystyle x}$ otherwise:

$${\displaystyle x^{+}={\begin{cases}0,&{\mbox{if }}x=0;\\x^{-1},&{\mbox{otherwise}}.\end{cases}}}$$

Vectors

The pseudoinverse of the null (all zero) vector is the transposed null vector. The pseudoinverse of a non-null vector is the conjugate transposed vector divided by its squared magnitude:

$${\displaystyle {\vec {x}}^{+}={\begin{cases}{\vec {0}}^{\operatorname {T} },&{\text{if }}{\vec {x}}={\vec {0}};\\{\dfrac {{\vec {x}}^{*}}{{\vec {x}}^{*}{\vec {x}}}},&{\text{otherwise}}.\end{cases}}}$$

Diagonal matrices

The pseudoinverse of a squared diagonal matrix is obtained by taking the reciprocal of the nonzero diagonal elements. Formally, if ${\displaystyle D}$ is a squared diagonal matrix with ${\displaystyle D={\tilde {D}}\oplus \mathbf {0} _{k\times k}}$ and ${\displaystyle {\tilde {D}}>0}$, then ${\displaystyle D^{+}={\tilde {D}}^{-1}\oplus \mathbf {0} _{k\times k}}$. More generally, if ${\displaystyle A}$ is any ${\displaystyle m\times n}$ rectangular matrix whose only nonzero elements are on the diagonal, meaning ${\displaystyle A_{ij}=\delta _{ij}a_{i}}$, ${\displaystyle a_{i}\in \mathbb {K} }$, then ${\displaystyle A^{+}}$ is a ${\displaystyle n\times m}$ rectangular matrix whose diagonal elements are the reciprocal of the original ones, that is, ${\displaystyle A_{ii}\neq 0\implies A_{ii}^{+}=1/A_{ii}}$.

Linearly independent columns

If the rank of ${\displaystyle A}$ is identical to its column rank, ${\displaystyle n}$, (for ${\displaystyle n\leq m}$,) there are ${\displaystyle n}$ linearly independent columns, and ${\displaystyle A^{*}A}$ is invertible. In this case, an explicit formula is:^[16]

$${\displaystyle A^{+}=\left(A^{*}A\right)^{-1}A^{*}.}$$

It follows that ${\displaystyle A^{+}}$ is then a left inverse of ${\displaystyle A}$:   ${\displaystyle A^{+}A=I_{n}}$.

Linearly independent rows

If the rank of ${\displaystyle A}$ is identical to its row rank, ${\displaystyle m}$, (for ${\displaystyle m\leq n}$,) there are ${\displaystyle m}$ linearly independent rows, and ${\displaystyle AA^{*}}$ is invertible. In this case, an explicit formula is:

$${\displaystyle A^{+}=A^{*}\left(AA^{*}\right)^{-1}.}$$

It follows that ${\displaystyle A^{+}}$ is a right inverse of ${\displaystyle A}$:   ${\displaystyle AA^{+}=I_{m}}$.

Orthonormal columns or rows

This is a special case of either full column rank or full row rank (treated above). If ${\displaystyle A}$ has orthonormal columns (${\displaystyle A^{*}A=I_{n}}$) or orthonormal rows (${\displaystyle AA^{*}=I_{m}}$), then:

$${\displaystyle A^{+}=A^{*}.}$$

Normal matrices

If ${\displaystyle A}$ is normal, that is, it commutes with its conjugate transpose, then its pseudoinverse can be computed by diagonalizing it, mapping all nonzero eigenvalues to their inverses, and mapping zero eigenvalues to zero. A corollary is that ${\displaystyle A}$ commuting with its transpose implies that it commutes with its pseudoinverse.

Orthogonal projection matrices

This is a special case of a normal matrix with eigenvalues 0 and 1. If ${\displaystyle A}$ is an orthogonal projection matrix, that is, ${\displaystyle A=A^{*}}$ and ${\displaystyle A^{2}=A}$, then the pseudoinverse trivially coincides with the matrix itself:

$${\displaystyle A^{+}=A.}$$

Circulant matrices

For a circulant matrix ${\displaystyle C}$, the singular value decomposition is given by the Fourier transform, that is, the singular values are the Fourier coefficients. Let ${\displaystyle {\mathcal {F}}}$ be the Discrete Fourier Transform (DFT) matrix; then^[17]

$${\displaystyle {\begin{aligned}C&={\mathcal {F}}\cdot \Sigma \cdot {\mathcal {F}}^{*},\\C^{+}&={\mathcal {F}}\cdot \Sigma ^{+}\cdot {\mathcal {F}}^{*}.\end{aligned}}}$$

Rank decomposition

Let ${\displaystyle r\leq \min(m,n)}$ denote the rank of ${\displaystyle A\in \mathbb {K} ^{m\times n}}$. Then ${\displaystyle A}$ can be (rank) decomposed as ${\displaystyle A=BC}$ where ${\displaystyle B\in \mathbb {K} ^{m\times r}}$ and ${\displaystyle C\in \mathbb {K} ^{r\times n}}$ are of rank ${\displaystyle r}$. Then ${\displaystyle A^{+}=C^{+}B^{+}=C^{*}\left(CC^{*}\right)^{-1}\left(B^{*}B\right)^{-1}B^{*}}$.

The QR method

For ${\displaystyle \mathbb {K} \in \{\mathbb {R} ,\mathbb {C} \}}$ computing the product ${\displaystyle AA^{*}}$ or ${\displaystyle A^{*}A}$ and their inverses explicitly is often a source of numerical rounding errors and computational cost in practice. An alternative approach using the QR decomposition of ${\displaystyle A}$ may be used instead.

Consider the case when ${\displaystyle A}$ is of full column rank, so that ${\displaystyle A^{+}=\left(A^{*}A\right)^{-1}A^{*}}$. Then the Cholesky decomposition ${\displaystyle A^{*}A=R^{*}R}$, where ${\displaystyle R}$ is an upper triangular matrix, may be used. Multiplication by the inverse is then done easily by solving a system with multiple right-hand sides,

$${\displaystyle A^{+}=\left(A^{*}A\right)^{-1}A^{*}\quad \Leftrightarrow \quad \left(A^{*}A\right)A^{+}=A^{*}\quad \Leftrightarrow \quad R^{*}RA^{+}=A^{*}}$$

which may be solved by forward substitution followed by back substitution.

The Cholesky decomposition may be computed without forming ${\displaystyle A^{*}A}$ explicitly, by alternatively using the QR decomposition of ${\displaystyle A=QR}$, where ${\displaystyle Q}$ has orthonormal columns, ${\displaystyle Q^{*}Q=I}$, and ${\displaystyle R}$ is upper triangular. Then

$${\displaystyle A^{*}A\,=\,(QR)^{*}(QR)\,=\,R^{*}Q^{*}QR\,=\,R^{*}R,}$$

so ${\displaystyle R}$ is the Cholesky factor of ${\displaystyle A^{*}A}$.

The case of full row rank is treated similarly by using the formula ${\displaystyle A^{+}=A^{*}\left(AA^{*}\right)^{-1}}$ and using a similar argument, swapping the roles of ${\displaystyle A}$ and ${\displaystyle A^{*}}$.

Singular value decomposition (SVD)

A computationally simple and accurate way to compute the pseudoinverse is by using the singular value decomposition.^[16][8][18] If ${\displaystyle A=U\Sigma V^{*}}$ is the singular value decomposition of ${\displaystyle A}$, then ${\displaystyle A^{+}=V\Sigma ^{+}U^{*}}$. For a rectangular diagonal matrix such as ${\displaystyle \Sigma }$, we get the pseudoinverse by taking the reciprocal of each non-zero element on the diagonal, leaving the zeros in place, and then transposing the matrix. In numerical computation, only elements larger than some small tolerance are taken to be nonzero, and the others are replaced by zeros. For example, in the MATLAB or GNU Octave function pinv, the tolerance is taken to be t = ε⋅max(m, n)⋅max(Σ), where ε is the machine epsilon.

The computational cost of this method is dominated by the cost of computing the SVD, which is several times higher than matrix–matrix multiplication, even if a state-of-the art implementation (such as that of LAPACK) is used.

The above procedure shows why taking the pseudoinverse is not a continuous operation: if the original matrix ${\displaystyle A}$ has a singular value 0 (a diagonal entry of the matrix ${\displaystyle \Sigma }$ above), then modifying ${\displaystyle A}$ slightly may turn this zero into a tiny positive number, thereby affecting the pseudoinverse dramatically as we now have to take the reciprocal of a tiny number.

Block matrices

Optimized approaches exist for calculating the pseudoinverse of block-structured matrices.

The iterative method of Ben-Israel and Cohen

Another method for computing the pseudoinverse (cf. Drazin inverse) uses the recursion

$${\displaystyle A_{i+1}=2A_{i}-A_{i}AA_{i},}$$

which is sometimes referred to as hyper-power sequence. This recursion produces a sequence converging quadratically to the pseudoinverse of ${\displaystyle A}$ if it is started with an appropriate ${\displaystyle A_{0}}$ satisfying ${\displaystyle A_{0}A=\left(A_{0}A\right)^{*}}$. The choice ${\displaystyle A_{0}=\alpha A^{*}}$ (where ${\displaystyle 0<\alpha <2/\sigma _{1}^{2}(A)}$, with ${\displaystyle \sigma _{1}(A)}$ denoting the largest singular value of ${\displaystyle A}$)^[19] has been argued not to be competitive to the method using the SVD mentioned above, because even for moderately ill-conditioned matrices it takes a long time before ${\displaystyle A_{i}}$ enters the region of quadratic convergence.^[20] However, if started with ${\displaystyle A_{0}}$ already close to the Moore–Penrose inverse and ${\displaystyle A_{0}A=\left(A_{0}A\right)^{*}}$, for example ${\displaystyle A_{0}:=\left(A^{*}A+\delta I\right)^{-1}A^{*}}$, convergence is fast (quadratic).

Updating the pseudoinverse

For the cases where ${\displaystyle A}$ has full row or column rank, and the inverse of the correlation matrix (${\displaystyle AA^{*}}$ for ${\displaystyle A}$ with full row rank or ${\displaystyle A^{*}A}$ for full column rank) is already known, the pseudoinverse for matrices related to ${\displaystyle A}$ can be computed by applying the Sherman–Morrison–Woodbury formula to update the inverse of the correlation matrix, which may need less work. In particular, if the related matrix differs from the original one by only a changed, added or deleted row or column, incremental algorithms exist that exploit the relationship.^[21][22]

Similarly, it is possible to update the Cholesky factor when a row or column is added, without creating the inverse of the correlation matrix explicitly. However, updating the pseudoinverse in the general rank-deficient case is much more complicated.^[23][24]

Software libraries

High-quality implementations of SVD, QR, and back substitution are available in standard libraries, such as LAPACK. Writing one's own implementation of SVD is a major programming project that requires a significant numerical expertise. In special circumstances, such as parallel computing or embedded computing, however, alternative implementations by QR or even the use of an explicit inverse might be preferable, and custom implementations may be unavoidable.

The Python package NumPy provides a pseudoinverse calculation through its functions matrix.I and linalg.pinv; its pinv uses the SVD-based algorithm. SciPy adds a function scipy.linalg.pinv that uses a least-squares solver.

The MASS package for R provides a calculation of the Moore–Penrose inverse through the ginv function.^[25] The ginv function calculates a pseudoinverse using the singular value decomposition provided by the svd function in the base R package. An alternative is to employ the pinv function available in the pracma package.

The Octave programming language provides a pseudoinverse through the standard package function pinv and the pseudo_inverse() method.

In Julia (programming language), the LinearAlgebra package of the standard library provides an implementation of the Moore–Penrose inverse pinv() implemented via singular-value decomposition.^[26]

Linear least-squares

The pseudoinverse provides a least squares solution to a system of linear equations.^[27] For ${\displaystyle A\in \mathbb {K} ^{m\times n}}$, given a system of linear equations

$${\displaystyle Ax=b,}$$

in general, a vector ${\displaystyle x}$ that solves the system may not exist, or if one does exist, it may not be unique. More specifically, a solution exists if and only if ${\displaystyle b}$ is in the image of ${\displaystyle A}$, and is unique if and only if ${\displaystyle A}$ is injective. The pseudoinverse solves the "least-squares" problem as follows:

• ${\displaystyle \forall x\in \mathbb {K} ^{n}}$, we have ${\displaystyle \left\|Ax-b\right\|_{2}\geq \left\|Az-b\right\|_{2}}$ where ${\displaystyle z=A^{+}b}$ and ${\displaystyle \|\cdot \|_{2}}$ denotes the Euclidean norm. This weak inequality holds with equality if and only if ${\displaystyle x=A^{+}b+\left(I-A^{+}A\right)w}$ for any vector ${\displaystyle w}$; this provides an infinitude of minimizing solutions unless ${\displaystyle A}$ has full column rank, in which case ${\displaystyle \left(I-A^{+}A\right)}$ is a zero matrix.^[28] The solution with minimum Euclidean norm is ${\displaystyle z.}$^[28]

This result is easily extended to systems with multiple right-hand sides, when the Euclidean norm is replaced by the Frobenius norm. Let ${\displaystyle B\in \mathbb {K} ^{m\times p}}$.

• ${\displaystyle \forall X\in \mathbb {K} ^{n\times p}}$, we have ${\displaystyle \|AX-B\|_{\mathrm {F} }\geq \|AZ-B\|_{\mathrm {F} }}$ where ${\displaystyle Z=A^{+}B}$ and ${\displaystyle \|\cdot \|_{\mathrm {F} }}$ denotes the Frobenius norm.

Obtaining all solutions of a linear system

If the linear system

$${\displaystyle Ax=b}$$

has any solutions, they are all given by^[29]

$${\displaystyle x=A^{+}b+\left[I-A^{+}A\right]w}$$

for arbitrary vector ${\displaystyle w}$. Solution(s) exist if and only if ${\displaystyle AA^{+}b=b}$.^[29] If the latter holds, then the solution is unique if and only if ${\displaystyle A}$ has full column rank, in which case ${\displaystyle I-A^{+}A}$ is a zero matrix. If solutions exist but ${\displaystyle A}$ does not have full column rank, then we have an indeterminate system, all of whose infinitude of solutions are given by this last equation.

Minimum norm solution to a linear system

For linear systems ${\displaystyle Ax=b,}$ with non-unique solutions (such as under-determined systems), the pseudoinverse may be used to construct the solution of minimum Euclidean norm ${\displaystyle \|x\|_{2}}$ among all solutions.

• If ${\displaystyle Ax=b}$ is satisfiable, the vector ${\displaystyle z=A^{+}b}$ is a solution, and satisfies ${\displaystyle \|z\|_{2}\leq \|x\|_{2}}$ for all solutions.

This result is easily extended to systems with multiple right-hand sides, when the Euclidean norm is replaced by the Frobenius norm. Let ${\displaystyle B\in \mathbb {K} ^{m\times p}}$.

• If ${\displaystyle AX=B}$ is satisfiable, the matrix ${\displaystyle Z=A^{+}B}$ is a solution, and satisfies ${\displaystyle \|Z\|_{\mathrm {F} }\leq \|X\|_{\mathrm {F} }}$ for all solutions.

Condition number

Using the pseudoinverse and a matrix norm, one can define a condition number for any matrix:

$${\displaystyle {\mbox{cond}}(A)=\|A\|\left\|A^{+}\right\|.}$$

A large condition number implies that the problem of finding least-squares solutions to the corresponding system of linear equations is ill-conditioned in the sense that small errors in the entries of ${\displaystyle A}$ can lead to huge errors in the entries of the solution.^[30]