NAG CL Interface
F01 (Matop)
Matrix Operations, Including Inversion

1 Scope of the Chapter

2 Background to the Problems

2.1 Matrix Inversion

1. (i)Nonsingular square matrices of order $n$.
   If $A$, a square matrix of order $n$, is nonsingular (has rank $n$), then its inverse $X$ exists and satisfies the equations $AX=XA=I$ (the identity or unit matrix).
   It is worth noting that if $AX-I=R$, so that $R$ is the ‘residual’ matrix, then a bound on the relative error is given by $‖R‖$, i.e.,

   $$\frac{‖X-A^{-1}‖}{‖A^{-1}‖}\le ‖R‖\text{.}$$

2. (ii)General real rectangular matrices.
   A real matrix $A$ has no inverse if it is square ($n$ by $n$) and singular (has rank $<n$), or if it is of shape ($m$ by $n$) with $m\ne n$, but there is a Generalized or Pseudo-inverse $A^{+}$ which satisfies the equations

   $$A{A}^{+}A=A\text{, \hspace{1em}}{A}^{+}A{A}^{+}=A^{+}\text{, \hspace{1em}}{\left(A{A}^{+}\right)}^{\mathrm{T}}=A{A}^{+}\text{, \hspace{1em}}{\left(A^{+}A\right)}^{\mathrm{T}}=A^{+}A$$

   (which of course are also satisfied by the inverse $X$ of $A$ if $A$ is square and nonsingular).

   1. (a)if $m\ge n$ and $\mathrm{rank}\left(A\right)=n$ then $A$ can be factorized using a $\mathit{Q}\mathit{R}$ factorization, given by

      $$A=Q\left(\begin{array}{l}R\\ 0\end{array}\right)\text{,}$$

      where $Q$ is an $m$ by $m$ orthogonal matrix and $R$ is an $n$ by $n$, nonsingular, upper triangular matrix. The pseudo-inverse of $A$ is then given by

      $$A^{+}=R^{-1}{\tilde{Q}}^{\mathrm{T}}\text{,}$$

      where $\tilde{Q}$ consists of the first $n$ columns of $Q$.
   2. (b)if $m\le n$ and $\mathrm{rank}\left(A\right)=m$ then $A$ can be factorized using an RQ factorization, given by

      $$A=\left(R0\right)Q^{\mathrm{T}}$$

      where $Q$ is an $n$ by $n$ orthogonal matrix and $R$ is an $m$ by $m$, nonsingular, upper triangular matrix. The pseudo-inverse of $A$ is then given by

      $$A^{+}=\tilde{Q}{R}^{-1}\text{,}$$

      where $\tilde{Q}$ consists of the first $m$ columns of $Q$.
   3. (c)if $m\ge n$ and $\mathrm{rank}\left(A\right)=r\le n$ then $A$ can be factorized using a $QR$ factorization, with column interchanges, as

      $$A=Q\left(\begin{array}{l}R\\ 0\end{array}\right)P^{\mathrm{T}}\text{,}$$

      where $Q$ is an $m$ by $m$ orthogonal matrix, $R$ is an $r$ by $n$ upper trapezoidal matrix and $P$ is an $n$ by $n$ permutation matrix. The pseudo-inverse of $A$ is then given by

      $$A^{+}=P{R}^{\mathrm{T}}{\left(R{R}^{\mathrm{T}}\right)}^{-1}{\tilde{Q}}^{\mathrm{T}}\text{,}$$

      where $\tilde{Q}$ consists of the first $r$ columns of $Q$.
   4. (d)if $\mathrm{rank}\left(A\right)=r\le k=\min \left(m,n\right)$, then $A$ can be factorized as the singular value decomposition

      $$A=U\Sigma V^{\mathrm{T}}\text{,}$$

      where $U$ is an $m$ by $m$ orthogonal matrix, $V$ is an $n$ by $n$ orthogonal matrix and $\Sigma$ is an $m$ by $n$ diagonal matrix with non-negative diagonal elements $\sigma$. The first $k$ columns of $U$ and $V$ are the left- and right-hand singular vectors of $A$ respectively and the $k$ diagonal elements of $\Sigma$ are the singular values of $A$. $\Sigma$ may be chosen so that

      $${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_k\ge 0$$

      and in this case if $\mathrm{rank}\left(A\right)=r$ then

      $${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0\text{, \hspace{1em}}{\sigma }_{r+1}=\cdots ={\sigma }_k=0\text{.}$$

      If $\tilde{U}$ and $\tilde{V}$ consist of the first $r$ columns of $U$ and $V$ respectively and $\tilde{\Sigma }$ is an $r$ by $r$ diagonal matrix with diagonal elements ${\sigma }_1,{\sigma }_2,\dots ,{\sigma }_r$ then $A$ is given by

      $$A=\tilde{U}\tilde{\Sigma }{\tilde{V}}^{\mathrm{T}}$$

      and the pseudo-inverse of $A$ is given by

      $$A^{+}=\tilde{V}{\tilde{\Sigma }}^{-1}{\tilde{U}}^{\mathrm{T}}\text{.}$$

      Notice that

      $$A^{\mathrm{T}}A=V\left({\Sigma }^{\mathrm{T}}\Sigma \right)V^{\mathrm{T}}$$

      which is the classical eigenvalue (spectral) factorization of $A^{\mathrm{T}}A$.
   5. (e)if $A$ is complex then the above relationships are still true if we use ‘unitary’ in place of ‘orthogonal’ and conjugate transpose in place of transpose. For example, the singular value decomposition of $A$ is

      $$A=U\Sigma V^{\mathrm{H}}\text{,}$$

      where $U$ and $V$ are unitary, $V^{\mathrm{H}}$ the conjugate transpose of $V$ and $\Sigma$ is as in (ii)(d) above.

2.2 Matrix Factorizations

2.3 Matrix Arithmetic and Manipulation

2.4 Matrix Functions

3 Recommendations on Choice and Use of Available Functions

3.1 Matrix Inversion

3.2 Matrix Factorizations

3.3 Matrix Arithmetic and Manipulation

3.3.1 NAG Names and LAPACK Names

3.4 Matrix Functions

4 Decision Trees

Tree 1

Tree 2: Inverse of a real n by n matrix of full rank

Tree 3: Inverse of a complex n by n matrix of full rank

Tree 4: Pseudo-inverses

Tree 5: Matrix functions $f\left(A\right)$ of an n by n real matrix $A$

Tree 6: Matrix functions $f\left(A\right)$ of an n by n complex matrix $A$

5 Functionality Index

6 Auxiliary Functions Associated with Library Function Arguments

7 Withdrawn or Deprecated Functions

8 References