Moore–Penrose pseudoinverse

In mathematics, and in particular linear algebra, a pseudoinverse A+ of a matrix A is a generalization of the inverse matrix.[1] The most widely known type of matrix pseudoinverse is the Moore–Penrose pseudoinverse, which was independently described by E. H. Moore[2] in 1920, Arne Bjerhammar[3] in 1951 and Roger Penrose[4] in 1955. Earlier, Fredholm had introduced the concept of a pseudoinverse of integral operators in 1903. When referring to a matrix, the term pseudoinverse, without further specification, is often used to indicate the Moore–Penrose pseudoinverse. The term generalized inverse is sometimes used as a synonym for pseudoinverse.

A common use of the pseudoinverse is to compute a 'best fit' (least squares) solution to a system of linear equations that lacks a unique 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 and unique for all matrices whose entries are real or complex numbers. It can be computed using the singular value decomposition.

Notation

In the following discussion, the following conventions are adopted.

Definition

For A\in \mathrm {M} (m,n;K), a pseudoinverse of A is defined as a matrix A^{+}\in \mathrm {M} (n,m;K) satisfying all of the following four criteria:[4][5]

  1. AA^{+}A=A\,\!       (AA+ need not be the general identity matrix, but it maps all column vectors of A to themselves);
  2. A^{+}AA^{+}=A^{+}\,\!       (A+ is a weak inverse for the multiplicative semigroup);
  3. (AA^{+})^{*}=AA^{+}\,\!       (AA+ is Hermitian); and
  4. (A^{+}A)^{*}=A^{+}A\,\!       (A+A is also Hermitian).

A^{+} exists for any matrix, A, but when the latter has full rank, A^{+} can be expressed as a simple algebraic formula.

In particular, when A has linearly independent columns (and thus matrix A^{*}A is invertible), A^{+} can be computed as:

A^{+}=(A^{*}A)^{-1}A^{*}\,.

This particular pseudoinverse constitutes a left inverse, since, in this case, A^{+}A=I.

When A has linearly independent rows (matrix AA^{*} is invertible), A^{+} can be computed as:

A^{+}=A^{*}(AA^{*})^{-1}\,.

This is a right inverse, as AA^{+}=I.

Properties

Proofs for some of these facts may be found on a separate page here.

Existence and uniqueness

A matrix satisfying the first condition of the definition is known as a generalized inverse. If the matrix also satisfies the second definition, 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

(A^{\mathrm {T} })^{+}=(A^{+})^{\mathrm {T} },~~(\,{\overline {A}}\,)^{+}={\overline {A^{+}}},~~(A^{*})^{+}=(A^{+})^{*}.\,\!
(\alpha A)^{+}=\alpha ^{-1}A^{+}\,\! for \alpha \neq 0.

Identities

The following identities can be used to cancel certain subexpressions or expand expressions involving pseudoinverses. Proofs for these properties can be found in the proofs subpage.

{\begin{array}{lclll}A^{+}&=&A^{+}&A^{+*}&A^{*}\\A^{+}&=&A^{*}&A^{+*}&A^{+}\\A&=&A^{+*}&A^{*}&A\\A&=&A&A^{*}&A^{+*}\\A^{*}&=&A^{*}&A&A^{+}\\A^{*}&=&A^{+}&A&A^{*}\\\end{array}}

Reduction to Hermitian case

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

as A^{*}A and AA^{*} are obviously Hermitian.

Products

If A\in \mathrm {M} (m,n;K),~B\in \mathrm {M} (n,p;K)\, and either,

then

(AB)^{+}\equiv B^{+}A^{+}\,\!.

The last property yields the equivalences:

{\begin{aligned}(AA^{*})^{+}&\equiv A^{+*}A^{+}\\(A^{*}A)^{+}&\equiv A^{+}A^{+*}\end{aligned}}

Projectors

P=AA^{+}\,\! 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:

Geometric construction

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

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

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

Subspaces

{\begin{aligned}\operatorname {ker} (A^{+})&=\operatorname {ker} (A^{*})\\\operatorname {im} (A^{+})&=\operatorname {im} (A^{*})\end{aligned}}

Limit relations

A^{+}=\lim _{\delta \searrow 0}(A^{*}A+\delta I)^{-1}A^{*}=\lim _{\delta \searrow 0}A^{*}(AA^{*}+\delta I)^{-1}
(see Tikhonov regularization). These limits exist even if (AA^{*})^{-1}\,\! or (A^{*}A)^{-1}\,\! do not exist.[5]:263

Continuity

Derivative

The derivative of a real valued pseudoinverse matrix which has constant rank at a point x may be calculated in terms of the derivative of the original matrix:[8]

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

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 x is zero if x is zero and the reciprocal of x otherwise:

x^{+}=\left\{{\begin{matrix}0,&{\mbox{if }}x=0;\\x^{-1},&{\mbox{otherwise}}.\end{matrix}}\right.

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:

x^{+}=\left\{{\begin{matrix}0^{\mathrm {T} },&{\mbox{if }}x=0;\\{x^{*} \over x^{*}x},&{\mbox{otherwise}}.\end{matrix}}\right.

Linearly independent columns

If the columns of A\,\! are linearly independent (so that m\geq n), then A^{*}A\,\! is invertible. In this case, an explicit formula is:[1]

A^{+}=(A^{*}A)^{-1}A^{*}\,\!.

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

Linearly independent rows

If the rows of A\,\! are linearly independent (so that m\leq n), then AA^{*} is invertible. In this case, an explicit formula is:

A^{+}=A^{*}(AA^{*})^{-1}\,\!.

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

Orthonormal columns or rows

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

Circulant matrices

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

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

Construction

Rank decomposition

Let r\leq \min(m,n) denote the rank of A\in \mathrm {M} (m,n;K)\,\!. Then A\,\! can be (rank) decomposed as A=BC\,\! where B\in \mathrm {M} (m,r;K)\,\! and C\in \mathrm {M} (r,n;K)\,\! are of rank r. Then A^{+}=C^{+}B^{+}=C^{*}(CC^{*})^{-1}(B^{*}B)^{-1}B^{*}\,\!.

The QR method

For K=\mathbb {R} \,\! or K=\mathbb {C} \,\! computing the product AA^{*} or 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 A\,\! may be used instead.

Considering the case when A\,\! is of full column rank, so that A^{+}=(A^{*}A)^{-1}A^{*}\,\!. Then the Cholesky decomposition A^{*}A=R^{*}R\,\!, where 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,

A^{+}=(A^{*}A)^{-1}A^{*}\quad \Leftrightarrow \quad (A^{*}A)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 A^{*}A\,\! explicitly, by alternatively using the QR decomposition of A=QR\,\!, where Q\,\,\! has orthonormal columns, Q^{*}Q=I, and R\,\! is upper triangular. Then

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

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

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

Singular value decomposition (SVD)

A computationally simple and accurate way to compute the pseudo inverse is by using the singular value decomposition.[1][5][10] If A=U\Sigma V^{*} is the singular value decomposition of A, then A^{+}=V\Sigma ^{+}U^{*}. For a rectangular diagonal matrix such as \Sigma , we get the pseudo inverse 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, GNU Octave, or NumPy 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 pseudo inverse is not a continuous operation: if the original matrix A has a singular value 0 (a diagonal entry of the matrix \Sigma above), then modifying A slightly may turn this zero into a tiny positive number, thereby affecting the pseudo inverse 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 uses the recursion

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 A if it is started with an appropriate A_{0} satisfying A_{0}A=(A_{0}A)^{*}. The choice A_{0}=\alpha A^{*} (where 0<\alpha <2/\sigma _{1}^{2}(A), with \sigma _{1}(A) denoting the largest singular value of A) [11] 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 A_{i} enters the region of quadratic convergence.[12] However, if started with A_{0} already close to the Moore–Penrose pseudoinverse and A_{0}A=(A_{0}A)^{*}, for example A_{0}:=(A^{*}A+\delta I)^{-1}A^{*}, convergence is fast (quadratic).

Updating the pseudoinverse

For the cases where A has full row or column rank, and the inverse of the correlation matrix (AA^{*} for A with full row rank or A^{*}A for full column rank) is already known, the pseudoinverse for matrices related to 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[13][14] 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.[15][16]

Software libraries

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

Applications

Linear least-squares

The pseudoinverse provides a least squares solution to a system of linear equations.[17] For A\in \mathrm {M} (m,n;K)\,\!, given a system of linear equations

Ax=b,\,

in general, a vector x that solves the system may not exist, or if one does exist, it may not be unique. The pseudoinverse solves the "least-squares" problem as follows:

This result is easily extended to systems with multiple right-hand sides, when the Euclidean norm is replaced by the Frobenius norm. Let B\in \mathrm {M} (m,p;K).

Obtaining all solutions of a linear system

If the linear system

Ax=b\,

has any solutions, they are all given by[19]

x=A^{+}b+[I-A^{+}A]w

for arbitrary vector w. Solution(s) exist if and only if AA^{+}b=b.[19] If the latter holds, then the solution is unique if and only if A has full column rank, in which case [I-A^{+}A] is a zero matrix. If solutions exist but 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. This solution is deeply connected to the Udwadia–Kalaba equation of classical mechanics to forces of constraint that do not obey D'Alembert's principle.

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 \|x\|_{2} among 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 B\in \mathrm {M} (m,p;K)\,\!.

Condition number

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

{\mbox{cond}}(A)=\|A\|\|A^{+}\|.\

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 A can lead to huge errors in the entries of the solution.[20]

Generalizations

In order to solve more general least-squares problems, one can define Moore–Penrose pseudoinverses for all continuous linear operators A : H1H2 between two Hilbert spaces H1 and H2, 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.[20] Those that do are precisely the ones whose range is closed in H2.

In abstract algebra, a Moore–Penrose pseudoinverse may be defined on a *-regular semigroup. This abstract definition coincides with the one in linear algebra.

See also

References

  1. 1 2 3 Ben-Israel, Adi; Thomas N.E. Greville (2003). Generalized Inverses. Springer-Verlag. ISBN 0-387-00293-6.
  2. Moore, E. H. (1920). "On the reciprocal of the general algebraic matrix". Bulletin of the American Mathematical Society 26 (9): 394–395. doi:10.1090/S0002-9904-1920-03322-7.
  3. Bjerhammar, Arne (1951). "Application of calculus of matrices to method of least squares; with special references to geodetic calculations". Trans. Roy. Inst. Tech. Stockholm 49.
  4. 1 2 Penrose, Roger (1955). "A generalized inverse for matrices". Proceedings of the Cambridge Philosophical Society 51: 406–413. doi:10.1017/S0305004100030401.
  5. 1 2 3 4 5 Golub, Gene H.; Charles F. Van Loan (1996). Matrix computations (3rd ed.). Baltimore: Johns Hopkins. pp. 257–258. ISBN 0-8018-5414-8.
  6. 1 2 3 Stoer, Josef; Bulirsch, Roland (2002). Introduction to Numerical Analysis (3rd ed.). Berlin, New York: Springer-Verlag. ISBN 978-0-387-95452-3..
  7. Rakočević, Vladimir (1997). "On continuity of the Moore–Penrose and Drazin inverses" (PDF). Matematički Vesnik 49: 163–172.
  8. http://www.jstor.org/stable/2156365
  9. Stallings, W. T.; Boullion, T. L. (1972). "The Pseudoinverse of an r-Circulant Matrix". Proceedings of the American Mathematical Society 34: 385–388. doi:10.2307/2038377.
  10. Linear Systems & Pseudo-Inverse
  11. Ben-Israel, Adi; Cohen, Dan (1966). "On Iterative Computation of Generalized Inverses and Associated Projections". SIAM Journal on Numerical Analysis 3: 410–419. doi:10.1137/0703035. JSTOR 2949637.pdf
  12. Söderström, Torsten; Stewart, G. W. (1974). "On the Numerical Properties of an Iterative Method for Computing the Moore–Penrose Generalized Inverse". SIAM Journal on Numerical Analysis 11: 61–74. doi:10.1137/0711008. JSTOR 2156431.
  13. Tino Gramß (1992). "Worterkennung mit einem künstlichen neuronalen Netzwerk". Georg-August-Universität zu Göttingen.
  14. , Mohammad Emtiyaz, "Updating Inverse of a Matrix When a Column is Added/Removed"
  15. Meyer, Carl D., Jr. Generalized inverses and ranks of block matrices. SIAM J. Appl. Math. 25 (1973), 597–602
  16. Meyer, Carl D., Jr. Generalized inversion of modified matrices. SIAM J. Appl. Math. 24 (1973), 315–323
  17. Penrose, Roger (1956). "On best approximate solution of linear matrix equations". Proceedings of the Cambridge Philosophical Society 52: 17–19. doi:10.1017/S0305004100030929.
  18. 1 2 Planitz, M., "Inconsistent systems of linear equations", Mathematical Gazette 63, October 1979, 181–185.
  19. 1 2 James, M., "The generalised inverse", Mathematical Gazette 62, June 1978, 109–114.
  20. 1 2 Roland Hagen, Steffen Roch, Bernd Silbermann. C*-algebras and Numerical Analysis, CRC Press, 2001. Section 2.1.2.

External links

This article is issued from Wikipedia - version of the Wednesday, November 25, 2015. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.