Non-standard matrix inversion was introduced into the statistics literature by Koopman(1997). Non-standard analysis was introduced into mathematics in the 1970's using results from ``model theory'' to embed the real numbers within a framework that includes ``infinite'' and ``infinitesimal'' numbers. It was hoped that this would allow simpler descriptions of limit calculations in real analysis, but never really caught on. See Blass(1978) for a brief survey of the ideas behind this. The goal is to compute an exact (limit) inverse of an input (symmetric) matrix of the form

\begin{equation} \mathbf{A}+\kappa \mathbf{B}\text{ as }\kappa \rightarrow \infty \end{equation}

The result will be the matrix of the form

\begin{equation} {\mathbf{C}} + {\mathbf{D}}\kappa ^{ - 1} + {\mathbf{E}}\kappa ^{ - 2} \label{eq:FormalExpansion} \end{equation}

Note that the \(\kappa ^{ - 2}\) term isn't needed in many applications. With the help of this expansion, it's possible to compute exact limits as \(\kappa \rightarrow \infty\) for calculations like:

\begin{equation} \left( {{\mathbf{F}} + {\mathbf{G}}\kappa } \right)({\mathbf{A}} + {\mathbf{B}}\kappa )^{ - 1} \approx \left( {{\mathbf{F}} + {\mathbf{G}}\kappa } \right)\left( {{\mathbf{C}} + {\mathbf{D}}\kappa ^{ - 1} } \right) \to {\mathbf{FC}} + {\mathbf{GD}} \end{equation}

You might think that you could just use a guess matrix of the form \eqref{eq:FormalExpansion} multiply, match terms and be done quickly. Unfortunately, because matrices generally don't commute, it's very easy to get a left inverse or a right inverse, but not a true inverse. The matching terms idea is correct, but it has to be done carefully. We will start with a case that's more manageable. Note that this derivation is simpler than that in Koopman's paper. Here, the "infinite" matrix is isolated into an identity block, and we look only at the expansion through \(\kappa ^{ - 1}\).

\begin{align} \left( \left[ \begin{array}{cc} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{12}^{\prime } & \mathbf{A}_{22}% \end{array}% \right] +\kappa \left[ \begin{array}{cc} \mathbf{I}_{r} & 0 \\ 0 & 0% \end{array}% \right] \right) &&\times \label{eq:ExactSimple} \\ \left( \left[ \begin{array}{cc} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime } & \mathbf{C}_{22}% \end{array}% \right] +\kappa ^{-1}\left[ \begin{array}{cc} \mathbf{D}_{11} & \mathbf{D}_{12} \\ \mathbf{D}_{12}^{\prime } & \mathbf{D}_{22}% \end{array}% \right] \right) &=&\left[ \begin{array}{cc} \mathbf{I}_{r} & 0 \\ 0 & \mathbf{I}_{n-r}% \end{array}% \right] +\kappa ^{-1}Rem \notag \end{align}

Since the \(\kappa \) term in the product has to zero out, \(\mathbf{C}_{11}=0\) and \(\mathbf{C}_{12}=0\). Using that and collecting the \(O(1)\) terms, we get

\begin{equation} \left( \left[ \begin{array}{cc} \mathbf{D}_{11} & \mathbf{A}_{12}\mathbf{C}_{22}+\mathbf{D}_{12} \\ 0 & \mathbf{A}_{22}\mathbf{C}_{22}% \end{array}% \right] \right) =\left[ \begin{array}{cc} \mathbf{I}_{r} & 0 \\ 0 & \mathbf{I}_{n-r}% \end{array}% \right] \end{equation}

So we get \(\mathbf{D}_{11}=\mathbf{I}_{r},\mathbf{C}_{22}=\mathbf{A}_{22}^{-1}\) and \(\mathbf{D}_{12}=-\mathbf{A}_{12}\mathbf{A}_{22}^{-1}\). \(\mathbf{D}_{22}\) is arbitrary; it just ends up in the remainder, so for simplicity we can make it zero. If we check the reverse multiplication

\begin{equation} \left( \left[ \begin{array}{cc} 0 & 0 \\ 0 & \mathbf{C}_{22}% \end{array}% \right] +\kappa ^{-1}\left[ \begin{array}{cc} \mathbf{I}_{r} & \mathbf{D}_{12} \\ \mathbf{D}_{12}^{\prime } & 0% \end{array}% \right] \right) \left( \left[ \begin{array}{cc} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{12}^{\prime } & \mathbf{A}_{22}% \end{array}% \right] +\kappa \left[ \begin{array}{cc} \mathbf{I}_{r} & 0 \\ 0 & 0% \end{array}% \right] \right) \end{equation}

we can verify that it also will be the identity + a term in \(\kappa ^{-1}\).

 

In general, we won't have an input matrix with the nice form in \eqref{eq:ExactSimple}. However, for a p.s.d. symmetric matrix \(\mathbf{B}\), we can always find a non-singular matrix \(\mathbf{T}\) such that \(\mathbf{T}\mathbf{B}\mathbf{T}'\) has that form. (The simplest to compute is based upon a modified version of the Cholesky factorization.) So, in general, we can compute the inverse with:

\begin{equation} \left( \mathbf{A}+\kappa \mathbf{B}\right) ^{-1}=\mathbf{T'}\left( \mathbf{T}\mathbf{A}\mathbf{T'}+\kappa \mathbf{T}\mathbf{B}\mathbf{T'}\right) ^{-1}\mathbf{T} \label{eq:GeneralExactInverse} \end{equation}

For an input matrix pair \(\left\{ \mathbf{A},\mathbf{B}\right\}\), this will produce an output pair \(\left\{\mathbf{C},\mathbf{D}\right\} \).

 

To derive the more accurate expansion in the case with the simpler form of \(\mathbf{B}\), our test inverse is

\begin{equation} \left( \left[ \begin{array}{cc} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime } & \mathbf{C}_{22}% \end{array}% \right] +\kappa ^{-1}\left[ \begin{array}{cc} \mathbf{D}_{11} & \mathbf{D}_{12} \\ \mathbf{D}_{12}^{\prime } & \mathbf{D}_{22}% \end{array}% \right] +\kappa ^{-2}\left[ \begin{array}{cc} \mathbf{E}_{11} & \mathbf{E}_{12} \\ \mathbf{E}_{12}^{\prime } & \mathbf{E}_{22}% \end{array}% \right] \right) \end{equation}

Everything from above goes through, except we now can't make \(\mathbf{D}_{22} \) arbitrary. When we multiply out, the \(\kappa^{-1}\) terms are

\begin{equation} \left[ \begin{array}{cc} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{12}^{\prime } & \mathbf{A}_{22}% \end{array}% \right] \left[ \begin{array}{cc} \mathbf{D}_{11} & \mathbf{D}_{12} \\ \mathbf{D}_{12}^{\prime } & \mathbf{D}_{22}% \end{array}% \right] +\left[ \begin{array}{cc} \mathbf{I}_{r} & 0 \\ 0 & 0% \end{array}% \right] \left[ \begin{array}{cc} \mathbf{E}_{11} & \mathbf{E}_{12} \\ \mathbf{E}_{12}^{\prime } & \mathbf{E}_{22}% \end{array}% \right] \end{equation}

or

\begin{equation} \left[ \begin{array}{cc} \mathbf{A}_{11}\mathbf{D}_{11}+\mathbf{A}_{12}\mathbf{D}_{12}^{\prime } & \mathbf{A}_{11}\mathbf{D}_{12}+\mathbf{A}_{12}\mathbf{D}_{22} \\ \mathbf{A}_{12}^{\prime }\mathbf{D}_{11}+\mathbf{A}_{22}\mathbf{D}_{12}^{\prime } & \mathbf{A}_{12}^{\prime }\mathbf{D}_{12}+\mathbf{A}_{22}\mathbf{D}_{22} \end{array}% \right] +\left[ \begin{array}{cc} \mathbf{E}_{11} & \mathbf{E}_{12} \\ 0 & 0% \end{array}% \right] \end{equation}

The bottom left element in the \(\mathbf{AD}\) matrix is zero because of the first order solution. Since \(\mathbf{D}_{22}\) was arbitrary from before, we can now solve for it as

\begin{equation} {{\bf{D}}_{22}} = - {\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}{{\bf{D}}_{12}} = {\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}{{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}. \end{equation}

With that, we also get

\begin{equation} {{\bf{E}}_{11}} = -{{{\bf{A}}_{11}}{{\bf{D}}_{11}} - {{\bf{A}}_{12}}{\bf{D}_{12}^{\prime }}} = - {{\bf{A}}_{11}} + {{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }} \end{equation}

and

\begin{equation} {{\bf{E}}_{12}} = -{{{\bf{A}}_{11}}{{\bf{D}}_{12}} - {{\bf{A}}_{12}}{{\bf{D}}_{22}}} = \left( {{{\bf{A}}_{11}} - {{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}} \right){{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1} \end{equation}

\(\mathbf{E}_{22}\) is now arbitrary. In terms of the input matrix, this is:

\begin{align} \left[ {\begin{array}{*{20}{c}} 0 & 0 \\ 0 & {{\bf{A}}_{22}^{ - 1}} \\ \end{array}} \right] + {\kappa ^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\bf{I}} & {{{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}} \\ {{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}} & {{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}{{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}} \notag \\ \end{array}} \right] + \\ {\kappa ^{ - 2}}\left[ {\begin{array}{*{20}{c}} { - {{\bf{A}}_{11}} + {{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}} & {\left( {{{\bf{A}}_{11}} - {{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}} \right){{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}} \\ {{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime}}\left( {{{\bf{A}}_{11}} - {{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}{\bf{A}_{12}^{\prime }}} \right)} & 0 \\ \end{array}} \right] \label{eq:InverseFinalForm} \end{align}

The extension of this to the more general \(\mathbf{B}\) matrix is as before. In the typical situation, the transforming \(\mathbf{T}\) matrix will take the form

\begin{equation} \left[ {\begin{array}{*{20}{c}} {{{\bf{T}}_1}} \\ {{{\bf{T}}_2}} \\ \end{array}} \right] \end{equation}

where \({{{\bf{T}}_1}}\) is \(r \times n\) and \({{{\bf{T}}_2}}\) is \((n-r) \times n\). If that's the case, then the component matrices in the formula \eqref{eq:InverseFinalForm} are \({{\bf{A}}_{11}} = {{\bf{T}}_1}{\bf{A}}{{{\bf{T'}}}_1}\), \({{\bf{A}}_{12}} = {{\bf{T}}_1}{\bf{AT'}}\) and \({{\bf{A}}_{22}} = {{\bf{T}}_2}{\bf{A}}{{{\bf{T'}}}_2}\). The expanded inverse will be \eqref{eq:InverseFinalForm} premultiplied by \(\mathbf{T}'\) and postmultiplied by \(\mathbf{T}\).

 

Copyright © 2025 Thomas A. Doan