Moore–Penrose inverse
   HOME

TheInfoList



OR:

In
mathematics Mathematics is a field of study that discovers and organizes methods, Mathematical theory, theories and theorems that are developed and Mathematical proof, proved for the needs of empirical sciences and mathematics itself. There are many ar ...
, and in particular
linear algebra Linear algebra is the branch of mathematics concerning linear equations such as :a_1x_1+\cdots +a_nx_n=b, linear maps such as :(x_1, \ldots, x_n) \mapsto a_1x_1+\cdots +a_nx_n, and their representations in vector spaces and through matrix (mathemat ...
, the Moore–Penrose inverse of a
matrix Matrix (: matrices or matrixes) or MATRIX may refer to: Science and mathematics * Matrix (mathematics), a rectangular array of numbers, symbols or expressions * Matrix (logic), part of a formula in prenex normal form * Matrix (biology), the m ...
, often called the pseudoinverse, is the most widely known generalization of the inverse matrix. It was independently described by E. H. Moore in 1920, Arne Bjerhammar in 1951, and
Roger Penrose Sir Roger Penrose (born 8 August 1931) is an English mathematician, mathematical physicist, Philosophy of science, philosopher of science and Nobel Prize in Physics, Nobel Laureate in Physics. He is Emeritus Rouse Ball Professor of Mathematics i ...
in 1955. Earlier, Erik Ivar Fredholm had introduced the concept of a pseudoinverse of integral operators in 1903. The terms ''pseudoinverse'' and ''
generalized inverse In mathematics, and in particular, algebra, a generalized inverse (or, g-inverse) of an element ''x'' is an element ''y'' that has some properties of an inverse element but not necessarily all of them. The purpose of constructing a generalized inv ...
'' are sometimes used as synonyms for the Moore–Penrose inverse of a matrix, but sometimes applied to other elements of algebraic structures which share some but not all properties expected for an
inverse element In mathematics, the concept of an inverse element generalises the concepts of opposite () and reciprocal () of numbers. Given an operation denoted here , and an identity element denoted , if , one says that is a left inverse of , and that ...
. A common use of the pseudoinverse is to compute a "best fit" ( least squares) approximate solution to a
system of linear equations In mathematics, a system of linear equations (or linear system) is a collection of two or more linear equations involving the same variable (math), variables. For example, : \begin 3x+2y-z=1\\ 2x-2y+4z=-2\\ -x+\fracy-z=0 \end is a system of th ...
that lacks an exact solution (see below under § Applications). Another use is to find the minimum ( Euclidean) norm solution to a system of linear equations with multiple solutions. The pseudoinverse facilitates the statement and proof of results in linear algebra. The pseudoinverse is defined for all rectangular matrices whose entries are real or
complex Complex commonly refers to: * Complexity, the behaviour of a system whose components interact in multiple ways so possible interactions are difficult to describe ** Complex system, a system composed of many components which may interact with each ...
numbers. Given a rectangular matrix with real or complex entries, its pseudoinverse is unique. It can be computed using the
singular value decomposition In linear algebra, the singular value decomposition (SVD) is a Matrix decomposition, factorization of a real number, real or complex number, complex matrix (mathematics), matrix into a rotation, followed by a rescaling followed by another rota ...
. In the special case where is a
normal matrix In mathematics, a complex square matrix is normal if it commutes with its conjugate transpose : :A \text \iff A^*A = AA^* . The concept of normal matrices can be extended to normal operators on infinite-dimensional normed spaces and to nor ...
(for example, a Hermitian matrix), the pseudoinverse annihilates the kernel of and acts as a traditional inverse of on the subspace
orthogonal In mathematics, orthogonality (mathematics), orthogonality is the generalization of the geometric notion of ''perpendicularity''. Although many authors use the two terms ''perpendicular'' and ''orthogonal'' interchangeably, the term ''perpendic ...
to the kernel.


Notation

In the following discussion, the following conventions are adopted. * will denote one of the fields of real or complex numbers, denoted , , respectively. The vector space of matrices over is denoted by . * For , the transpose is denoted and the Hermitian transpose (also called conjugate transpose) is denoted . If \mathbb = \mathbb, then A^* = A^\mathsf. * For , (standing for " range") denotes the column space (
image An image or picture is a visual representation. An image can be Two-dimensional space, two-dimensional, such as a drawing, painting, or photograph, or Three-dimensional space, three-dimensional, such as a carving or sculpture. Images may be di ...
) of (the space spanned by the column vectors of ) and denotes the kernel (null space) of . * For any positive integer , the
identity matrix In linear algebra, the identity matrix of size n is the n\times n square matrix with ones on the main diagonal and zeros elsewhere. It has unique properties, for example when the identity matrix represents a geometric transformation, the obje ...
is denoted .


Definition

For A \in \mathbb^, a pseudoinverse of is defined as a matrix satisfying all of the following four criteria, known as the Moore–Penrose conditions: # need not be the general identity matrix, but it maps all column vectors of to themselves: A A^+ A = \; A. # acts like a weak inverse: A^+ A A^+ = \; A^+. # is Hermitian: \left(A A^+\right)^* = \; A A^+. # is also Hermitian: \left(A^+ A\right)^* = \; A^+ A. Note that A^+A and AA^+ are idempotent operators, as follows from (AA^+)^2=A A^+ and (A^+ A)^2=A^+ A. More specifically, A^+A projects onto the image of A^T (equivalently, the span of the rows of A), and AA^+ projects onto the image of A (equivalently, the span of the columns of A). In fact, the above four conditions are fully equivalent to A^+A and AA^+ being such orthogonal projections: AA^+ projecting onto the image of A implies (A A^+)A=A, and A^+A projecting onto the image of A^T implies (A^+A)A^T=A^T. The pseudoinverse A^+ exists for any matrix A \in \mathbb^. If furthermore A is full rank, that is, its rank is , then can be given a particularly simple algebraic expression. In particular: * When has linearly independent columns (equivalently, is injective, and thus is invertible), can be computed asA^+ = \left(A^* A\right)^ A^*.This particular pseudoinverse is a ''left inverse'', that is, A^+A = I. * If, on the other hand, A has linearly independent rows (equivalently, A is surjective, and thus is invertible), can be computed asA^+ = A^* \left(A A^*\right)^.This is a ''right inverse'', as A A^+ = I. In the more general case, the pseudoinverse can be expressed leveraging the
singular value decomposition In linear algebra, the singular value decomposition (SVD) is a Matrix decomposition, factorization of a real number, real or complex number, complex matrix (mathematics), matrix into a rotation, followed by a rescaling followed by another rota ...
. Any matrix can be decomposed as A=UDV^* for some isometries U,V and diagonal nonnegative real matrix D. The pseudoinverse can then be written as A^+=V D^ U^*, where D^ is the pseudoinverse of D and can be obtained by transposing the matrix and replacing the nonzero values with their multiplicative inverses. That this matrix satisfies the above requirement is directly verified observing that AA^+=UU^* and A^+ A=VV^*, which are the projections onto image and support of A, respectively.


Properties


Existence and uniqueness

As discussed above, for any matrix there is one and only one pseudoinverse . A matrix satisfying only the first of the conditions given above, namely A A^+ A = A, is known as a generalized inverse. If the matrix also satisfies the second condition, namely A^+ A A^+ = 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 has real entries, then so does . * If is
invertible In mathematics, the concept of an inverse element generalises the concepts of opposite () and reciprocal () of numbers. Given an operation denoted here , and an identity element denoted , if , one says that is a left inverse of , and that ...
, its pseudoinverse is its inverse. That is, A^+ = A^.. * The pseudoinverse of the pseudoinverse is the original matrix: \left(A^+\right)^+ = A. * Pseudoinversion commutes with transposition, complex conjugation, and taking the conjugate transpose: \left(A^\mathsf\right)^+ = \left(A^+\right)^\mathsf, \quad \left(\overline\right)^+ = \overline, \quad \left(A^*\right)^+ = \left(A^+\right)^* . * The pseudoinverse of a scalar multiple of is the reciprocal multiple of :\left(\alpha A\right)^+ = \alpha^ A^+ for ; otherwise, \left(0 A\right)^+ = 0 A^+ = 0 A^\mathsf, or 0^+=0^\mathsf. * The kernel and image of the pseudoinverse coincide with those of the conjugate transpose: \ker\left(A^+\right) = \ker\left(A^*\right) and \operatorname\left(A^+\right) = \operatorname\left(A^*\right).


Identities

The following identity formula can be used to cancel or expand certain subexpressions involving pseudoinverses: A = AA^*A^ = A^A^*A. Equivalently, substituting A^+ for A gives A^+ =A^+A^A^* = A^*A^A^+, while substituting A^* for A gives A^* =A^*AA^+=A^+AA^*.


Reduction to Hermitian case

The computation of the pseudoinverse is reducible to its construction in the Hermitian case. This is possible through the equivalences: A^+ = \left(A^*A\right)^+ A^*, A^+ = A^* \left(A A^*\right)^+, as and are Hermitian.


Pseudoinverse of products

The equality does not hold in general. Rather, suppose . Then the following are equivalent: # (AB)^+ = B^+ A^+ # A^+ A BB^* A^* = BB^* A^* and BB^+ A^* A B = A^* A B # \left(A^+ A BB^*\right)^* = A^+ A BB^* and \left(A^* A BB^+\right)^* = A^* A BB^+ # A^+ A BB^* A^* A BB^+ = BB^* A^* A # A^+ A B = B (AB)^+ AB and BB^+ A^* = A^* A B (AB)^+. The following are sufficient conditions for : # has orthonormal columns (then A^*A = A^+ A = I_n), or # has orthonormal rows (then BB^* = BB^+ = I_n), or # has linearly independent columns (then A^+ A = I ) and has linearly independent rows (then BB^+ = I),   or # B = A^*, or # B = A^+. The following is a necessary condition for : # (A^+ A) (BB^+) = (BB^+) (A^+ A) The fourth sufficient condition yields the equalities \begin \left(A A^*\right)^+ &= A^ A^+, \\ \left(A^* A\right)^+ &= A^+ A^. \end Here is a counterexample where : \Biggl( \begin 1 & 1 \\ 0 & 0 \end \begin 0 & 0 \\ 1 & 1 \end \Biggr)^+ = \begin 1 & 1 \\ 0 & 0 \end^+ = \begin \tfrac12 & 0 \\ \tfrac12 & 0 \end \quad \neq \quad \begin \tfrac14 & 0 \\ \tfrac14 & 0 \end = \begin 0 & \tfrac12 \\ 0 & \tfrac12 \end \begin \tfrac12 & 0 \\ \tfrac12 & 0 \end = \begin 0 & 0 \\ 1 & 1 \end^+ \begin 1 & 1 \\ 0 & 0 \end^+


Projectors

P = A A^+ and Q = A^+A are orthogonal projection operators, that is, they are Hermitian (P = P^*, Q = Q^*) and idempotent (P^2 = P and Q^2 = Q). The following hold: * PA = AQ = A and A^+ P = QA^+ = A^+ * is the orthogonal projector onto the range of (which equals the
orthogonal complement In the mathematical fields of linear algebra and functional analysis, the orthogonal complement of a subspace W of a vector space V equipped with a bilinear form B is the set W^\perp of all vectors in V that are orthogonal to every vector in W. I ...
of the kernel of ). * is the orthogonal projector onto the range of (which equals the orthogonal complement of the kernel of ). * I - Q = I - A^+A is the orthogonal projector onto the kernel of . * I - P = I - A A^+ is the orthogonal projector onto the kernel of . The last two properties imply the following identities: * A\,\ \left(I - A^+ A\right)= \left(I - A A^+\right)A\ \ = 0 * A^*\left(I - A A^+\right) = \left(I - A^+A\right)A^* = 0 Another property is the following: if is Hermitian and idempotent (true if and only if it represents an orthogonal projection), then, for any matrix the following equation holds: A(BA)^+ = (BA)^+ This can be proven by defining matrices C = BA, D = A(BA)^+, and checking that is indeed a pseudoinverse for by verifying that the defining properties of the pseudoinverse hold, when is Hermitian and idempotent. From the last property it follows that, if is Hermitian and idempotent, for any matrix (AB)^+A = (AB)^+ Finally, if is an orthogonal projection matrix, then its pseudoinverse trivially coincides with the matrix itself, that is, A^+ = A.


Geometric construction

If we view the matrix as a linear map over the field then can be decomposed as follows. We write for the
direct sum The direct sum is an operation between structures in abstract algebra, a branch of mathematics. It is defined differently but analogously for different kinds of structures. As an example, the direct sum of two abelian groups A and B is anothe ...
, for the
orthogonal complement In the mathematical fields of linear algebra and functional analysis, the orthogonal complement of a subspace W of a vector space V equipped with a bilinear form B is the set W^\perp of all vectors in V that are orthogonal to every vector in W. I ...
, for the kernel of a map, and for the image of a map. Notice that \mathbb^n = \left(\ker A\right)^\perp \oplus \ker A and \mathbb^m = \operatorname A \oplus \left(\operatorname A\right)^\perp. The restriction A: \left(\ker A\right)^\perp \to \operatorname A is then an isomorphism. This implies that on is the inverse of this isomorphism, and is zero on \left(\operatorname A\right)^\perp . In other words: To find for given in , first project orthogonally onto the range of , finding a point in the range. Then form , that is, find those vectors in that sends to . This will be an affine subspace of parallel to the kernel of . The element of this subspace that has the smallest length (that is, is closest to the origin) is the answer we are looking for. It can be found by taking an arbitrary member of and projecting it orthogonally onto the orthogonal complement of the kernel of . This description is closely related to the minimum-norm solution to a linear system.


Limit relations

The pseudoinverse are limits: A^+ = \lim_ \left(A^* A + \delta I\right)^ A^* = \lim_ A^* \left(A A^* + \delta I\right)^ (see Tikhonov regularization). These limits exist even if or do not exist.


Continuity

In contrast to ordinary matrix inversion, the process of taking pseudoinverses is not continuous: if the sequence converges to the matrix (in the maximum norm or Frobenius norm, say), then need not converge to . However, if all the matrices have the same rank as , will converge to .


Derivative

Let x \mapsto A(x) be a real-valued differentiable matrix function with constant rank in a neighborhood of a point . The derivative of x \mapsto A^+(x) at x_0 may be calculated in terms of the derivative of A at x_0: \left.\frac\_ A^+ = -A^+ \left( \frac \right) A^+ ~+~ A^+ A^ \left(\frac \right) \left(I - A A^+\right) ~+~ \left(I - A^+ A\right) \left(\frac \right) A^ A^+, where the functions A, A^+ and derivatives on the right side are evaluated at x_0 (that is, A := A(x_0), A^+ := A^+(x_0), etc.). For a complex matrix, the transpose is replaced with the conjugate transpose. For a real-valued symmetric matrix, the Magnus-Neudecker derivative is established.


Examples

Since for invertible matrices the pseudoinverse equals the usual inverse, only examples of non-invertible matrices are considered below. * For A = \begin 0 & 0 \\ 0 & 0 \end, the pseudoinverse is A^+ = \begin 0 & 0 \\ 0 & 0 \end. The uniqueness of this pseudoinverse can be seen from the requirement A^+ = A^+ A A^+, since multiplication by a zero matrix would always produce a zero matrix. * For A = \begin 1 & 0 \\ 1 & 0 \end, the pseudoinverse is A^+ = \begin \frac & \frac \\ 0 & 0 \end. : Indeed, A\,A^+ = \begin \frac & \frac \\ \frac & \frac \end, and thus A\,A^+ A = \begin 1 & 0 \\ 1 & 0\end = A. Similarly, A^+A = \begin 1 & 0 \\ 0 & 0 \end, and thus A^+A\,A^+ = \begin \frac & \frac \\ 0 & 0 \end = A^+. : Note that is neither injective nor surjective, and thus the pseudoinverse cannot be computed via A^+ = \left(A^* A\right)^ A^* nor A^+ = A^* \left( A A^*\right)^, as A^* A and A A^* are both singular, and furthermore A^+ is neither a left nor a right inverse. : Nonetheless, the pseudoinverse can be computed via SVD observing that A=\sqrt2 \left(\frac\right) \mathbf e_1^*, and thus A^+=\frac \,\mathbf e_1 \left(\frac\right)^*. * For A = \begin 1 & 0 \\ -1 & 0 \end, A^+ = \begin \frac & -\frac \\ 0 & 0 \end. * For A = \begin 1 & 0 \\ 2 & 0 \end, A^+ = \begin \frac & \frac \\ 0 & 0 \end. The denominators are here 5 = 1^2 + 2^2. * For A = \begin 1 & 1 \\ 1 & 1 \end, A^+ = \begin \frac & \frac \\ \frac & \frac \end. * For A = \begin 1 & 0 \\ 0 & 1 \\ 0 & 1 \end, the pseudoinverse is A^+ = \begin 1 & 0 & 0 \\ 0 & \frac & \frac \end. : For this matrix, the left inverse exists and thus equals A^+, indeed, A^+A = \begin 1 & 0 \\ 0 & 1 \end.


Special cases


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 is zero if is zero and the reciprocal of otherwise: x^+ = \begin 0, & \mboxx = 0; \\ x^, & \mbox. \end


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: \vec^+ = \begin \vec^\mathsf, & \text \vec = \vec; \\ pt\dfrac, & \text. \end


Diagonal matrices

The pseudoinverse of a squared diagonal matrix is obtained by taking the reciprocal of the nonzero diagonal elements. Formally, if D is a squared diagonal matrix with D=\tilde D\oplus \mathbf 0_ and \tilde D>0, then D^+=\tilde D^\oplus \mathbf 0_. More generally, if A is any m\times n rectangular matrix whose only nonzero elements are on the diagonal, meaning A_=\delta_ a_i, a_i\in\mathbb K, then A^+ is a n\times m rectangular matrix whose diagonal elements are the reciprocal of the original ones, that is, A_\neq 0\implies A^+_=\frac.


Linearly independent columns

If the rank of is identical to the number of columns, , (for ,) there are
linearly independent In the theory of vector spaces, a set of vectors is said to be if there exists no nontrivial linear combination of the vectors that equals the zero vector. If such a linear combination exists, then the vectors are said to be . These concep ...
columns, and is invertible. In this case, an explicit formula is: A^+ = \left(A^*A\right)^A^*. It follows that is then a left inverse of :   A^+ A = I_n.


Linearly independent rows

If the rank of is identical to the number of rows, , (for ,) there are
linearly independent In the theory of vector spaces, a set of vectors is said to be if there exists no nontrivial linear combination of the vectors that equals the zero vector. If such a linear combination exists, then the vectors are said to be . These concep ...
rows, and is invertible. In this case, an explicit formula is: A^+ = A^*\left(A A^*\right)^. It follows that is a right inverse of :   A A^+ = I_m.


Orthonormal columns or rows

This is a special case of either full column rank or full row rank (treated above). If has orthonormal columns (A^*A = I_n) or orthonormal rows (A A^* = I_m), then: A^+ = A^* .


Normal matrices

If 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 commuting with its transpose implies that it commutes with its pseudoinverse.


EP matrices

A (square) matrix is said to be an EP matrix if it commutes with its pseudoinverse. In such cases (and only in such cases), it is possible to obtain the pseudoinverse as a polynomial in . A polynomial p(t) such that A^+=p(A) can be easily obtained from the characteristic polynomial of or, more generally, from any annihilating polynomial of .


Orthogonal projection matrices

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


Circulant matrices

For a circulant matrix , the singular value decomposition is given by the
Fourier transform In mathematics, the Fourier transform (FT) is an integral transform that takes a function as input then outputs another function that describes the extent to which various frequencies are present in the original function. The output of the tr ...
, that is, the singular values are the Fourier coefficients. Let be the Discrete Fourier Transform (DFT) matrix; then \begin C &= \mathcal\cdot\Sigma\cdot\mathcal^*,\\ C^+ &= \mathcal\cdot\Sigma^+\cdot\mathcal^*. \end


Construction


Rank decomposition

Let denote the rank of . Then can be (rank) decomposed as A = BC where and are of rank . Then A^+ = C^+B^+ = C^*\left(CC^*\right)^\left(B^*B\right)^B^*.


The QR method

For \mathbb \in \ computing the product or 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 may be used instead. Consider the case when is of full column rank, so that A^+ = \left(A^*A\right)^A^*. Then the Cholesky decomposition A^*A = R^*R, where 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, A^+ = \left(A^*A\right)^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 explicitly, by alternatively using the QR decomposition of A = Q R, where Q has orthonormal columns, Q^*Q = I, and is upper triangular. Then A^*A\, =\, (Q R)^*(Q R) \,=\, R^*Q^*Q R \,=\, R^*R , so is the Cholesky factor of . The case of full row rank is treated similarly by using the formula A^+ = A^*\left(A A^*\right)^ and using a similar argument, swapping the roles of and .


Using polynomials in matrices

For an arbitrary , one has that A^*A is normal and, as a consequence, an EP matrix. One can then find a polynomial p(t) such that (A^*A)^+=p(A^*A). In this case one has that the pseudoinverse of is given by A^+= p(A^*A)A^*= A^*p(AA^*).


Singular value decomposition (SVD)

A computationally simple and accurate way to compute the pseudoinverse is by using the
singular value decomposition In linear algebra, the singular value decomposition (SVD) is a Matrix decomposition, factorization of a real number, real or complex number, complex matrix (mathematics), matrix into a rotation, followed by a rescaling followed by another rota ...
.Linear Systems & Pseudo-Inverse
/ref> If A = U\Sigma V^* is the singular value decomposition of , then A^+ = V\Sigma^+ U^*. For a rectangular diagonal matrix such as , we get the pseudoinverse by taking the reciprocal of each non-zero element on the diagonal, leaving the zeros in place. 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 MATLAB (an abbreviation of "MATrix LABoratory") is a proprietary multi-paradigm programming language and numeric computing environment developed by MathWorks. MATLAB allows matrix manipulations, plotting of functions and data, implementat ...
or
GNU Octave GNU Octave is a scientific programming language for scientific computing and numerical computation. Octave helps in solving linear and nonlinear problems numerically, and for performing other numerical experiments using a language that is mostly ...
function , the tolerance is taken to be , 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 has a singular value 0 (a diagonal entry of the matrix above), then modifying 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 A_ = 2A_i - A_i A A_i, which is sometimes referred to as hyper-power sequence. This recursion produces a sequence converging quadratically to the pseudoinverse of if it is started with an appropriate satisfying A_0 A = \left(A_0 A\right)^*. The choice A_0 = \alpha A^* (where 0 < \alpha < 2/\sigma^2_1(A), with denoting the largest singular value of ) 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 enters the region of quadratic convergence. However, if started with already close to the Moore–Penrose inverse and A_0 A = \left(A_0 A\right)^*, for example A_0 := \left(A^* A + \delta I\right)^ A^*, convergence is fast (quadratic).


Updating the pseudoinverse

For the cases where has full row or column rank, and the inverse of the correlation matrix ( for with full row rank or for full column rank) is already known, the pseudoinverse for matrices related to 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. 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.


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 Parallel computing is a type of computing, computation in which many calculations or Process (computing), processes are carried out simultaneously. Large problems can often be divided into smaller ones, which can then be solved at the same time. ...
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. 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) Julia is a high-level programming language, high-level, general-purpose programming language, general-purpose dynamic programming language, dynamic programming language, designed to be fast and productive, for e.g. data science, artificial intel ...
, the LinearAlgebra package of the standard library provides an implementation of the Moore–Penrose inverse pinv() implemented via singular-value decomposition.


Applications


Linear least-squares

The pseudoinverse provides a least squares solution to a
system of linear equations In mathematics, a system of linear equations (or linear system) is a collection of two or more linear equations involving the same variable (math), variables. For example, : \begin 3x+2y-z=1\\ 2x-2y+4z=-2\\ -x+\fracy-z=0 \end is a system of th ...
. For , given a system of linear equations A x = b, in general, a vector 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 b is in the image of A, and is unique if and only if A is injective. The pseudoinverse solves the "least-squares" problem as follows: * , we have \left\, Ax - b\right\, _2 \ge \left\, Az - b\right\, _2 where z = A^+b and \, \cdot\, _2 denotes the
Euclidean norm Euclidean space is the fundamental space of geometry, intended to represent physical space. Originally, in Euclid's ''Elements'', it was the three-dimensional space of Euclidean geometry, but in modern mathematics there are ''Euclidean spaces'' ...
. This weak inequality holds with equality if and only if x = A^+b + \left(I - A^+A\right)w for any vector ; this provides an infinitude of minimizing solutions unless has full column rank, in which case is a zero matrix. The solution with minimum Euclidean norm is This result is easily extended to systems with multiple right-hand sides, when the Euclidean norm is replaced by the Frobenius norm. Let . * , we have \, AX - B\, _ \ge \, AZ -B\, _ where Z = A^+B and \, \cdot\, _ denotes the
Frobenius norm In the field of mathematics, norms are defined for elements within a vector space. Specifically, when the vector space comprises matrices, such norms are referred to as matrix norms. Matrix norms differ from vector norms in that they must also ...
.


Obtaining all solutions of a linear system

If the linear system A x = b has any solutions, they are all given by x = A^+ b + \left - A^+ A\right for arbitrary vector . Solution(s) exist if and only if A A^+ b = b. If the latter holds, then the solution is unique if and only if has full column rank, in which case is a zero matrix. If solutions exist but 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 Ax = b, with non-unique solutions (such as under-determined systems), the pseudoinverse may be used to construct the solution of minimum
Euclidean norm Euclidean space is the fundamental space of geometry, intended to represent physical space. Originally, in Euclid's ''Elements'', it was the three-dimensional space of Euclidean geometry, but in modern mathematics there are ''Euclidean spaces'' ...
\, x\, _2 among all solutions. * If Ax = b is satisfiable, the vector z = A^+b is a solution, and satisfies \, z\, _2 \le \, 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 . * If AX = B is satisfiable, the matrix Z = A^+B is a solution, and satisfies \, Z\, _ \le \, X\, _ for all solutions.


Condition number

Using the pseudoinverse and a
matrix norm In the field of mathematics, norms are defined for elements within a vector space. Specifically, when the vector space comprises matrices, such norms are referred to as matrix norms. Matrix norms differ from vector norms in that they must also ...
, one can define a condition number for any matrix: \mbox(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 can lead to huge errors in the entries of the solution.


Generalizations

The weighted pseudoinverse generalizes the Moore-Penrose inverse between metric spaces with weight matrices in the domain and range. These weights are the identity for the standard Moore-Penrose inverse, which assumes an orthonormal basis in both spaces. In order to solve more general least-squares problems, one can define Moore–Penrose inverses for all continuous linear operators between two
Hilbert space In mathematics, a Hilbert space is a real number, real or complex number, complex inner product space that is also a complete metric space with respect to the metric induced by the inner product. It generalizes the notion of Euclidean space. The ...
s and , using the same four conditions as in our definition above. It turns out that not every continuous linear operator has a continuous linear pseudoinverse in this sense. Those that do are precisely the ones whose range is closed in . A notion of pseudoinverse exists for matrices over an arbitrary field equipped with an arbitrary involutive
automorphism In mathematics, an automorphism is an isomorphism from a mathematical object to itself. It is, in some sense, a symmetry of the object, and a way of mapping the object to itself while preserving all of its structure. The set of all automorphism ...
. In this more general setting, a given matrix doesn't always have a pseudoinverse. The necessary and sufficient condition for a pseudoinverse to exist is that \operatorname(A) = \operatorname\left(A^* A\right) = \operatorname\left(A A^*\right), where A^* denotes the result of applying the involution operation to the transpose of A. When it does exist, it is unique. Example: Consider the field of complex numbers equipped with the identity involution (as opposed to the involution considered elsewhere in the article); do there exist matrices that fail to have pseudoinverses in this sense? Consider the matrix A = \begin1 & i\end^\mathsf. Observe that \operatorname\left(A A^\mathsf\right) = 1 while \operatorname\left(A^\mathsf A\right) = 0. So this matrix doesn't have a pseudoinverse in this sense. In
abstract algebra In mathematics, more specifically algebra, abstract algebra or modern algebra is the study of algebraic structures, which are set (mathematics), sets with specific operation (mathematics), operations acting on their elements. Algebraic structur ...
, a Moore–Penrose inverse may be defined on a *-regular semigroup. This abstract definition coincides with the one in linear algebra.


See also

* Drazin inverse * Hat matrix *
Inverse element In mathematics, the concept of an inverse element generalises the concepts of opposite () and reciprocal () of numbers. Given an operation denoted here , and an identity element denoted , if , one says that is a left inverse of , and that ...
* Linear least squares (mathematics) * Pseudo-determinant *
Von Neumann regular ring In mathematics, a von Neumann regular ring is a ring ''R'' (associative, with 1, not necessarily commutative) such that for every element ''a'' in ''R'' there exists an ''x'' in ''R'' with . One may think of ''x'' as a "weak inverse" of the eleme ...


Notes


References

* * * *


External links

*
Interactive program & tutorial of Moore–Penrose Pseudoinverse
* * *
The Moore–Penrose Pseudoinverse. A Tutorial Review of the Theory


{{DEFAULTSORT:Moore-Penrose inverse Matrix theory Singular value decomposition Numerical linear algebra Roger Penrose