# Inverting a complex matrix

Zhen Dai  ,  Lek-Heng Lim  Computational and Applied Mathematics Initiative,
University of Chicago, Chicago, IL 60637-1514<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Ke Ye  KLMM, Academy of Mathematics and Systems
Science, Chinese Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>
###### Abstract.

We analyze a complex matrix inversion algorithm first proposed by Frobenius,
but largely forgotten:
$(A+iB)^{-1}=(A+BA^{-1}B)^{-1}-iA^{-1}B(A+BA^{-1}B)^{-1}$ when $A$ is
invertible and $(A+iB)^{-1}=B^{-1}A(AB^{-1}A+B)^{-1}-i(AB^{-1}A+B)^{-1}$ when
$B$ is invertible. This may be viewed as an inversion analogue of the
aforementioned Gauss multiplication. We proved that Frobenius inversion is
optimal — it uses the least number of real matrix multiplications and
inversions among all complex matrix inversion algorithms. We also showed that
Frobenius inversion runs faster than the standard method based on LU
decomposition if and only if the ratio of the running time for real matrix
inversion to that for real matrix multiplication is greater than $5/4$. We
corroborate this theoretical result with extensive numerical experiments,
applying Frobenius inversion to evaluate matrix sign function, solve Sylvester
equation, and compute polar decomposition, concluding that for these problems,
Frobenius inversion is more efficient than LU decomposition with nearly no
loss in accuracy.

##  1\. Introduction

We analyze a complex matrix inversion algorithm mentioned in [7, 13, 15, 36,
40, 52, 55, 59, 64]. Given a complex matrix $A+iB\in\mathbb{C}^{n\times n}$
with $A,B\in\mathbb{R}^{n\times n}$, we compute its inverse by

$Z^{-1}=\begin{cases}(A+BA^{-1}B)^{-1}-iA^{-1}B(A+BA^{-1}B)^{-1}&\text{if $A$
is invertible},\\\ B^{-1}A(AB^{-1}A+B)^{-1}-i(AB^{-1}A+B)^{-1}&\text{if $B$ is
invertible}.\end{cases}$

We call this algorithm “Frobenius inversion”. Although this algorithm is
known, its numerical properties have not been studied, as far as we know. In
this paper, we give a thorough analysis of Frobenius inversion. We prove that
Frobenius inversion is optimal in the sense of least number of matrix
multiplications and inversions over $\mathbb{R}$. Then, we show that Frobenius
inversion is faster than MATLAB’s inversion algorithm. In addition, we provide
several applications where Frobenius inversion is more efficient than MATLAB’s
inversion algorithm.

Matrix inversion is not preferred in most numerical linear algebra problems.
For instance, when solving a system of linear equations $Ax=b$, where $A$ is
an $n$ by $n$ invertible matrix, we should not compute $A^{-1}$ first and then
multiplying it to $b$. Instead, we could use a LU factorization approach.
First, we compute a LU factorization of $A=LU$. Then, we solve the systems
$Ly=b$ and $Ux=y$ by forward substitution and backward substitution. This
approach is both faster and more accurate than the one using matrix inversion
[26]. In terms of speed, the LU factorization approach takes $2n^{3}/3$ flops
while the matrix inversion approach takes $2n^{3}$ flops [26]. To compare the
accuracy, let $\widehat{x}_{LU}$ and $\widehat{x}_{\operatorname{inv}}$ be the
computed values of $x$ using LU factorization and matrix inversion
respectively. Let $\widehat{L}$ and $\widehat{U}$ be the computed LU factors
of $A$. As shown in [26],

$|b-A\hat{x}_{\operatorname{inv}}|\leq
n|A||A^{-1}||b|\mathsf{u}+O(\mathsf{u}^{2})\quad\textrm{and}\quad|b-A\hat{x}_{LU}|\leq
3n|\widehat{L}||\widehat{U}||\hat{x}_{LU}|\mathsf{u}+O(\mathsf{u}^{2}),$

where $\mathsf{u}$ is the unit roundoff, $n$ is the dimension of $b$, and
$|\cdot|$ is applied componentwise. Usually,
$\lVert|\widehat{L}||\widehat{U}|\rVert_{\infty}\approx\lVert
A\rVert_{\infty}$ [26]. Thus, $\hat{x}_{LU}$ is likely to be much more
accurate than $\hat{x}_{\operatorname{inv}}$ when $\lVert
x\rVert_{\infty}\ll\lVert|A^{-1}||b|\rVert_{\infty}$ [26].

Although matrix inversion is generally not preferred in numerical linear
algebra, it is necessary to do matrix inversion in certain circumstances. In
some eigenvalue-related problems, matrix inversion is necessary in iterations
of certain algorithms [1, 9, 28, 51]. For instance, Newton’s iteration of
solving Sylvester equation has the following update rule [51]:

$X_{t+1}=\frac{1}{2}(X_{t}+X_{t}^{-1}),\quad t\in\mathbb{N},$

where $X_{t}$ is some matrix related to the inputs of Sylvester equation. We
will discuss this in detail in section 7. In superconductivity computations,
matrix inversion is needed in some numerical integrations [22]. Specifically,
in the KKR-CPA algorithm, we need to compute the KKR matrix inverse and then
integrate it over the first Brillouin zone [22]. In some MIMO radios, matrix
inverse must be applied in hardware [12, 14, 58]. In statistics, sometimes
matrix inverse reveals important statistical properties [4, 41, 43]. For
instance, in linear regression, we want to learn a model $Y=X\hat{\beta}$,
where $X$ is the design matrix, $Y$ are the observed values of the dependent
variable, and $\hat{\beta}$ are the coefficients of least square regression
[43]. Then, the covariance matrix of $\hat{\beta}$ is
$\sigma^{2}\cdot(X^{\scriptscriptstyle\mathsf{T}}X)^{-1}$, where $\sigma^{2}$
is the variance of the dependent variable [42]. Thus, we need to compute the
inverse of $X^{\scriptscriptstyle\mathsf{T}}X$ in order to understand the
statistical properties of $\hat{\beta}$. Moreover, if we want to estimate the
accuracy of the prediction at a new point $x$, the variance of the predicted
value $x^{T}\hat{\beta}$ is
$\sigma^{2}x^{\scriptscriptstyle\mathsf{T}}(X^{\scriptscriptstyle\mathsf{T}}X)^{-1}x$.
In order to build a confidence interval for this new prediction, we need to
compute
$\sigma^{2}x^{\scriptscriptstyle\mathsf{T}}(X^{\scriptscriptstyle\mathsf{T}}X)^{-1}x$.
Since we are interested in predicting the value of $Y$ at many different
points $x$, it is preferable to compute the inverse of
$X^{\scriptscriptstyle\mathsf{T}}X$. In addition, by Cramer-Rao lower bound
[10, 49], the inverse of the Fisher information matrix is an asymptotic lower
bound for the covariance matrix of an unbiased estimator. In some Gaussian
process, this lower bound could be attained [33]. Thus, we need to compute the
inverse of Fisher information matrix in order to understand certain
statistical problems. In this case, matrix inversion is unavoidable since the
final result is the inverse of some matrix. In graph theory, the inverses of
the adjacency matrix and forward adjacency matrix of a graph $G$ reveal
important combinatorial properties of $G$ [44, 46, 47, 63]. For instance, the
inverse of the adjacency matrix of $G$ is the adjacency matrix of the inverse
graph of $G$.

Furthermore, matrix inversion could even be beneficial in certain numerical
linear algebra tasks. We first reproduce an example in [11]. Consider the
scenario where we want to solve a sequence of linear systems
$Ax^{(k)}=b^{(k)}$, for $k=1,\dots,K$, with the same coefficient matrix $A$.
Suppose that $b^{(k+1)}$ depends on $b^{(k)}$ so that we have to solve the
linear systems one by one. Then there are two natural approaches to solve this
problem. One approach is to compute $A^{-1}$ first and then multiply it to
$b^{(k)}$ one at a time. Another approach is to compute a LU factorization of
$A=LU$ and then solve $LUx^{(k)}=b^{(k)}$ by forward substitution and backward
substitution for each $k=1,\dots,K$. When $K\gg n$, both algorithms takes
roughly the same number of flops. However, the matrix inversion approach is
faster than the LU factorization approach in practice due to some data
structure considerations [11]. In terms of accuracy, the forward errors
$\lVert\hat{x}^{(k)}_{\operatorname{inv}}-x^{(k)}\rVert$ and
$\lVert\hat{x}^{(k)}_{LU}-x^{(k)}\rVert$ are close to each other under mild
conditions [12], where $\hat{x}^{(k)}_{\operatorname{inv}}$ and
$\hat{x}^{(k)}_{LU}$ are solutions of $Ax^{(k)}=b^{(k)}$ computed by matrix
inversion and LU factorization respectively for each $k=1,\dots,K$. In
addition, matrix inversion could also be beneficial for rank-one updates.
Suppose that we have computed $B^{-1}$ for some matrix $B$. Then, for any
vectors $u$ and $v$, we can compute the inverse of
$A=B+uv^{\scriptscriptstyle\mathsf{T}}$ by the reinforcement method [17]:

$A^{-1}=B^{-1}-\frac{1}{1+v^{\scriptscriptstyle\mathsf{T}}B^{-1}u}B^{-1}uv^{\scriptscriptstyle\mathsf{T}}B^{-1},$

assuming that $1+v^{\scriptscriptstyle\mathsf{T}}B^{-1}u\neq 0$. Besides this
example, computing the inverse of a matrix $A$ could also be beneficial if $A$
has special structures [18, 32, 35]. For instance, if $A$ is block circulant,
then $A$ is diagonalizable by some $FFT$ type matrices, which makes it fast to
compute the inverse of $A$ [18].

###  1.1. Related work

Frobenius inversion was first discovered by Frobenius and Schur [59]. It was
then derived from Frobenius-Schur’s relation in [15, pp. 217–219]. Frobenius
inversion was also derived in [40, pp. 137–138]. In [7], the authors proposed
a way to invert $A+iB$ by inverting the real matrix $\begin{bmatrix}A&B\\\
-B&A\end{bmatrix}$. In [13], the authors compared several ways to invert a
complex matrix using real matrix inversions with directly inverting a complex
matrix in complex arithmetic. They proposed that inverting a complex matrix in
complex arithmetic is better than inverting a complex matrix using real
inversions. The methods they considered include Frobenius inversion and the
one proposed in [7]. However, in this paper, we show that Frobenius inversion
could actually be advantageous in certain circumstances. In [36], the authors
proposed a way to invert $A+iB$ using real matrix inversions when both $A$ and
$B$ are singular. Their algorithm uses Frobenius inversion as an intermediate
step. In [55], the authors proposed another way to invert $A+iB$ using real
matrix inversions without assuming $A$, $B$, $A+B$ or $A-B$ is invertible. In
[52, 64], Frobenius inversion was derived again.

##  2\. Multiplications in quadratic field extensions

Let $\Bbbk$ be a field and let $\mathbb{F}$ be a quadratic field extension of
$\Bbbk$. Thus there exists some monic irreducible quadratic polynomial
$f(x)\in\Bbbk[x]$ such that

$\mathbb{F}\simeq\Bbbk[x]/(f(x)).$

We explicitly write $f(x)=x^{2}+\beta x+\tau$ for some $\beta,c\in\Bbbk$. Up
to an isomorphism, we are able to write $f$ in a normal form:

  * •

$\operatorname{char}(\Bbbk)\neq 2$: $\beta=0$ and $-\tau$ is not a complete
square in $\Bbbk$.

  * •

$\operatorname{char}(\Bbbk)=2$: either $\beta=0$ and $-\tau$ is not a complete
square in $\Bbbk$, or $\beta=1$ and $x^{2}+x+\tau$ has no solution in $\Bbbk$.

Let $\xi$ be a root of $f(x)$ in $\overline{\Bbbk}$. We also have
$\mathbb{F}\simeq\Bbbk[\xi]$. In particular, an element in $\mathbb{F}$ can be
written as $a_{1}+a_{2}\xi$. The multiplicaiton on $\mathbb{F}$ is given by
the lemma that follows.

######  Lemma 2.1.

Let $\Bbbk,f,\xi$ be as above. Given two elements
$x_{1}+x_{2}\xi,y_{1}+y_{2}\xi$ in $\Bbbk[\xi]$, their product is

$(x_{1}+x_{2}\xi)(y_{1}+y_{2}\xi)=\begin{cases}(x_{1}y_{1}-\tau
x_{2}y_{2})+(x_{1}y_{2}+x_{2}y_{1})\xi,~{}\text{if}~{}f(x)=x^{2}+\tau,\\\
(x_{1}y_{1}-cx_{2}y_{2})+(x_{1}y_{2}+x_{2}y_{1}-x_{2}y_{2})\xi,~{}\text{if}~{}f(x)=x^{2}+x+\tau.\end{cases}$
(1)

######  Proposition 2.2.

Let $\Bbbk,f,\tau,\xi$ be as above. There exists an algorithm for the
multiplication in $\Bbbk[\xi]$, which costs three multiplications in $\Bbbk$.
Moreover, the algorithm is optimal in the sense of bilinear complexity, i.e.,
the minimum number of multiplications in $\Bbbk$.

###### Proof.

We first assume that $f(x)=x^{2}+\tau$. In this case, the multiplication
$(x_{1}+x_{2}\xi)(y_{1}+y_{2}\xi)$ can be computed by three multiplications:
$M_{1}=(x_{1}-x_{2})(y_{1}+\tau y_{2}),M_{2}=x_{1}y_{2},M_{3}=x_{2}y_{1}$,
since

$x_{1}y_{1}-\tau x_{2}y_{2}=M_{1}-\tau M_{2}+M_{3},\quad
x_{1}y_{2}+x_{2}y_{1}=M_{2}+M_{3}.$ (2)

Next we assume that $f(x)=x^{2}+x+\tau$. We consider: $M_{1}=x_{1}y_{1}$,
$M_{2}=x_{2}y_{2}$ and $M_{3}=(x_{1}-x_{2})(y_{1}-y_{2})$. It is
straightforward to verify that

$x_{1}y_{1}-cx_{2}y_{2}=M_{1}-cM_{2},\quad
x_{1}y_{2}+x_{2}y_{1}-x_{2}y_{2}=M_{1}-M_{3}.$ (3)

Now we assume that in both cases there exists an algorithm for the
multiplication in $\Bbbk[\xi]$, which costs two multiplications
$M_{1}^{\prime}$ and $M_{2}^{\prime}$. Then we have

$x_{1}y_{1}-\tau x_{2}y_{2},x_{1}y_{2}+x_{2}y_{1}-\delta
x_{2}y_{2}\in\operatorname{span}\\{M^{\prime}_{1},M^{\prime}_{2}\\},$

where $\delta=0$ if $f(x)=x^{2}+\tau$ and $\delta=1$ if $f(x)=x^{2}+x+\tau$.
Clearly $x_{1}y_{1}-\tau x_{2}y_{2}$ and $x_{1}y_{2}+x_{2}y_{1}-\delta
x_{2}y_{2}$ are not collinear, thus

$M^{\prime}_{1},M^{\prime}_{2}\in\operatorname{span}\\{x_{1}y_{1}-\tau
x_{2}y_{2},x_{1}y_{2}+x_{2}y_{1}-\delta x_{2}y_{2}\\}.$

Therefore, we may find constants $a,b,d,e\in\Bbbk$ such that

$\displaystyle M^{\prime}_{1}$ $\displaystyle=a(x_{1}y_{1}-\tau
x_{2}y_{2})+b(x_{1}y_{2}+x_{2}y_{1}-\delta
x_{2}y_{2})=ax_{1}y_{1}+bx_{1}y_{2}+bx_{2}y_{1}+(-\tau a-\delta b)x_{2}y_{2},$
$\displaystyle M^{\prime}_{2}$ $\displaystyle=c(x_{1}y_{1}-\tau
x_{2}y_{2})+d(x_{1}y_{2}+x_{2}y_{1}-\delta
x_{2}y_{2})=cx_{1}y_{1}+dx_{1}y_{2}+dx_{2}y_{1}+(-\tau c-\delta d)x_{2}y_{2}.$

Clearly, $ad-bc\neq 0$ and in particular, at least one of $a,b,c,d$ is
nonzero. We observe that $M_{1}^{\prime}$ is a multiplication, thus it can be
written as $(rx_{1}+sx_{2})(uy_{1}+vy_{2})$. This implies that

$\displaystyle a(-\tau a-\delta b)$ $\displaystyle=b^{2},$ (4) $\displaystyle
c(-\tau c-\delta d)$ $\displaystyle=d^{2}.$ (5)

If $f(x)=x^{2}+\tau$, then (4) reduces to $\tau a^{2}+b^{2}=0$. This implies
that $a=b=0$ since $-\tau$ is not a complete square in $\Bbbk$. Similarly, we
also have $c=d=0$. This contradicts the assumption that at least one of
$a,b,c,d$ is nonzero.

If $f(x)=x^{2}+x+\tau$ then (4) is $\tau a^{2}+ab+b^{2}=0$. We remark that
$a\neq 0$. Otherwise, we may derive $b=0$ which contradicts the fact that $ad-
bc\neq 0$. By the substitution $b^{\prime}=b/a$, we have
${b^{\prime}}^{2}+b^{\prime}+\tau=0$. This contradicts to the assumption that
the equation $x^{2}+x+\tau=0$ has no solution in $\Bbbk$. ∎

It is worthy to remark that we may regard the multiplication in $\Bbbk[\xi]$
as the bilinear map

$m:\Bbbk[\xi]\times\Bbbk[\xi]\to\Bbbk[\xi],\quad(x_{1}+x_{2}\xi,y_{1}+y_{2}\xi)\mapsto(x_{1}+x_{2}\xi)(y_{1}+y_{2}\xi).$

A multiplication in $\Bbbk$ means the multiplication between two
indeterminates from different inputs. For example, we count $x_{1}y_{1}$ as
one multiplication while we do not count $\tau x_{1}$ as a multiplication
since $\tau$ is a constant. Moreover, we do not allow multiplications between
indeterminates from the same input, such as $x_{1}x_{2}$. In terms of tensor
rank, Proposition 2.2 means the rank of the structure tensor of $m$ is exactly
three. According to [19, 37], it is known that every tensor in
$\mathbb{C}^{2}\otimes\mathbb{C}^{2}\otimes\mathbb{C}^{2}$ has rank at most
three. As a comparison, Proposition 2.2 deals with a particular tensor over an
arbitrary field $\Bbbk$. More importantly, it determines the exact value of
the rank of that tensor, together with a rank decomposition. An important
special case of Proposition 2.2 is $\Bbbk=\mathbb{R}$ and $f(x)=x^{2}+1$. In
this case, $\xi=i$ and $\Bbbk[\xi]=\mathbb{C}$. The algorithm presented in (2)
is the celebrated Gauss algorithm [34] for the multiplication of complex
numbers, whose optimality is proved in [45, 60].

##  3\. Gauss matrix multiplication

Let $n$ be a positive integer and let $\mathbb{L}$ be a field. We denote by
$M_{n}(\mathbb{L})$ the $\mathbb{L}$-algebra consisting of $n\times n$
matrices over $\mathbb{L}$. We denote by $\operatorname{GL}_{n}(\mathbb{L})$
the group of invertible $n\times n$ matrices over $\mathbb{L}$ and we consider
the multiplication map

$m_{n,\mathbb{L}}:\operatorname{M}_{n}(\mathbb{L})\times\operatorname{M}_{n}(\mathbb{L})\to\operatorname{M}_{n}(\mathbb{L}),\quad
m_{n,\mathbb{L}}(Z,W)=ZW$

and the inversion map

$\operatorname{inv}_{n,\mathbb{L}}:\operatorname{GL}_{n}(\mathbb{L})\to\operatorname{GL}_{n}(\mathbb{L}),\quad\operatorname{inv}_{n,\mathbb{L}}(Z)=Z^{-1}.$

Assume that $\mathbb{F}\simeq\Bbbk[x]/(f(x))\simeq\Bbbk[\xi]$ is a quadratic
field extension of $\Bbbk$, where $\xi$ is a root of $f\in\Bbbk[x]$ in
$\overline{\Bbbk}$. Here $f(x)$ is either $x^{2}+\tau$ or $x^{2}+x+\tau$. We
notice that $M_{n}(\mathbb{F})=M_{n}(\Bbbk)\otimes_{\Bbbk}\mathbb{F}$. Thus an
element in $Z\in M_{n}(\mathbb{F})$ can be written as $Z=A+\xi B$ where
$A,B\in M_{n}(\Bbbk)$. As a consequence of Proposition 2.2, we have the
following.

######  Proposition 3.1 (Gauss matrix multiplication).

Let $\Bbbk,\mathbb{F},n,f,\tau,\xi$ be as above. Suppose that $Z=A+\xi B$ and
$W=C+\xi D$ are two elements in $M_{n}(\mathbb{F})$ where $A,B,C,D\in
M_{n}(\Bbbk)$. If $f(x)=x^{2}+\tau$, then one can compute the product $ZW$ by

$\displaystyle M_{1}$ $\displaystyle=(A-B)(C+\tau D),\quad
M_{2}=AD,\quad,M_{3}=BC,$ $\displaystyle N_{1}$ $\displaystyle=M_{1}-\tau
M_{2}+M_{3},\quad N_{2}=M_{2}+M_{3},$ (6) $\displaystyle ZW$
$\displaystyle=N_{1}+\xi N_{2}.$

If $f(x)=x^{2}+x+\tau$, then one can compute the product $ZW$ by

$\displaystyle M_{1}$ $\displaystyle=AC,\quad
M_{2}=BD,\quad,M_{3}=(A-B)(C-D),$ $\displaystyle N_{1}$
$\displaystyle=M_{1}-\tau M_{2},\quad N_{2}=M_{1}-M_{3},$ (7) $\displaystyle
ZW$ $\displaystyle=N_{1}+\xi N_{2}.$

Moreover, (3.1) and (3.1) are optimal in the sense of minimum number of
multiplications in $M_{n}(\Bbbk)$.

We remind the readers that if we denote by $\mathcal{E}$ the collection of
algorithms for $m_{n,\mathbb{F}}$ using multiplications and additions of
indeterminates, together with scalar multiplications in $M_{n}(\Bbbk)$, then
the algorithm presented in Proposition 3.1 is an minimal element with respect
to the partial ordering on $\mathcal{E}$ induced by the number of
multiplications of indeterminates in $M_{n}(\Bbbk)$.

###### Proof.

It is straightforward to verify that algorithms in (3.1) and (3.1) indeed
compute $ZW$. To see the optimality, one may repeat the argument in the proof
of Proposition 2.2. ∎

##  4\. Frobenius matrix inversion

In this section, we discuss the matrix inversion on $M_{n}(\mathbb{F})$. To do
that, we observe that for $Z=A+\xi B,W=C+\xi
D\in\operatorname{GL}_{n}(\mathbb{F})$ where $A,B,C,D\in M_{n}(\Bbbk)$, $W$ is
the inverse of $Z$ if and only if

$(A+\xi B)(C+\xi D)=ZW=I_{n}.$ (8)

Thus we may have two equations for $C$ and $D$. According to Lemma 2.1, the
multiplication formula in $\mathbb{F}$ thus in $M_{n}(\mathbb{F})$ depends on
the form of $f$. In the following, we separate our discussions with respect to
the two normal forms of $f$.

###  4.1. First case: $f(x)=x^{2}+\tau$

We define

$\displaystyle\mathcal{S}_{1}$ $\displaystyle\coloneqq\\{Z=A+\xi
B\in\operatorname{GL}_{n}(\mathbb{F}):A,A+\tau
BA^{-1}B\in\operatorname{GL}_{n}(\Bbbk)\\},$ $\displaystyle\mathcal{S}_{2}$
$\displaystyle\coloneqq\\{Z=A+\xi B\in\operatorname{GL}_{n}(\mathbb{F}):B,\tau
B+AB^{-1}A\in\operatorname{GL}_{n}(\Bbbk)\\}.$

It is worth mentioning that the complement of
$\mathcal{S}_{1}\cup\mathcal{S}_{2}$ consists of $n\times n$
$\mathbb{F}$-matrices $Z=A+\xi B$ where $\det(A(A+\tau BA^{-1}B))=\det(B(\tau
B+AB^{-1}A))=0$. Therefore, $\mathcal{S}_{1}\cup\mathcal{S}_{2}$ is an open
dense (in Zariski topology) subset of $\mathbb{F}^{n\times
n}\simeq\Bbbk^{2\times n\times n}$. Thus with respect to any reasonable
probability measure on $\mathbb{\Bbbk}^{2\times n\times n}$, a random complex
matrix $Z$ lies in $\mathcal{S}_{1}\cup\mathcal{S}_{2}$ with probability one.
The restriction of $\operatorname{inv}_{n,\mathbb{F}}$ on $\mathcal{S}_{1}$
and $\mathcal{S}_{2}$ can be easily expressed in terms of
$\operatorname{inv}_{n,\Bbbk}$ and $m_{n,\Bbbk}$. Indeed, according to (8), we
have

$AC-\tau BD=I_{n},\quad AD+BC=0.$

Solving the equation for $C$ and $D$, we obtain the lemma that follows.

######  Lemma 4.1 (Frobenius inversion 1).

Let $\mathbb{F}=\Bbbk[x]/(f(x))$ where $f(x)=x^{2}+\tau$ is an irreducible
polynomial over $\Bbbk$. For each $Z=A+\xi
B\in\operatorname{GL}_{n}(\mathbb{F})$, where
$A,B\in\operatorname{M}_{n}(\Bbbk)$, we have

  1. (a)

If $Z\in\mathcal{S}_{1}$, then

$Z^{-1}=(A+\tau BA^{-1}B)^{-1}-\xi A^{-1}B(A+\tau BA^{-1}B)^{-1}.$ (9)

  2. (b)

If $Z\in\mathcal{S}_{2}$, then

$Z^{-1}=B^{-1}A(AB^{-1}A+\tau B)^{-1}-\xi(AB^{-1}A+\tau B)^{-1}.$ (10)

In particular, $\operatorname{inv}_{n,\mathbb{F}}$ on
$\mathcal{S}_{1}\cup\mathcal{S}_{2}$ can be evaluated by performing
$\operatorname{inv}_{n,\Bbbk}$ twice and $m_{n,\Bbbk}$ thrice.

We remark that if $\Bbbk=\mathbb{R},\mathbb{F}=\mathbb{C}$, then (9) and (10)
are the well-known inversion formulae for complex matrices first discovered by
Frobenius [7, 13, 15, 36, 40, 52, 55, 59, 64]. It is easy to turn (9) and (10)
into Algorithm 1.

Algorithm 1 Frobenius inversion 1

1:$Z=A+\xi B\in\mathcal{S}_{1}\cup\mathcal{S}_{2}$

2:inverse of $Z$

3:if $Z\in\mathcal{S}_{1}$ then

4: set $X=A,Y=B,\tau_{1}=1,\tau_{2}=\tau$;

5:else if $Z\in\mathcal{S}_{2}$ then

6: set $X=B,Y=A,\tau_{1}=\tau,\tau_{2}=1$;

7:end if

8:compute $X^{-1}$;

9:compute $X^{-1}Y$;

10:compute $YX^{-1}Y$;

11:compute $\tau_{1}X+\tau_{2}YX^{-1}Y$;

12:compute $J=(\tau_{1}X+\tau_{2}YX^{-1}Y)^{-1}$;

13:compute $K=X^{-1}Y(\tau_{1}X+\tau_{2}YX^{-1}Y)^{-1}$;

14:if $Z\in\mathcal{S}_{1}$ then return $Z^{-1}=J-\xi K$;

15:else if $Z\in\mathcal{S}_{2}$ then return $Z^{-1}=K-\xi J$;

16:end if

Note that Algorithm 1 reduces to the algorithm in [8, p. 15] and [38] when
$n=1$, $\Bbbk=\mathbb{R}$ and $\mathbb{F}=\mathbb{C}$.

One may expect that (9) and (10) are consequences of any of the
Sherman–Morrison–Woodbury (SMW) identities [20, 21, 23, 26, 50, 61, 62]:

$\displaystyle\left(A+B\right)^{-1}$
$\displaystyle=A^{-1}-A^{-1}(B^{-1}+A^{-1})^{-1}A^{-1}$
$\displaystyle=A^{-1}-A^{-1}\left(AB^{-1}+I_{n}\right)^{-1}$
$\displaystyle=A^{-1}-\left(A+AB^{-1}A\right)^{-1}$
$\displaystyle=A^{-1}-A^{-1}B\left(A+B\right)^{-1}$

However, this is not the case since SMW identities all involve the inversion
of a matrix over $\mathbb{F}$ while (9) and (10) only require inversions of
matrices over $\Bbbk$.

From the algebraic perspective, $\operatorname{inv}_{n,\mathbb{F}}$ is an
operation on $M_{n}(\mathbb{F})$, which is an $M_{n}(\Bbbk)$-bimodule. Hence
it is natural to discuss the computational complexity of
$\operatorname{inv}_{n,\mathbb{F}}$ over $M_{n}(\Bbbk)$. We observe that
Algorithm 1 performs $m_{n,\Bbbk}$ three times, $\operatorname{inv}_{n,\Bbbk}$
twice and real matrix addition once. In the following we prove the optimality
of Algorithm 1.

######  Theorem 4.2 (optimality 1).

Let $n\geq 2$ be a positive integer. Algorithm 1 is optimal in the sense of
least number of multiplications, inversions and additions in $M_{n}(\Bbbk)$.

###### Proof.

It is clear that one matrix addition over $\Bbbk$ is necessary to compute
$\operatorname{inv}_{n,\mathbb{F}}$. Let $\Gamma$ be a straight-line program
computing $Z^{-1}$ for $Z\in\mathcal{S}_{1}\cup\mathcal{S}_{2}$ over the
algebra $\operatorname{M}_{n}(\Bbbk)$. We denote by $L(\Gamma)$ the total
number of inversions and multiplications in $\operatorname{M}_{n}(\Bbbk)$
performed by $\Gamma$. If we denote by $\Gamma_{0}$ the straight-line program
for Algorithm 1, then it is clear that $L(\Gamma_{0})=5$. We assume that there
exists some straight-line program $\Gamma_{1}$ with $L(\Gamma_{1})=4$
computing $Z^{-1}$.

We first prove that $\Gamma_{1}$ performs $\operatorname{inv}_{n,\Bbbk}$ at
least twice. Otherwise, $\Gamma_{1}$ can compute $Z^{-1}$ by only one
$\Bbbk$-matrix inversion. Without loss of generality, we may assume that
$Z\in\mathcal{S}_{1}$ so that $(A+\tau BA^{-1}B)^{-1}$ can be computed by at
most one $\Bbbk$-matrix inversion. Furthermore, we may also assume that $B$ is
invertible. Now we have

$A(A+\tau BA^{-1}B)^{-1}=(I_{n}+\tau(BA^{-1})^{2})^{-1}.$

We consider the case where $B=\Lambda A$ where
$\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{n})$ is a diagonal
matrix. Then $(I_{n}+\tau(BA^{-1})^{2})^{-1}$ is the diagonal matrix whose
diagonal elements are $1/(1+\tau\lambda_{j}^{2}),j=1,\dots,n$. This indicates
that $\Gamma_{1}$ requires at least two inversions in $M_{n}(\Bbbk)$ to
compute $(I_{n}+\tau(BA^{-1})^{2})^{-1}$, since $A,B$ are not necessarily
commuting.

Next, we prove that $\Gamma_{1}$ requires at least three $\Bbbk$-matrix
multiplications. Again, we assume that $Z=A+\xi B\in\mathcal{S}_{1}$. Then the
output of $\Gamma_{1}$ is the pair

$\left((A+\tau BA^{-1}B)^{-1},\quad A^{-1}B(A+\tau BA^{-1}B)^{-1}\right).$

We proceed by contradiction. Suppose that $\Gamma_{1}$ costs at most two
$\Bbbk$-matrix multiplications. Since $\Gamma_{1}$ is a straight-line program,
we may assume without loss of generality, that $\Gamma_{1}$ first computes
$(A+\tau BA^{-1}B)^{-1}$ and then it computes $A^{-1}B(A+\tau BA^{-1}B)^{-1}$.
We claim that $\Gamma_{1}$ already requires two $\Bbbk$-matrix multiplications
to compute $(A+\tau BA^{-1}B)^{-1}$. To that end, we notice that

$BA^{-1}B=\tau^{-1}\left(\left((A+\tau BA^{-1}B)^{-1}\right)^{-1}-A\right).$

This implies that one can compute $BA^{-1}B$ using the same number of
$\Bbbk$-matrix multiplications as $(A+\tau BA^{-1}B)^{-1}$. However, it is
clear that computing $BA^{-1}B$ by one multiplication is impossible and this
proves our claim. Thus $\Gamma_{1}$ computes $A^{-1}B(A+\tau BA^{-1}B)^{-1}$
as a linear combination of $(A+\tau BA^{-1}B)^{-1}$ and intermediate outputs.
However, this is not possible. Indeed, we may again consider the case where
$B=\Lambda A$ and $\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{n})$
is a diagonal matrix. Therefore, we have

$\left((A+\tau BA^{-1}B)^{-1},\quad A^{-1}B(A+\tau
BA^{-1}B)^{-1}\right)=(A^{-1}(I_{n}+\tau\Lambda^{2})^{-1},A^{-1}(I_{n}+\tau\Lambda^{2})^{-1}\Lambda).$

Moreover, $A^{-1}(I_{n}+\tau\Lambda^{2})^{-1}(I_{n}-\Lambda)$ can be computed
as a linear combination of intermediate outputs. In particular, one can
compute both $A^{-1}(I_{n}+\tau\Lambda^{2})^{-1}$ and
$A^{-1}(I_{n}+\tau\Lambda^{2})^{-1}\Lambda$ by two multiplications in
$M_{n}(\Bbbk)$, which is ridiculous. ∎

###  4.2. Second case: $f(x)=x^{2}+x+\tau$

We define

$\displaystyle\mathcal{T}_{1}$ $\displaystyle\coloneqq\left\\{Z=A+\xi B\in
M_{n}(\mathbb{F}):A,B\in\Bbbk,A-B,A+\tau
B(A-B)^{-1}B\in\operatorname{GL}_{n}(\Bbbk)\right\\},$
$\displaystyle\mathcal{T}_{2}$ $\displaystyle\coloneqq\left\\{Z=A+\xi B\in
M_{n}(\mathbb{F}):A,B\in\Bbbk,B,AB^{-1}A-A+\tau
B\in\operatorname{GL}_{n}(\Bbbk)\right\\}.$

Again, $\mathcal{S}_{1}\cup\mathcal{S}_{2}$ is an open dense subset of
$\mathbb{F}^{n\times n}\simeq\Bbbk^{2\times n\times n}$ in Zariski topology.
According to (8), we have

$AC-\tau BD=I_{n},\quad AD+BC-BD=0.$

Solving the equation for $C$ and $D$, we have the following

######  Lemma 4.3 (Frobenius inversion 2).

Let $\Bbbk$ be a field of characteristic $2$ and let
$\mathbb{F}=\Bbbk[x]/(f(x))$ where $f(x)=x^{2}+x+\tau$ is an irreducible
polynomial over $\Bbbk$. For each $Z=A+\xi
B\in\operatorname{GL}_{n}(\mathbb{F})$ where $A,B\in M_{n}(\Bbbk)$, we have

  1. (a)

If $Z\in\mathcal{T}_{1}$, then

$Z^{-1}=(A+\tau B(A-B)^{-1}B)^{-1}-\xi(A-B)^{-1}B(A+\tau B(A-B)^{-1}B)^{-1}.$
(11)

  2. (b)

If $Z\in\mathcal{T}_{2}$, then

$Z^{-1}=(B^{-1}A-I_{n})(AB^{-1}A-A+\tau B)^{-1}-\xi(AB^{-1}A-A+\tau B)^{-1}.$
(12)

As a direct consequence of (11) and (12), we obtain Algorithm 2.

Algorithm 2 Frobenius inversion 2

1:$Z=A+\xi B\in\mathcal{T}_{1}\cup\mathcal{T}_{2}$

2:inverse of $Z$

3:if $Z\in\mathcal{T}_{1}$ then

4: compute $X_{1}=(A-B)^{-1}$;

5: compute $X_{2}=X_{1}B$;

6: compute $X_{3}=A+\tau BX_{2}$;

7: compute $X_{4}=X_{3}^{-1}$;

8: compute $X_{5}=X_{2}X_{4}$; return $Z^{-1}=X_{4}-\xi X_{5}$

9:else if $Z\in\mathcal{T}_{2}$ then

10: compute $X_{1}=B^{-1}$;

11: compute $X_{2}=X_{1}A$;

12: compute $X_{3}=AX_{2}-A+\tau B$;

13: compute $X_{4}=X_{3}^{-1}$;

14: compute $X_{5}=(X_{2}-I_{n})X_{4}$; return $Z^{-1}=X_{5}-\xi X_{4}$

15:end if

It is straightforward to verify that Algorithm 2 costs two matrix inversions
and three matrix multiplications over $\Bbbk$. In fact, by a similar argument
for Theorem 4.2, we can prove that Algorithm 2 is optimal.

######  Theorem 4.4 (optimality 2).

Let $n\geq 2$ be a positive integer. Algorithm 2 is optimal in the sense of
least number of multiplications, inversions and additions in $M_{n}(\Bbbk)$.

###  4.3. an application

To conclude this section, we discuss an application of Algorithms 1 and 2. We
consider a tower of quadratic field extensions:

$\Bbbk=\mathbb{K}_{0}\subsetneq\mathbb{K}_{1}\subsetneq\cdots\subsetneq\mathbb{K}_{m},$
(13)

where $[\mathbb{K}_{j}:\mathbb{K}_{j-1}]=2,j=1,\dots,m$. For each
$j=1,\dots,m$ there exists some $\xi_{j}\in\mathbb{K}_{j}$ such that
$\mathbb{K}_{j}=\mathbb{K}_{j-1}[\xi_{j}]$. We denote by $f_{j}\in\Bbbk[x]$
the minimal polynomial of $\xi_{j}$. Then $f_{j}$ is a monic irreducible
quadratic polynomial. We may also assume that $f_{j}$ is in the normal form:

$f_{j}(x)=x^{2}+\tau_{j}~{}\text{or}~{}f_{j}(x)=x^{2}+x+\tau_{j},\quad
j=1,\dots,m.$

It is obvious that $[\mathbb{K}_{m}:\Bbbk]=2^{m}$ and an element $x$ in
$\mathbb{K}_{m}$ can be written as

$x=\sum_{\alpha\in\\{0,1\\}^{m}}c_{\alpha}\xi^{\alpha},$

where for each $\alpha=(\alpha_{1},\dots,\alpha_{m})\in\\{0,1\\}^{m}$,
$\xi^{\alpha}=\xi_{1}^{\alpha_{1}}\cdots\xi_{m}^{\alpha_{m}}$ and
$c_{\alpha}\in\Bbbk$. Moreover, we may regard $\mathbb{K}_{m}$ as a quotient
ring of the polynomial ring in $m$ variables, i.e.,

$\mathbb{K}_{m}\simeq\Bbbk[x_{1},\dots,x_{m}]/(f_{1},\dots,f_{m})=\bigotimes_{j=1}^{m}\left(\Bbbk[x]/(f_{j})\right).$
(14)

A particularly interesting example of $\mathbb{K}_{m}$ is $\mathbb{Q}[S]$
where $S=\\{\sqrt{a_{j}}:a_{j}\in\mathbb{Q},1\leq j\leq m\\}$. It is proved in
[6] that

$\mathbb{Q}\subsetneq\mathbb{Q}[\sqrt{a_{1}}]\subsetneq\mathbb{Q}[\sqrt{a_{1}},\sqrt{a_{2}}]\subsetneq\cdots\subsetneq\mathbb{Q}[S]$

is a tower of quadratic field extensions if the product of any nonempty subset
of $S$ is not in $\mathbb{Q}$. In this case, we have
$\mathbb{K}_{j}=\mathbb{Q}[\sqrt{a}_{1},\dots,\sqrt{a}_{j}]$ and
$f_{j}(x)=x^{2}-a_{j},j=1,\dots,m$. Another commonly seen example of the tower
(13) is

$\mathbb{Q}\subsetneq\mathbb{Q}[a^{\frac{1}{2}}]\subsetneq\cdots\subsetneq\mathbb{Q}[a^{\frac{1}{2^{m}}}],$

where $a\in\mathbb{Q}$ is not a complete square.

We also remark that (13) is analogous to the tower of multicomplex number
systems [48]:

$\mathbb{R}\subsetneq\mathbb{C}_{1}\subsetneq\cdots\subsetneq\mathbb{C}_{m},$

where $\mathbb{C}_{j},1\leq j\leq m$ is defined inductively as

$\displaystyle\mathbb{C}_{1}$ $\displaystyle=\mathbb{C},$
$\displaystyle\mathbb{C}_{k+1}$
$\displaystyle=\\{a+bi_{k+1}:a,b\in\mathbb{C}_{k}\\},\quad k\geq 1$

where $i_{k+1}^{2}=-1$ and $i_{k}i_{l}=i_{l}i_{k}$ for any $k,l\geq 1$.
However, since $\mathbb{C}_{1}=\mathbb{C}$ is algebraic closed,
$\mathbb{C}_{j}$ is not even a field extension of $\mathbb{C}_{1}$ for $j\geq
2$. We refer interested readers to [53, 48] for detailed discussions of
analytic properties of $\mathbb{C}_{m}$. Another example i

We observe that
$M_{n}(\mathbb{K}_{m})=M_{n}(\Bbbk)\otimes_{\Bbbk}\mathbb{K}_{m}$, which
implies that a matrix $Z\in M_{n}(\mathbb{K}_{m})$ can be written as

$Z=\sum_{\alpha\in\\{0,1\\}^{m}}C_{\alpha}\xi^{\alpha}.$ (15)

On the other hand, we also have

$M_{n}(\Bbbk)\otimes_{\Bbbk}\mathbb{K}_{m}=M_{n}(\Bbbk)\otimes_{\Bbbk}\mathbb{K}_{1}\otimes_{\mathbb{K}_{1}}\cdots\otimes_{\mathbb{K}_{m-1}}\mathbb{K}_{m}.$

Thus $Z\in\mathbb{K}_{m}$ can also be inductively constructed as follows:

$\displaystyle Z$ $\displaystyle=A_{0}+\xi_{m}A_{1},$ $\displaystyle
A_{\beta}$
$\displaystyle=A_{0,\beta}+\xi_{m-j}A_{1,\beta},\beta\in\\{0,1\\}^{j},j=1,\dots,m-1,$
(16)

where $A_{\beta}\in\operatorname{M}_{n}(\mathbb{K}_{m-|\beta|})$. We record in
the next lemma the relation between the two expressions of $Z$ in (15) and
(16).

######  Lemma 4.5.

Let $C_{\alpha}$ and $A_{\beta}$ be as in (15) and (16) respectively, where
$\alpha\in\\{0,1\\}^{m}$ and $\beta\in\bigcup_{j=1}^{m}\\{0,1\\}^{j}$. For
each $1\leq j\leq m$,

$Z=\sum_{\beta\in\\{0,1\\}^{j}}A_{\beta}\xi_{m-j+1}^{\beta_{1}}\cdots\xi_{m}^{\beta_{j}}$

Moreover, for each $\beta\in\\{0,1\\}^{j}$, we have

$A_{\beta}=\sum_{\beta^{\prime}\in\\{0,1\\}^{m-j}}C_{\beta^{\prime},\beta}\xi_{1}^{\beta^{\prime}_{1}}\cdots\xi_{m-j}^{\beta^{\prime}_{m-j}},\quad\beta^{\prime}\in\\{0,1\\}^{m-j}.$

In particular, $C_{\alpha}=A_{\alpha}$.

###### Proof.

We proceed by induction on $j$. Clearly the formula holds for $j=1$ by (16).
Assume that the formula holds for $j=s$, i.e.,

$Z=\sum_{\beta\in\\{0,1\\}^{s}}A_{\beta}\xi_{m-s+1}^{\beta_{1}}\cdots\xi_{m}^{\beta_{s}}.$

To prove the formula for $j=s+1$, we notice that
$A_{\beta}=A_{0,\beta}+\xi_{m-s}A_{1,\beta}$, thus

$Z=\sum_{\beta\in\\{0,1\\}^{s}}(A_{0,\beta}+\xi_{m-s}A_{1,\beta})\xi_{m-s+1}^{\beta_{1}}\cdots\xi_{m}^{\beta_{s}}=\sum_{\gamma\in\\{0,1\\}^{s+1}}A_{\gamma}\xi_{m-s}^{\gamma_{1}}\cdots\xi_{m}^{\gamma_{s+1}}$

and this completes the induction. The moreover part follows easily by
comparing (16) with (15). ∎

Lemma 4.5 provides us a method to compute the multiplication and inversion in
$M_{n}(\mathbb{K}_{m})$.

######  Proposition 4.6 (multiplication in $M_{n}(\mathbb{K}_{m})$).

There exists an algorithm for the multiplication in $M_{n}(\mathbb{K}_{m})$,
which costs $3^{m}$ multiplications in $M_{n}(\Bbbk)$.

###### Proof.

Let $Z,W\in M_{n}(\mathbb{K}_{m})$. According to (16), we may write

$Z=A_{0}+\xi_{m}A_{1},\quad W=B_{0}+\xi_{m}B_{1}.$

Thus one can compute $ZW$ in terms of $A_{0},A_{1},B_{0},B_{1}$ by three
multiplications in $M_{n}(\mathbb{K}_{m-1})$ by Proposition 3.1. Each
multiplication in $M_{n}(\mathbb{K}_{m-1})$ costs three multiplications in
$M_{n}(\mathbb{K}_{m-2})$ again by Proposition 3.1. Repeat the above process
until we arrive at multiplications in $M_{n}(\mathbb{K}_{0})=M_{n}(\Bbbk)$.
Thus the total cost of multiplications in $M_{n}(\Bbbk)$ is $3^{m}$. ∎

According to (15), we can also write $Z$ and $W$ as
$Z=\sum_{\alpha\in\\{0,1\\}^{m}}C_{\alpha}\xi^{\alpha}$ and
$W=\sum_{\beta\in\\{0,1\\}^{m}}D_{\beta}\xi^{\alpha}$ respectively. Here
$C_{\alpha},D_{\beta}\in M_{n}(\Bbbk)$ for $\alpha,\beta\in\\{0,1\\}^{m}$.
Thus one can compute $ZW$ via the formula

$ZW=\sum_{\gamma\in\\{0,2\\}^{m}}\left(\sum_{\begin{subarray}{c}\alpha+\beta=\gamma\\\
\alpha,\beta\in\\{0,1\\}^{m}\end{subarray}}C_{\alpha}D_{\beta}\right)\xi^{\gamma}.$

However, it is obvious that the above formula costs $4^{m}$ multiplications in
$M_{n}(\Bbbk)$. Hence the algorithm presented in Proposition 4.6 reduces the
complexity of evaluating $m_{n,\mathbb{K}_{m}}$ over $M_{n}(\Bbbk)$ from
$O(N^{2})$ to $O(N^{\log_{2}3})$, where $N=2^{m}$.

Due to Proposition 3.1, one can recognize the algorithm in Proposition 4.6 as
an analogue of the Karatsuba algorithm [30] for fast integer multiplication.
By (14), we may also regard the algorithm in Proposition 4.6 as an analogue of
the multidimensional fast Fourier transform [54].

We also remark that the algorithm presented in Proposition 4.6 relies on the
technique called divide and conquer, which is employed by Strassen to design
the first algorithm [57] for $n\times n$ matrix multiplication whose
complexity is smaller than $O(n^{3})$. In our case, however, the size $n$ of
matrices is fixed, while the level $m$ of the tower (13) varies.

######  Proposition 4.7 (inversion in $M_{n}(\mathbb{K}_{m})$).

There exists an algorithm for the inversion of a generic element in
$M_{n}(\mathbb{K}_{m})$, which costs $3(3^{m}-2^{m})$ multiplications and
$2^{m}$ inversions in $M_{n}(\Bbbk)$.

###### Proof.

For $1\leq j\leq m$ and a generic $Z\in M_{n}(\mathbb{K}_{j})$. We write
$Z=A_{0}+\xi_{m}A_{1}$ where $A_{0},A_{1}\in M_{n}(\mathbb{K}_{j-1})$. Since
$Z$ is generic, Algorithm 1 (resp. Algorithm 2) is applicable for
$f_{j}(x)=x^{2}+\tau_{j}$ (resp. $f_{j}(x)=x^{2}+x+\tau_{j}$). This costs
three multiplications and two inversions in $M_{n}(\mathbb{K}_{j-1})$. Thus we
have

$\operatorname{inv}_{n,\mathbb{K}_{j}}=3m_{n,\mathbb{K}_{j-1}}+2\operatorname{inv}_{n,\mathbb{K}_{j-1}},1\leq
j\leq m.$

Moreover, Proposition 3.1 implies
$m_{n,\mathbb{K}_{j}}=3m_{n,\mathbb{K}_{j-1}}$. Inductively, we may derive
that

$\operatorname{inv}_{n,\mathbb{K}_{m}}=3(3^{m}-2^{m})m_{n,\Bbbk}+2^{m}\operatorname{inv}_{n,\Bbbk}.$

∎

##  5\. General matrix inversion

In Section 4, we investigate properties of the Frobenius inversion from the
symbolic perspective. In practice, the most important quadratic field
extension is $\mathbb{R}\subseteq\mathbb{C}$. The rest of this paper is
devoted to the discussion of numerical properties of the Frobenius inversion.
To this end, we consider the case where $\Bbbk=\mathbb{R}$ and
$\mathbb{F}=\mathbb{C}$, which is probably the most important quadratic field
extension.

In this section, we compare the computational complexity of Algorithm 1 with
the usual complex matrix inversion algorithm based on LU decomposition, which
is widely employed in various main stream platforms for numerical computing
such as MATLAB, Maple, Julia and Python.

###  5.1. Frobenius inversion v.s. invertion via LU decomposition

To begin with, we first recall the algorithm for matrix inversion via LU
decomposition [26], which we reproduce in Algorithm 3 for ease of reference.
We notice that the main idea behind Algorithm 3 is that $B=A^{-1}$ if and only
if $BA=I_{n}$. It is obviously true that we may replace $BA=I_{n}$ by
$AB=I_{n}$ in the above. Accordingly, we obtain Algorithm 4 which is slightly
different from Algorithm 3.

Algorithm 3 matrix inversion via LU decomposition

1:$A\in\operatorname{GL}_{n}(\Bbbk)$

2:inverse of $A$

3:compute LU factorization of $A=LU$;

4:compute $U^{-1}$;

5:solve for $X$ from $XL=U^{-1}$;

6:return $X$;

Algorithm 4 matrix inversion via LU decomposition

1:$A\in\operatorname{GL}_{n}(\Bbbk)$

2:inverse of $A$

3:compute LU factorization of $A=LU$;

4:compute $L^{-1}$;

5:solve for $X$ from $UX=L^{-1}$;

6:return $X$;

Before we proceed, we fix some notations. Let $\mathcal{A}$ be an algorithm
for real matrix multiplication. We denote by
$T_{\operatorname{mult}}^{\mathcal{A}}(n)$ the average running time of
$\mathcal{A}$ on pairs of $n\times n$ real matrices. In addition, we let
$T_{\operatorname{inv}}^{\mathcal{A}}(n)$ be the average running time of
Algorithm 3 on invertible $n\times n$ real matrices, in which real matrix
multiplications are computed by $\mathcal{A}$. By symmetry, the average
running time of Algorithm 4 on invertible $n\times n$ real matrices is also
$T_{\operatorname{inv}}^{\mathcal{A}}(n)$. Now with these notations, we are
ready to present our threshold theorem.

######  Theorem 5.1 (threshold).

Let $\mathcal{A}$ be an algorithm for real matrix multiplication. Assume that
the running time of $\mathcal{A}$ on pairs of $n\times n$ matrices of which at
least one is upper or lower triangular is $\lambda
T_{\operatorname{mult}}^{\mathcal{A}}(n)$ for some $0<\lambda\leq 1$. Then
Algorithm 1 is asymptotically faster than Algorithm 3 over $\mathbb{C}$ if and
only if
$\lim_{n\to\infty}\left(T_{\operatorname{inv}}^{\mathcal{A}}(n)/T_{\operatorname{mult}}^{\mathcal{A}}(n)\right)>1+\lambda/2$.
In particular, if $\mathcal{A}$ is the usual matrix multiplication algorithm,
then Algorithm 1 is asymptotically faster than Algorithm 3 over $\mathbb{C}$
if and only if
$\lim_{n\to\infty}\left(T_{\operatorname{inv}}^{\mathcal{A}}(n)/T_{\operatorname{mult}}^{\mathcal{A}}(n)\right)>5/4$.

###### Proof.

We first show that the running time of Algorithm 1 is dominated by
$2T_{\operatorname{inv}}^{\mathcal{A}}(n)+5T_{\operatorname{mult}}^{\mathcal{A}}(n)/2$.
In fact, we notice that the first two steps in Algorithm 1 ( i.e., computing
$X^{-1}$ and $X^{-1}Y$) can be combined into one step by solving for $F$ from
$XF=Y$. To that end, we compute an LU- decomposition $X=LU$. After that, we
compute $L^{-1}$ and $L^{-1}Y$. Finally, we solve for $F$ from $UF=L^{-1}Y$.
As a comparison, we note that the only step that is not included in Algorithm
4 is the computation of $L^{-1}Y$.

Since $L^{-1}$ is lower triangular, multiplying $L^{-1}$ with $Y$ takes
$\lambda T_{\operatorname{mult}}^{\mathcal{A}}(n)$ operations. Therefore, the
first two steps in Algorithm 1 takes
$T_{\operatorname{inv}}^{\mathcal{A}}(n)+\lambda
T_{\operatorname{mult}}^{\mathcal{A}}(n)$ time. Then, computing $YX^{-1}Y$
requires one matrix multiplication, computing $X+YX^{-1}Y$ requires one matrix
addition, computing $J=(X+YX^{-1}Y)^{-1}$ requires one matrix inversion, and
computing $K=X^{-1}YJ$ requires one matrix multiplication. Since matrix
addition takes $O(n^{2})$ flops, it is not the dominant term in the
computation time of Algorithm 1 and we can omit that in this analysis. To sum
up, the running time of Algorithm 1 is dominated by
$2T_{\operatorname{inv}}^{\mathcal{A}}(n)+(2+\lambda)T_{\operatorname{mult}}^{\mathcal{A}}(n)$.

Next, we consider the running time of Algorithm 3 for complex matrices. We
prove that the running time of Algorithm 3 is dominated by
$4T_{\operatorname{inv}}^{\mathcal{A}}(n)$. Note that the complex addition
takes $2$ real flops and the complex multiplication takes $6$ real flops. In
Algorithm 3, there are “roughly” the same number of additions and
multiplications. Algorithm 3 contains roughly three operations: computing LU
factorization, computing forward substitution, and computing the inverse of an
upper triangular matrix. Computing LU factorization using Gaussian Elimination
requires roughly the same number of additions and multiplications. Computing
forward substitution also requires roughly the same number of additions and
multiplications. Finally, according to Method $1$ in page $263$ of [26],
inverting a triangular matrix can be done by a sequence of forward
substitutions. Thus, inverting a triangular matrix also requires roughly the
same number of additions and multiplications. Therefore, the running time of
Algorithm 3 over $\mathbb{C}$ is dominated by
$4T_{\operatorname{inv}}^{\mathcal{A}}(n)$.

Now, Algorithm 1 is faster than Algorithm 3 if and only if for $n$
sufficiently large,

$4T_{\operatorname{inv}}^{\mathcal{A}}(n)>2T_{\operatorname{inv}}^{\mathcal{A}}(n)+(2+\lambda)T_{\operatorname{mult}}^{\mathcal{A}}(n),$

which is equivalent to

$\lim_{n\to\infty}\left(\frac{T_{\operatorname{inv}}^{\mathcal{A}}(n)}{T_{\operatorname{mult}}^{\mathcal{A}}(n)}\right)>1+\frac{\lambda}{2}.\qed$

It is remarkable that the inequality
$\lim_{n\to\infty}\left(T_{\operatorname{inv}}^{\mathcal{A}}(n)/T_{\operatorname{mult}}^{\mathcal{A}}(n)\right)>5/4$
holds in MATLAB. In Section 7, we shall see by numerical examples that
Algorithm 1 is indeed faster than Algorithm 3 in MATLAB, which confirms
Theorem 5.1. Moreover, according to Theorem 5.1, for any $\mathcal{A}$,
Algorithm 1 is asymptotically faster than Algorithm 3 if
$\lim_{n\to\infty}\left(T_{\operatorname{inv}}^{\mathcal{A}}(n)/T_{\operatorname{mult}}^{\mathcal{A}}(n)\right)>3/2$.

We conclude this subsection by a remark on solving systems of linear equations

$(A+iB)(x+iy)=c+id$ (17)

by the Frobenius inversion, where $A,B\in\mathbb{R}^{n\times n}$, and
$x,y,c,d\in\mathbb{R}^{n}$. Namely, we can first compute $(A+iB)^{-1}$ by
Algorithm 1 and then compute $(A+iB)^{-1}(c+id)$. For a single linear system
(17), it is obviously more efficient to solve it by the LU decomposition
together with backward and forward substitutions. However, as pointed out in
[11], it is common in scientific computing that one needs to solve (17)
repeatedly with the same $(A+iB)$ but different $(c+di)$. To be more precise,
we have

$(A+iB)(x^{(k)}+iy^{(k)})=c^{(k)}+id^{(k)},\quad k=1,\dots,K,$

where $K$ is much larger than $n$. In this scenario, inverting $(A+iB)^{-1}$
by Algorithm 1 is more favourable in the sense of computational efficiency.

###  5.2. rounding error analysis

In this subsection, we provide a rounding error analysis for Algorithm 1. To
do that, we denote by $\mathsf{u}$ the unit roundoff. Given an $n\times n$
matrix $X$, we denote by $L_{X}$ (resp. $U_{X}$) the computed lower (resp.
upper) triangular factor in the LU decomposition of $X$, i.e., $X=L_{X}U_{X}$.
We also denote by $|X|$ the matrix whose entries are absolute values of
elements of $X$.

######  Theorem 5.2.

Let $Z=A+iB\in\mathcal{S}_{1}\cup\mathcal{S}_{2}$ be the input of Algorithm 1
and $W$ be the output of Algorithm 1. Let $X=A$ and $Y=B$ if
$Z\in\mathcal{S}_{1}$ and let $X=B$ and $Y=A$ if $Z\in\mathcal{S}_{2}$. Then
we have

$\begin{split}|Z^{-1}-W|&\leq
O(n)\left(|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\cdot|J|\right.\\\
&\left.+(|X^{-1}Y|+I)|J|(|X|+|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|+|L_{P}|\cdot|U_{P}|)|J|\right)\mathsf{u},\end{split}$

where $\leq$ between matrices holds componentwise and $P$ is the computed
value of $X+YX^{-1}Y$ by Algorithm 1 and $J=(X+YX^{-1}Y)^{-1}$.

The proof of Theorem 5.2 relies on the following lemma.

######  Lemma 5.3.

Let $A\in\operatorname{GL}_{n}(\mathbb{R})$ and $B\in\mathbb{R}^{n\times n}$
be such that $A+\mathsf{u}B\in\operatorname{GL}_{n}(\mathbb{R})$, where
$\mathsf{u}$ is the unit roundoff. Then,

$(A+\mathsf{u}B)^{-1}=A^{-1}-\mathsf{u}A^{-1}BA^{-1}+O(\mathsf{u}^{2}).$

###### Proof.

Note that

$\begin{split}(A+\mathsf{u}B)^{-1}&=(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}\\\
&\stackrel{{\scriptstyle(a)}}{{=}}(I_{n}-\mathsf{u}A^{-1}B+\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B)A^{-1}\\\
&=A^{-1}-\mathsf{u}A^{-1}BA^{-1}+O(\mathsf{u}^{2}),\end{split}$

where $(a)$ follows from the fact that

$\begin{split}&\quad(I_{n}+\mathsf{u}A^{-1}B)(I_{n}-\mathsf{u}A^{-1}B+\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B)\\\
&=I_{n}-\mathsf{u}A^{-1}B+\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B+\mathsf{u}A^{-1}B-\mathsf{u}^{2}A^{-1}BA^{-1}B\\\
&\quad+\mathsf{u}^{3}A^{-1}BA^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B\\\
&=I_{n}+\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B-\mathsf{u}^{2}A^{-1}BA^{-1}B+\mathsf{u}^{3}A^{-1}BA^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B\\\
&=I_{n}+\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B-\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B\\\
&\quad+\mathsf{u}^{3}A^{-1}BA^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B\\\
&=I_{n}+\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B-\mathsf{u}^{2}A^{-1}B(I_{n}+\mathsf{u}A^{-1}B)^{-1}A^{-1}B\\\
&=I_{n}.\end{split}$

∎

The proof of Theorem 5.2 also requires the following facts in [26]:

  * •

Let $A\in\operatorname{GL}_{n}(\mathbb{R})$ and let $\widehat{X}$ be the
computed inverse of $A$ by Algorithm 3. Then,

$\displaystyle\widehat{X}$
$\displaystyle=A^{-1}+O(n)|\widehat{X}|\cdot|L|\cdot|U|\cdot|A^{-1}|\mathsf{u}+O(\mathsf{u}^{2})$
$\displaystyle=A^{-1}+O(n)|A^{-1}|\cdot|L|\cdot|U|\cdot|A^{-1}|\mathsf{u}+O(\mathsf{u}^{2}),$
(18)

where $A=LU$ is the computed LU decomposition of $A$.

  * •

Let $A,B\in\mathbb{R}^{n\times n}$ and let $\widehat{C}$ be the computed
product of $A$ and $B$. Then,

$\widehat{C}=AB+n|A|\cdot|B|\mathsf{u}+O(\mathsf{u}^{2}).$ (19)

  * •

Let $a,b\in\mathbb{R}$ and let $\hat{c}$ be the computed sum of $a$ and $b$.
Then,

$\hat{c}=(a+b)+(|a|+|b|)\mathsf{u}+O(\mathsf{u}^{2}).$ (20)

  * •

Let $a,b\in\mathbb{R}$ and let $\hat{c}$ be the computed value of $a-b$. Then,

$\hat{c}=(a-b)+(|a|+|b|)\mathsf{u}+O(\mathsf{u}^{2}).$ (21)

Now, we are able to prove Theorem 5.2. In the following, we denote by
$\widehat{X}$ the the computed value of a matrix $X$.

###### Proof of Theorem 5.2.

Let $L_{X}$ and $U_{X}$ be the computed LU factors of $X$. Then according to
(• ‣ 5.2), we have

$\widehat{X^{-1}}=X^{-1}+O(n)|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\mathsf{u}+O(\mathsf{u}^{2}).$

Let $H_{1}=X^{-1}Y$ and let $\widehat{H}_{1}$ be the computed value of
$H_{1}$. By (• ‣ 5.2) and (19), we derive

$\begin{split}\widehat{H}_{1}&=X^{-1}Y+O(n)|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\mathsf{u}+O(n)|X^{-1}|\cdot|Y|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=H_{1}+O(n)|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\mathsf{u}+O(\mathsf{u}^{2}),\end{split}$
(22)

since
$|X^{-1}|\cdot|Y|=|I_{n}|\cdot|X^{-1}|\cdot|Y|=O(|X^{-1}L_{X}U_{X}|\cdot|X^{-1}||Y|)=O(|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|)$.

Let $H_{2}=YX^{-1}Y$ and $\widehat{H}_{2}$ be the computed value of
$H_{2}$.Then (• ‣ 5.2) and (19) again imply

$\begin{split}\widehat{H}_{2}&=YX^{-1}Y+O(n)|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\mathsf{u}+O(n)|Y|\cdot|X^{-1}|\cdot|Y|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=H_{2}+O(n)|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\mathsf{u}+O(\mathsf{u}^{2}),\end{split}$

since
$|Y|\cdot|X^{-1}|\cdot|Y|=|Y|O(|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|)=O(|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|)$.

Let $H_{3}=X+YX^{-1}Y$ and $\widehat{H_{3}}$ be the computed value of $H_{3}$.
Then according to (• ‣ 5.2) and (19), we have

$\widehat{H_{3}}=H_{3}+O(n)(|X|+|Y||X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|)\mathsf{u}+O(\mathsf{u}^{2}).$

We recall that $P$ is the computed value of $X+YX^{-1}Y$ and
$J=(X+YX^{-1}Y)^{-1}$, thus Lemma 5.3 indicates that

$\displaystyle\widehat{P^{-1}}$
$\displaystyle=H_{3}^{-1}-O(n)H_{3}^{-1}(|X|+|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|)H_{3}^{-1}\mathsf{u}$
$\displaystyle+O(n)H_{3}^{-1}|L_{P}|\cdot|U_{P}|\cdot|H_{3}^{-1}|\mathsf{u}+O(\mathsf{u}^{2})$
(23)
$\displaystyle=J+O(n)|H_{3}^{-1}|(|X|+|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|+|L_{P}|\cdot|U_{P}|)|H_{3}^{-1}|\mathsf{u}+O(\mathsf{u}^{2})$
$\displaystyle=J+O(n)|J|(|X|+|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|+|L_{P}|\cdot|U_{P}|)|J|\mathsf{u}+O(\mathsf{u}^{2}),$

where $L_{P}$ and $U_{P}$ are the computed LU factors of $P$.

Let $H_{4}=X^{-1}YP^{-1}$ and let $\widehat{H}_{4}$ be the computed value of
$H_{4}$. By (22) and (5.2), we have

$\begin{split}\widehat{H}_{4}&=X^{-1}YJ+O(n)|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\cdot|J|\mathsf{u}\\\
&\quad+O(n)|X^{-1}Y|\cdot|J|(|X|+|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|+|L_{P}|\cdot|U_{P}|)|J|\mathsf{u}\\\
&\quad+O(n)|X^{-1}Y|\cdot|J|\mathsf{u}+O(\mathsf{u}^{2})\\\
&\stackrel{{\scriptstyle(b)}}{{=}}K+O(n)\left(|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\cdot|J|\right.\\\
&\left.+|X^{-1}Y|\cdot|J|(|X|+|Y|\cdot|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|+|L_{P}|\cdot|U_{P}|)|J|\right)\mathsf{u}+O(\mathsf{u}^{2}),\end{split}$
(24)

where $K=X^{-1}YP^{-1}$ and $(b)$ follows from the fact that

$|X^{-1}Y|\cdot|J|=O(|X^{-1}|\cdot|L_{X}|\cdot|U_{X}|\cdot|X^{-1}|\cdot|Y|\cdot|J|).$

The proof is complete by rearranging the terms in (24) and adding the error
terms in (24) and (5.2). ∎

###  5.3. randomized Frobenius inversion

######  Lemma 5.4.

Let $Z=A+iB\in\mathbb{C}^{n\times n}$ be an invertible complex matrix where
$A,B\in\mathbb{R}^{n\times n}$. There exist at most $n$ values of
$\mu\in\mathbb{R}$ such that the real part of $(1+\mu i)Z$ is invertible.

###### Proof.

The real part of $(1+\mu i)(A+iB)$ is $A-\mu B$. We consider the matrix pencil
$D(t)\coloneqq A+tB$ where $t\in\mathbb{C}$. Since $\det(D(t))$ is a
polynomial in $t$ of degree at most $n$ and $\det(D(i))=\det(Z)\neq 0$, $D(t)$
is singular for at most $n$ values of $t\in\mathbb{C}$. In particular, $A-\mu
B$ must be invertible for all but at most $n$ values of $\mu\in\mathbb{R}$. ∎

We remark that for an arbitrary pair $(A,B)$ of matrices, it is possible that
the pencil $D(t)=A+tB$ is singular for all $t\in\mathbb{C}$. The essential
ingredient in the proof of Lemma 5.4 is that $A$ (resp. $B$) is the real
(resp. imaginary) part of an invertible complex matrix. Based on Lemma 5.4, we
obtain Algorithm 5, a randomized version of the Frobenius inversion formula.

Algorithm 5 randomized Frobenius inversion

1:$Z=A+iB\in\operatorname{GL}_{n}(\mathbb{C})$

2:inverse of $Z$

3:randomly generate a real number $\mu\in[0,1]$;

4:compute $X=A-\mu B$ and $Y=\mu A+B$;

5:compute $W=(X+iY)^{-1}$; $\triangleright$ by Algorithm 1

6:set $W_{1}=\operatorname{Re}(W)$ and $W_{2}=\operatorname{Im}(W)$;

7:return $Z^{-1}=(W_{1}-\mu W_{2})+i(\mu W_{1}+W_{2})$;

######  Proposition 5.5.

The running time of Algorithm 5 is $T_{n}+8n^{2}$ where $T_{n}$ denotes the
running time of the Frobenius inversion in Algorithm 1. If we sample
$\mu\in[0,1]$ with respect to a non-atomic probability measure, then Algorithm
5 outputs $Z^{-1}$ correctly, with probability one.

###### Proof.

The statement for the running time of Algorithm 5 is obvious since the running
time of Steps 4 and 6 is $4n^{2}$ respectively. It suffices to prove the
almost sure correctness of Algorithm 5. Since $\mu$ is randomly picked from
$[0,1]$ with a non-atomic probability measure, any finite subset of $[0,1]$ is
a null set. According to Lemma 5.4, $X=A-\mu B$ in Step 4 is invertible with
probability one. Thus Algorithm 1 is applicable to $X+iY$ and we have

$(W_{1}-\mu W_{2})+i(\mu W_{1}+W_{2})=(1+\mu i)W=(1+\mu i)(X+iY)^{-1}=(1+\mu
i)\left(Z(1+\mu i)\right)^{-1}=Z^{-1}.$

∎

It is clear that $Z^{-1}=(MZ)^{-1}M$ for any invertible matrix
$M\in\mathbb{C}^{n\times n}$. By suitably choosing $M$, we obtain variants of
the Frobenius inversion formula [64, 16]. Unfortunately, these variants are
more computationally expensive. Indeed, if we multiply real matrices by the
usual method. Then the one given in [16] costs $6n^{3}$ multiplications and
the one in [64] costs $7n^{3}+2n^{2}$ multiplications. As a comparison, the
Frobenius formula only costs $5n^{3}$ multiplications. On the other side, the
Frobenius inversion formula (and its variants) is not applicable to $Z$ in a
subvariety of $\operatorname{GL}_{n}(\mathbb{C})$. To resolve the issue, we
may randomly generate $M$ so that the inversion formula is applicable to $MZ$
with probability one. However, if $M$ is a dense matrix then computing $MZ$
and $(MZ)^{-1}M$ costs $n^{3}$ multiplications, respectively. Thus the
resulting randomized algorithm would cost $7n^{3}$ multiplications [64].
According to Proposition 5.5, it is sufficient to take $M=(1+\mu i)I_{n}$
where $\mu$ is randomly picked from $[0,1]$. More importantly, such a choice
of $M$ reduces the number of multiplications from $7n^{3}$ to $5n^{3}+8n^{2}$.

##  6\. Hermitian positive definite matrix inversion

In this section, we consider a special class of matrices that occur frequently
in practice. We assume that $Z=A+iB$ is Hermitian positive definite and
$Z\in\mathcal{S}_{1}$ According to Algorithm 1, we need to compute $A^{-1}B$,
$BA^{-1}B$, $(A+BA^{-1}B)^{-1}$ and $A^{-1}B(A+BA^{-1}B)^{-1}$. Since $Z$ is
Hermitian positive definite, we claim (I’ll add a proof later) that both $A$
and $A+BA^{-1}B$ are symmetric positive definite. Moreover, $B$ is skew-
symmetric. Therefore we may compute $A^{-1}B$ and $(A+BA^{-1}B)^{-1}$ by
Cholesky decomposition. Suppose that $A=U^{{\scriptscriptstyle\mathsf{T}}}U$
is the Cholesky decomposition of $A$. Then we have

$\begin{split}&A^{-1}B=U^{-1}(U^{-1})^{\scriptscriptstyle\mathsf{T}}B,\\\
&BA^{-1}B=BU^{-1}(U^{-1})^{\scriptscriptstyle\mathsf{T}}B=-\left[(U^{-1})^{{\scriptscriptstyle\mathsf{T}}}B\right]^{{\scriptscriptstyle\mathsf{T}}}\left[(U^{-1})^{\scriptscriptstyle\mathsf{T}}B\right],\\\
&(A+BA^{-1}B)^{-1}=\left(A-\left[(U^{-1})^{{\scriptscriptstyle\mathsf{T}}}B\right]^{{\scriptscriptstyle\mathsf{T}}}\left[(U^{-1})^{\scriptscriptstyle\mathsf{T}}B\right]\right)^{-1}=V^{-1}(V^{-1})^{\scriptscriptstyle\mathsf{T}},\\\
&A^{-1}B(A+BA^{-1}B)^{-1}=A^{-1}BV^{-1}(V^{-1})^{\scriptscriptstyle\mathsf{T}},\end{split}$
(25)

where $A+BA^{-1}B=V^{\scriptscriptstyle\mathsf{T}}V$ is the Cholesky
decomposition of
$A-\left[(U^{-1})^{{\scriptscriptstyle\mathsf{T}}}B\right]^{{\scriptscriptstyle\mathsf{T}}}\left[(U^{-1})^{\scriptscriptstyle\mathsf{T}}B\right]$.

######  Lemma 6.1.

Let $Z=A+iB$ be an $n\times n$ Hermitian positive definite matrix where
$A,B\in\mathbb{R}^{n\times n}$. Then the following properties hold:

  1. (a)

$A$ is symmetric positive definite and $B$ is skew-symmetric.

  2. (b)

$A+BA^{-1}B$ is positive definite.

In particular, if we denote by $\mathbb{H}_{n}^{++}$ the set of all $n\times
n$ Hermitian positive definite matrix. Then we must have

$\mathbb{H}_{n}^{++}\subseteq\mathcal{S}_{1}\coloneqq\\{Z=A+iB\in\operatorname{GL}_{n}(\mathbb{\mathbb{C}}):A,A+BA^{-1}B\in\operatorname{GL}_{n}(\mathbb{R})\\}.$

###### Proof.

We notice that $A=(Z+\overline{Z})/2$ and $B=(Z-\overline{Z})/2i$. Thus $A$
(resp. $B$) is symmetric (resp. skew-symmetric) since $Z$ is Hermitian. For
any $x\in\mathbb{R}^{n}$, we have

$x^{\scriptscriptstyle\mathsf{T}}Ax=\frac{x^{\ast}(Z+\overline{Z})x}{2}=\frac{x^{\ast}Zx}{2}+\overline{\left(\frac{x^{\ast}Zx}{2}\right)}=x^{\ast}Zx\geq
0$

and the equality holds if and only if $x=0$. This proves the positive
definiteness of $A$.

For each $z\in\mathbb{C}^{n}$, we have

$z^{\ast}\overline{Z}z=\overline{\overline{z}^{\ast}Z\overline{z}}\geq 0,$

and the equality holds if and only if $z=0$ since $Z$ is positive definite.
This implies that $\overline{Z}$ is Hermitian positive definite. We also
observe that
$A^{-\frac{1}{2}}ZA^{-\frac{1}{2}}=I_{n}+iA^{-\frac{1}{2}}BA^{-\frac{1}{2}}\succ
0$ and moreover

$A+BA^{-1}B=A^{\frac{1}{2}}\left(I_{n}+(A^{-\frac{1}{2}}BA^{-\frac{1}{2}})(A^{-\frac{1}{2}}BA^{-\frac{1}{2}})\right)A^{\frac{1}{2}}.$

Therefore, it is sufficient to prove (b) for $A=I_{n}$. In this case, we have

$I_{n}+iB\succ 0,\quad I_{n}-iB\succ 0,$

from which we may conclude that $I_{n}+B^{2}\succ 0$ since
$I_{n}+B^{2}=(I_{n}-iB)^{\frac{1}{2}}(I_{n}+iB)(I_{n}-iB)^{\frac{1}{2}}$ ∎

We remark that there is another way to see the positive definiteness of
$A+BA^{-1}B$ if we assume a priori it is invertible. Indeed, Lemma 4.1 impies
that $(A+BA^{-1}B)^{-1}$ is the real part of $Z^{-1}$. By assumption, $Z$ is
positive definite, so are $Z^{-1}$ and its real part $(A+BA^{-1}B)^{-1}$. This
implies that $A+BA^{-1}B$ must be positive definite as well.

We give two implementations of equation (25) which have the same time
complexity.

Algorithm 6 First Variant of Frobenius inversion 1

1:$Z=A+iB\in\mathcal{S}_{1}\cup\mathcal{S}_{2}$

2:inverse of $Z$

3:if $Z\in\mathcal{S}_{1}$ then

4: set $X=A,Y=B$;

5:else if $Z\in\mathcal{S}_{2}$ then

6: set $X=B,Y=A$;

7:end if

8:compute Cholesky decomposition of $X=U^{\scriptscriptstyle\mathsf{T}}U$;

9:compute $K_{1}=(U^{\scriptscriptstyle\mathsf{T}})^{-1}Y$;

10:compute $K_{2}=U^{-1}K_{1}$;

11:compute $K_{3}=K_{1}^{\scriptscriptstyle\mathsf{T}}K_{1}$;

12:compute $K_{4}=X-K_{3}$;

13:compute Cholesky decomposition of
$K_{4}=V^{\scriptscriptstyle\mathsf{T}}V$;

14:compute $K_{5}=V^{-1}$;

15:compute $K_{6}=K_{5}K_{5}^{\scriptscriptstyle\mathsf{T}}$;

16:compute $K_{7}=K_{2}K_{6}$;

17:if $Z\in\mathcal{S}_{1}$ then return $Z^{-1}=K_{6}-iK_{7}$;

18:else if $Z\in\mathcal{S}_{2}$ then return $Z^{-1}=K_{7}-iK_{6}$;

19:end if

Algorithm 7 Second Variant of Frobenius inversion 1

1:$Z=A+iB\in\mathcal{S}_{1}\cup\mathcal{S}_{2}$

2:inverse of $Z$

3:if $Z\in\mathcal{S}_{1}$ then

4: set $X=A,Y=B$;

5:else if $Z\in\mathcal{S}_{2}$ then

6: set $X=B,Y=A$;

7:end if

8:compute Cholesky decomposition of $X=U^{\scriptscriptstyle\mathsf{T}}U$;

9:compute $K_{1}=(U^{\scriptscriptstyle\mathsf{T}})^{-1}Y$;

10:compute $K_{2}=U^{-1}K_{1}$;

11:compute $K_{3}=K_{1}^{\scriptscriptstyle\mathsf{T}}K_{1}$;

12:compute $K_{4}=X-K_{3}$;

13:compute Cholesky decomposition of
$K_{4}=V^{\scriptscriptstyle\mathsf{T}}V$;

14:compute $K_{5}=(V^{\scriptscriptstyle\mathsf{T}})^{-1}$;

15:compute $K_{6}=V^{-1}K_{5}$;

16:compute $K_{7}=K_{2}K_{6}$;

17:if $Z\in\mathcal{S}_{1}$ then return $Z^{-1}=K_{6}-iK_{7}$;

18:else if $Z\in\mathcal{S}_{2}$ then return $Z^{-1}=K_{7}-iK_{6}$;

19:end if

Note that the only differences between these two implementations are line $12$
and line $13$. In the first implementation, computing $K_{5}=V^{-1}$ takes
$\Theta(n^{3})$ flops since $V$ is upper triangular and computing
$K_{6}=K_{5}K_{5}^{\scriptscriptstyle\mathsf{T}}$ takes $\Theta(n^{3})$ flops
since $K_{5}$ is upper triangular. In the second implementation, computing
$K_{5}=(V^{\scriptscriptstyle\mathsf{T}})^{-1}$ takes $\Theta(n^{3})$ flops
since $V^{\scriptscriptstyle\mathsf{T}}$ is lower triangular and computing
$K_{6}=V^{-1}K_{5}$ takes $\Theta(n^{3})$ flops since solving the triangular
system $VK_{6}=K_{5}$ takes $n$ back substitutions, each of which takes
$\Theta(n^{2})$ flops. Thus, Algorithm 6 and Algorithm 7 have the same flop
count. In practice, Algorithm 6 is more efficient than Algorithm 7 when
implemented on Matlab.

###  6.1. variant of Frobenius inversion vs matrix inversion via Cholesky
decomposition

In this section, we compare the variant of Frobenius inversion to matrix
inversion algorithm using Cholesky decomposition. When the input matrix $Z$ is
hermitian positive definite, Matlab’s inversion function exploits Cholesky
decomposition to compute $Z^{-1}$. Similar to what we did to the variant of
Frobenius inversion, we give two implementations of Cholesky decompositions
that have the same time complexity.

Algorithm 8 First matrix inversion via Cholesky decomposition

1:$A\in\operatorname{GL}_{n}(\Bbbk)$

2:inverse of $A$

3:compute Cholesky decomposition of $A=U^{\scriptscriptstyle\mathsf{T}}U$;

4:compute $K=U^{-1}$;

5:compute $X=KK^{\scriptscriptstyle\mathsf{T}}$;

6:return $X$;

Algorithm 9 Second matrix inversion via Cholesky decomposition

1:$A\in\operatorname{GL}_{n}(\Bbbk)$

2:inverse of $A$

3:compute Cholesky decomposition of $A=U^{\scriptscriptstyle\mathsf{T}}U$;

4:compute $(U^{\scriptscriptstyle\mathsf{T}})^{-1}$;

5:solve for $X$ from $UX=(U^{\scriptscriptstyle\mathsf{T}})^{-1}$;

6:return $X$;

These two implementations have the same flop counts by the same reasoning as
we did in the variant of Frobenius inversion case. In practice, Algorithm 8
runs faster than Algorithm 9.

For the sake of speed analysis, we will compare Algorithm 7 to Algorithm 9.
However, our results hold for Algorithm 6 versus Algorithm 8 as well since
Algorithm 7 has the same flop count as Algorithm 6 and Algorithm 9 has the
same flop count as Algorithm 8. In order to analyze the speed of Algorithm 7,
we present an algorithm for Cholesky decomposition in Algorithm 10 [56].

Algorithm 10 Cholesky decomposition

1:$Z\in\mathbb{C}^{n}$ is hermitian positive definite

2:inverse of $Z$

3:let $Z$ be the upper triangular part of $Z$;

4:for $k=1$ to $n$ do

5: $Z[1:k,k]=(Z[1:k-1,1:k-1]^{*})^{-1}Z[1:k-1,k]$;

6: $Z[k,k]=\sqrt{Z[k,k]-Z[1:k-1,k]^{*}Z[1:k-1,k]}$;

7:end for

8:return $Z$;

Note that Algorithm 10 uses roughly the same number of additions and
multiplications since its main operation is forward substitution: solve for
$X$ in $Z[1:k-1,1:k-1]^{*}X=Z[1:k-1,k]$, where $Z[1:k-1,1:k-1]$ is upper
triangular since we take $Z$ to be the upper triangular part of $Z$ in the
first line of Algorithm 10.

Before we proceed, we fix some notations. Let $\mathcal{A}$ be an algorithm
for real matrix multiplication. We denote by
$T_{\operatorname{mult}}^{\mathcal{A}}(n)$ the average running time of
$\mathcal{A}$ on pairs of $n\times n$ real matrices. In addition, we let
$T_{\operatorname{pinv}}^{\mathcal{A}}(n)$ be the average running time of
Algorithm 9 on invertible $n\times n$ real symmetric positive definite
matrices, in which real matrix multiplications are computed by $\mathcal{A}$.
Now with these notations, we are ready to present our threshold theorem.

######  Theorem 6.2 (threshold).

Let $\mathcal{A}$ be the usual algorithm for real matrix multiplication. Then
Algorithm 7 is asymptotically faster than Algorithm 9 over $\mathbb{C}$ if and
only if

$\lim_{n\to\infty}\left(T_{\operatorname{pinv}}^{\mathcal{A}}(n)/T_{\operatorname{mult}}^{\mathcal{A}}(n)\right)>7/6.$

###### Proof.

We first show that the running time of Algorithm 7 is dominated by
$2T_{\operatorname{pinv}}^{\mathcal{A}}(n)+7T_{\operatorname{mult}}^{\mathcal{A}}(n)/3$.

First, we show that the first three steps in Algorithm 7 take
$T_{\operatorname{pinv}}^{\mathcal{A}}(n)+T_{\operatorname{mult}}^{\mathcal{A}}(n)/3$
time. Recall that the first three steps in Algorithm 7 are computing the
Cholesky decomposition of $X=U^{\scriptscriptstyle\mathsf{T}}U$, computing
$K_{1}=(U^{\scriptscriptstyle\mathsf{T}})^{-1}Y,$ and computing
$K_{2}=U^{-1}K_{1}$. In terms of speed, the only difference between these
three steps and Algorithm 9 is in step two, where Algorithm 9 solve for $X$ in
$U^{\scriptscriptstyle\mathsf{T}}X=I$ whereas Algorithm 7 solve for $X$ in
$U^{\scriptscriptstyle\mathsf{T}}X=Y$. Solving for $X$ in
$U^{\scriptscriptstyle\mathsf{T}}X=I$ only takes $\Theta(n^{3}/3)$ flops
whereas solving for $X$ in $U^{\scriptscriptstyle\mathsf{T}}X=Y$ takes
$\Theta(n^{3})$ flops if $Y$ does not have special structures. Thus, Algorithm
7 takes $2\Theta(n^{3}/3)$ more flops which is
$T_{\operatorname{mult}}^{\mathcal{A}}(n)/3$ since the usual matrix
multiplication algorithm takes $\Theta(2n^{3})$ flops.

Then, note that the running time of the rest of Algorithm 7 is dominated by
two real matrix multiplication and one real matrix inversion. Thus, the
running time of Algorithm 7 is dominated by
$2T_{\operatorname{pinv}}^{\mathcal{A}}(n)+7T_{\operatorname{mult}}^{\mathcal{A}}(n)/3$.

Next, we consider the running time of Algorithm 9 for complex matrices. We
prove that the running time of Algorithm 9 is dominated by
$4T_{\operatorname{pinv}}^{\mathcal{A}}(n)$. Note that the complex addition
takes $2$ real flops and the complex multiplication takes $6$ real flops. In
Algorithm 9, there are “roughly” the same number of additions and
multiplications. Algorithm 9 contains roughly three operations: computing
Cholesky factorization, computing backward substitution, and computing the
inverse of an lower triangular matrix. As we have observed in Algorithm 10,
computing Cholesky factorization requires roughly the same number of additions
and multiplications. Computing backward (or forward) substitution also
requires roughly the same number of additions and multiplications. Finally,
according to Method $1$ in page $263$ of [26], inverting a triangular matrix
can be done by a sequence of forward substitutions. Thus, inverting a
triangular matrix also requires roughly the same number of additions and
multiplications. Therefore, the running time of Algorithm 9 over $\mathbb{C}$
is dominated by $4T_{\operatorname{pinv}}^{\mathcal{A}}(n)$.

Now, Algorithm 7 is faster than Algorithm 9 if and only if for $n$
sufficiently large,

$4T_{\operatorname{inv}}^{\mathcal{A}}(n)>2T_{\operatorname{inv}}^{\mathcal{A}}(n)+7T_{\operatorname{mult}}^{\mathcal{A}}(n)/3,$

which is equivalent to

$\lim_{n\to\infty}\left(\frac{T_{\operatorname{inv}}^{\mathcal{A}}(n)}{T_{\operatorname{mult}}^{\mathcal{A}}(n)}\right)>7/6.\qed$

###  6.2. rounding error analysis

In this section, we give a rounding error analysis of Algorithm 6. Before
giving the result, we first introduce some notations. For each matrix $X$, we
use $\widehat{X}$ to denote the computed value of $X$. Let

$\begin{split}\Delta_{1}^{\prime}&=n|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}|\cdot|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}\\\
\Delta_{2}^{\prime}&=n|\widehat{U}^{-1}|\cdot|\widehat{U}|\cdot|\widehat{U}^{-1}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}\\\
\Delta_{3}^{\prime}&=|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\Delta_{1}^{\prime}+\Delta_{1}^{\prime{\scriptscriptstyle\mathsf{T}}}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+n|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}\\\
\Delta_{4}^{\prime}&=|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}\Delta
XX^{-1}Y|+\Delta_{3}^{\prime}\\\
\Delta_{5}^{\prime}&=\Delta_{4}^{\prime}+(|X|+|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y|)\mathsf{u}\\\
\Delta_{6}^{\prime}&=O(n)|\widehat{V}^{\scriptscriptstyle\mathsf{T}}||\widehat{V}|\mathsf{u}\\\
\Delta_{7}^{\prime}&=O(n)|\widehat{V}^{-1}|\cdot(|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{V}^{\scriptscriptstyle\mathsf{T}}|+|\widehat{V}|\cdot|\widehat{V}^{-1}|+I_{n})\cdot|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}\\\
\Delta_{8}^{\prime}&=\Delta_{5}^{\prime}+\Delta_{6}^{\prime}\\\
\Delta_{9}^{\prime}&=|(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}|\Delta_{8}^{\prime}|(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}|+\Delta_{7}^{\prime}.\end{split}$

######  Theorem 6.3.

Let $Z=A+iB\in\mathcal{S}_{1}\cup\mathcal{S}_{2}$ be a hermitian positive
definite matrix and $W$ be the computed inverse of $Z$ using Algorithm 6. Let
$X=A$ and $Y=B$ if $Z\in\mathcal{S}_{1}$ and let $X=B$ and $Y=A$ if
$Z\in\mathcal{S}_{2}$. Then we have

$\begin{split}|Z^{-1}-W|&\leq(|K_{2}|+1)\Delta_{9}^{\prime}+\Delta_{2}^{\prime}|K_{6}|+n|K_{2}||K_{6}|\mathsf{u}+O(\mathsf{u}^{2}).\end{split}$
(26)

The proof of Theorem 6.3 relies on the following facts in [26]:

  * •

Let $A\in\mathbb{R}^{n\times n}$ be a symmetric positive definite matrix. Let
$\widehat{R}$ be the computed Cholesky factor of $A$. Then

$\widehat{R}^{\scriptscriptstyle\mathsf{T}}\widehat{R}=A+\Delta A,$ (27)

where $|\Delta
A|\leq(n+1)|\widehat{R}^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{R}|\mathsf{u}+O(\mathsf{u}^{2})$.

  * •

Let $T\in\mathbb{R}^{n\times n}$ be a triangular matrix and
$b\in\mathbb{R}^{n}$. Let $\hat{x}$ be the solution to $Tx=b$ obtained by
forward/backward substitution. Then

$(T+\Delta T)\hat{x}=b,$ (28)

where $|\Delta T|\leq n|T|\mathsf{u}+O(u^{2})$.

  * •

Let $T$ be a triangular matrix. Let $\widehat{X}$ be the computed inverse of
$T$. Then

$|\widehat{X}-T^{-1}|\leq
O(n)|T^{-1}|\cdot|T|\cdot|T^{-1}|\mathsf{u}+O(\mathsf{u}^{2}).$ (29)

Now, we are able to prove Theorem 6.3. In the following, we denote by
$\widehat{X}$ the the computed value of a matrix $X$.

###### proof of Theorem 6.3.

We introduce some error terms $\Delta_{i}$ in the proof such that
$|\Delta_{i}|\leq\Delta_{i}^{\prime}+O(\mathsf{u}^{2})$ for all $i$. Let
$\widehat{U}$ be the computed Cholesky factor of $X$. By equation (27),

$\widehat{U}^{\scriptscriptstyle\mathsf{T}}\widehat{U}=X+\Delta
X,\quad\quad|\Delta X|\leq
O(n)|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{U}|\mathsf{u}+O(\mathsf{u}^{2}).$
(30)

Then, by equation (28),

$(\widehat{U}^{\scriptscriptstyle\mathsf{T}}+\Delta\widehat{U}^{\scriptscriptstyle\mathsf{T}})\widehat{K}_{1}=Y,\quad\quad|\Delta\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\leq
n|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}+O(\mathsf{u}^{2}).$
(31)

This implies that

$\begin{split}\widehat{K}_{1}&=(\widehat{U}^{\scriptscriptstyle\mathsf{T}}+\Delta\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y\\\
&\stackrel{{\scriptstyle(a)}}{{=}}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y-(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}\Delta\widehat{U}^{\scriptscriptstyle\mathsf{T}}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+O(\mathsf{u}^{2})\\\
&=(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+\Delta_{1},\end{split}$
(32)

where $(a)$ follows from Lemma 5.3 and
$\Delta_{1}\vcentcolon=-(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}\Delta\widehat{U}^{\scriptscriptstyle\mathsf{T}}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+O(\mathsf{u}^{2})$.
Note that

$|\Delta_{1}|\leq
n|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}|\cdot|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}+O(\mathsf{u}^{2}),$
(33)

by equation (31). Similarly,

$(\widehat{U}+\Delta\widehat{U})\widehat{K}_{2}=\widehat{K}_{1},\quad\quad|\Delta\widehat{U}|\leq
n|\widehat{U}|\mathsf{u}+O(\mathsf{u}^{2}).$ (34)

This implies that

$\begin{split}\widehat{K}_{2}&=(\widehat{U}+\Delta\widehat{U})^{-1}\widehat{K}_{1}\\\
&\stackrel{{\scriptstyle(b)}}{{=}}\widehat{U}^{-1}\widehat{K}_{1}-\widehat{U}^{-1}\Delta\widehat{U}\widehat{U}^{-1}\widehat{K}_{1}+O(\mathsf{u}^{2})\\\
&\stackrel{{\scriptstyle(c)}}{{=}}\widehat{U}^{-1}\widehat{K}_{1}-\widehat{U}^{-1}\Delta\widehat{U}\widehat{U}^{-1}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+O(\mathsf{u}^{2})\\\
&=\widehat{U}^{-1}\widehat{K}_{1}+\Delta_{2},\end{split}$ (35)

where $(b)$ follows from Lemma 5.3, $(c)$ follows from equation (32), and

$\Delta_{2}\vcentcolon=-\widehat{U}^{-1}\Delta\widehat{U}\widehat{U}^{-1}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+O(\mathsf{u}^{2}).$

Note that

$|\Delta_{2}|\leq
n|\widehat{U}^{-1}|\cdot|\widehat{U}|\cdot|\widehat{U}^{-1}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}+O(\mathsf{u}^{2})$
(36)

by equation (34). Then, by equation (19),

$\begin{split}\widehat{K}_{3}&=\widehat{K}_{1}^{\scriptscriptstyle\mathsf{T}}\widehat{K}_{1}+n|\widehat{K}_{1}^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{K}_{1}|\mathsf{u}+O(\mathsf{u}^{2})\\\
&\stackrel{{\scriptstyle(d)}}{{=}}Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}\Delta_{1}+\Delta_{1}^{\scriptscriptstyle\mathsf{T}}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+n|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+\Delta_{3}\\\
&\stackrel{{\scriptstyle(e)}}{{=}}Y^{\scriptscriptstyle\mathsf{T}}(X+\Delta
X)^{-1}Y+\Delta_{3}\\\
&\stackrel{{\scriptstyle(f)}}{{=}}Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}\Delta
XX^{-1}Y+\Delta_{3}\\\
&=Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y+\Delta_{4},\end{split}$ (37)

where $(d)$ follows from equation (32), $(e)$ follows from equation (30),
$(f)$ follows from Lemma 5.3,

$\Delta_{3}\vcentcolon=Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}\Delta_{1}+\Delta_{1}^{\scriptscriptstyle\mathsf{T}}(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y+n|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}+O(\mathsf{u}^{2})$

and

$\Delta_{4}\vcentcolon=-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}\Delta
XX^{-1}Y+\Delta_{3}.$

Note that

$\begin{split}|\Delta_{4}|&\stackrel{{\scriptstyle(a)}}{{\leq}}O(n)|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}|\cdot|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{U}|\cdot|X^{-1}Y|\mathsf{u}+n|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}|\cdot|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}\\\
&\quad+n|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot|\widehat{U}|\cdot|\widehat{U}^{-1}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}+n|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=O(n)(|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}|\cdot|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{U}|\cdot|X^{-1}Y|+|Y^{\scriptscriptstyle\mathsf{T}}\widehat{U}^{-1}|\cdot(|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}|\cdot|\widehat{U}^{\scriptscriptstyle\mathsf{T}}|+|\widehat{U}|\cdot|\widehat{U}^{-1}|+I_{n})\cdot|(\widehat{U}^{\scriptscriptstyle\mathsf{T}})^{-1}Y|)\mathsf{u}+O(\mathsf{u}^{2}),\end{split}$
(38)

where $(a)$ follows from equation (30) and equation (33). Then, by equation
(21),

$\begin{split}\widehat{K}_{4}&=X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y+\Delta_{4}+(|X|+|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y|)\mathsf{u}+O(\mathsf{u}^{2})\\\
&=X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y+\Delta_{5},\end{split}$ (39)

where

$\Delta_{5}\vcentcolon=\Delta_{4}+(|X|+|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y|)\mathsf{u}+O(\mathsf{u}^{2}).$
(40)

Let $\widehat{V}$ be the computed Cholesky factor of $\widehat{K}_{4}$. By
equation (27),

$\widehat{V}^{\scriptscriptstyle\mathsf{T}}\widehat{V}=\widehat{K}_{4}+\Delta_{6},\quad\quad|\Delta_{6}|\leq
O(n)|\widehat{V}^{\scriptscriptstyle\mathsf{T}}||\widehat{V}|\mathsf{u}+O(\mathsf{u}^{2}).$
(41)

Then, by equation (29),

$\widehat{K}_{5}=\widehat{V}^{-1}+O(n)|\widehat{V}^{-1}|\cdot|\widehat{V}|\cdot|\widehat{V}^{-1}|\mathsf{u}+O(\mathsf{u}^{2}).$
(42)

Then, by equation (19),

$\begin{split}\widehat{K}_{6}&=\widehat{V}^{-1}(\widehat{V}^{\scriptscriptstyle\mathsf{T}})^{-1}+O(n)|\widehat{V}^{-1}|\cdot|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{V}^{\scriptscriptstyle\mathsf{T}}|\cdot|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}\\\
&\quad+O(n)|\widehat{V}^{-1}|\cdot|\widehat{V}|\cdot|\widehat{V}^{-1}||(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}+n|\widehat{V}^{-1}||(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=\widehat{V}^{-1}(\widehat{V}^{\scriptscriptstyle\mathsf{T}})^{-1}+O(n)|\widehat{V}^{-1}|\cdot(|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{V}^{\scriptscriptstyle\mathsf{T}}|+|\widehat{V}|\cdot|\widehat{V}^{-1}|+I_{n})\cdot|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=\widehat{V}^{-1}(\widehat{V}^{\scriptscriptstyle\mathsf{T}})^{-1}+\Delta_{7}\\\
&\stackrel{{\scriptstyle(a)}}{{=}}(\widehat{K}_{4}+\Delta_{6})^{-1}+\Delta_{7}\\\
&\stackrel{{\scriptstyle(b)}}{{=}}(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y+\Delta_{8})^{-1}+\Delta_{7}\\\
&\stackrel{{\scriptstyle(c)}}{{=}}(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}+(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}\Delta_{8}(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}+\Delta_{7}+O(\mathsf{u}^{2})\\\
&=K_{6}+\Delta_{9},\end{split}$ (43)

where $(a)$ follows from equation (41), $(b)$ follows from equation (39),
$(c)$ follows from Lemma 5.3, and

$\begin{split}\Delta_{7}&\vcentcolon=O(n)|\widehat{V}^{-1}|\cdot(|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{V}^{\scriptscriptstyle\mathsf{T}}|+|\widehat{V}|\cdot|\widehat{V}^{-1}|+I_{n})\cdot|(\widehat{V}^{-1})^{\scriptscriptstyle\mathsf{T}}|\mathsf{u}+O(\mathsf{u}^{2})\\\
\Delta_{8}&\vcentcolon=\Delta_{5}+\Delta_{6}\\\
\Delta_{9}&\vcentcolon=(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}\Delta_{8}(X-Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y)^{-1}+\Delta_{7}+O(\mathsf{u}^{2})\end{split}$

Note that

$\begin{split}|\Delta_{9}|&\leq|K_{6}|\cdot|\Delta_{8}|\cdot|K_{6}|+|\Delta_{7}|+O(\mathsf{u}^{2})\\\
&\leq|K_{6}|\cdot(|\Delta_{5}|+|\Delta_{6}|)\cdot|K_{6}|+|\Delta_{7}|+O(\mathsf{u}^{2})\\\
&\stackrel{{\scriptstyle(a)}}{{\leq}}|K_{6}|\cdot(|\Delta_{4}|+(|X|+|Y^{\scriptscriptstyle\mathsf{T}}X^{-1}Y|)\mathsf{u}+O(n)|\widehat{V}^{\scriptscriptstyle\mathsf{T}}|\cdot|\widehat{V}|\mathsf{u})\cdot|K_{6}|+|\Delta_{7}|+O(\mathsf{u}^{2}),\end{split}$
(44)

where $(a)$ follows from equation (41) and equation (40). Then, by equation
(19),

$\begin{split}\widehat{K}_{7}&=\widehat{K}_{2}\widehat{K}_{6}+n|\widehat{K}_{2}|\cdot|\widehat{K}_{6}|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=(K_{2}+\Delta_{2})(K_{6}+\Delta_{9})+n|K_{2}||K_{6}|\mathsf{u}+O(\mathsf{u}^{2})\\\
&=K_{7}+|K_{2}||\Delta_{9}|+|\Delta_{2}||K_{6}|+n|K_{2}||K_{6}|\mathsf{u}+O(\mathsf{u}^{2}).\end{split}$
(45)

The result then follows from equation (45) and equation (43). ∎

##  7\. Experiments

In this section, we conduct some experiments to compare the Frobenius
inversion (Algorithm 1) with the inversion via LU decomposition (Algorithm 3).
We compare both the speed and accuracy of the two methods by randomly
generated examples. Moreover, we test Algorithms 1 and 3 on matrix sign
function, Sylvester equations, Lyapunov equations and polar decomposition.
These experiments show that the Frobenius inversion is more efficient than the
inversion via LU decomposition, confirming the threshold Theorem 5.1. It is
also clear from these experiments that the Frobenius inversion is a bit less
accuracy than the inversion via LU decomposition. However, comparing with the
increase of efficiency, such a decrease of accuracy is ignorable. Finally, we
compare the variant of Frobenius inversion (Algorithm 6) to matrix inversion
using Cholesky decomposition (Algorithm 8) on hermitian positive matrices. We
show that the variant of Frobenius inversion is more efficient than matrix
inversion using Cholesky decomposition with little loss in accuracy. This in
turn leads to a more efficient procedure to process MIMO radios.

###  7.1. efficiency

Let $3600\leq n\leq 6000$ be a positive integer. We first generate two real
matrices $A,B\in\mathbb{R}^{n\times n}$ where elements of $A$ and $B$ are
generated uniformly and randomly from $(0,1)$. Then we take $X=A+iB$. For each
dimension $n$, we generate $10$ random matrices by the above procedure.

For each randomly generated matrix $X$, we compute the inverse of $X$ by
Algorithm 1 and Algorithm 3111In MATLAB, this is simply the command
$X\backslash\operatorname{eye}(n)$). respectively. For each fixed $n$, we
average the computation time for the $10$ random instances. In Figure 1, we
exhibit the computation time of the two algorithms versus the logarithmic
dimension of matrices. It is clear from Figure 1 that Frobenius inversion is
faster than the inversion algorithm via LU decomposition.

Figure 1. comparison of efficiency

.

###  7.2. accuracy

Next we compare the accuracy of the two algorithms. For any complex invertible
matrix $X=A+iB$, where $A,B\in\mathbb{R}^{n\times n}$, we let $\widehat{X}$ be
the computed inverse of $X$. We assess the accuracy of $\widehat{X}$ in terms
of its left and right relative residual defined as

$\operatorname{res}_{L}(X,\widehat{X})\vcentcolon=\frac{\lVert\widehat{X}X-I\rVert_{\max}}{\lVert
X\rVert_{\max}\lVert\widehat{X}\rVert_{\max}}\quad\textrm{and}\quad\operatorname{res}_{R}(X,\widehat{X})\vcentcolon=\frac{\lVert
X\widehat{X}-I\rVert_{\max}}{\lVert
X\rVert_{\max}\lVert\widehat{X}\rVert_{\max}},$ (46)

where $\lVert\cdot\rVert_{\max}$ is defined as

$\lVert A+iB\rVert_{\max}\vcentcolon=\max(\lVert A\rVert_{\max},\lVert
B\rVert_{\max})\vcentcolon=\max\biggl{(}\max_{i,j\in[n]}|a_{ij}|,\max_{i,j\in[n]}|b_{ij}|\biggr{)},$
(47)

where $a_{ij}$ and $b_{ij}$ are the $i,j$th entries of $A$ and $B$
respectively. To avoid numerical issues, we only test the two algorithms on
well-conditioned matrices ( i.e., matrices whose condition number is $10$)
whose dimension goes from $2$ to $4096$. We generate a random matrices whose
condition number is $\kappa$ by the procedure that follows. We first generate
$(n-2)$ integers from $1$ to $(\kappa-1)$ uniformly and randomly. Then, we
consider the $n$ by $n$ diagonal matrix $D$ whose diagonal elements are the
$(n-2)$ integers we just generate together with $1$ and $\kappa$. We next
multiply the Hadamard matrix $H$ of order $n$ to the left of $D$ and the
transpose of $H$ to the right of $D$ to obtain
$HDH^{\scriptscriptstyle\mathsf{T}}$. Lastly, we scale this matrix by some
constant so that its Frobenius norm is one. Since the smallest and largest
singular values of this matrix are $c$ and $\kappa c$ for some
$c\in\mathbb{R}$, the condition number of this matrix is $\kappa$. In our
experiment, we take $\kappa=10$. We note that there is no specific sparsity
patterns on the matrices we generate so that our experiment is not biased
toward a certain class of matrices with special sparsity patterns.

In Figure 2, we compare the left and right relative residuals of $\widehat{X}$
computed by the two algorithms. It is obvious from the right plot in Figure 2
that the Frobenius inversion is slightly less accurate than the inversion via
LU decomposition. However, the difference is so small that one can ignore it
for practical purposes.

Figure 2. comparison of relative residuals

###  7.3. matrix sign function

In this subsection, we apply matrix inversion algorithms to compute the matrix
sign function. The tatrix sign function can be used to solve a wide range of
problems such as algebraic Riccati equation [51], Sylvester equation [25, 51],
polar decomposition [25], and spectral decomposition [2, 3, 5, 29, 39]. We
first recall the definition of the matrix sign function. Given a matrix
$A\in\mathbb{C}^{n\times n}$, we write $A=ZJZ^{-1}$ where
$J=\begin{bmatrix}J_{1}&0\\\ 0&J_{2}\end{bmatrix}$ is the Jordan canonical
form of $A$ such that eigenvalues of $A$ in the diagonal of $J_{1}$ (resp.
$J_{2}$) have negative (resp. positive) real parts. The _matrix sign function_
is defined to be

$\operatorname{sign}(A)=Z\begin{bmatrix}-I&0\\\ 0&I\end{bmatrix}Z^{-1}.$

The matrix sign function $\operatorname{sign}(A)$ can be computed by Newton
iterations [27, 51]

$X_{t+1}=\frac{1}{2}(X_{t}+X_{t}^{-1}),\quad t\in\mathbb{N},$

with $X_{0}=A$. In each iteration, we compute $X_{t}^{-1}$ either by Frobenius
inversion or by inversion via LU decomposition. This gives us two computed
value of $\operatorname{sign}(A)$ for each matrix $A$. We measure the progress
in each Newton iteration by the relative change in $X_{t}$:

$\delta_{t}\vcentcolon=\frac{\lVert X_{t}-X_{t-1}\rVert_{\max}}{\lVert
X_{t}\rVert_{\max}},$

where $\lVert\cdot\rVert_{\max}$ is defined in (47). We stop the Newton
iterations when either $\delta_{t}\leq 10^{-3}$ or $t\geq 100$ is achieved.

We compute sign function for random square matrices whose dimension $n$ is an
even number between $2100$ and $4000$. To generate such a random matrix, we
first generate an $n\times n$ diagonal matrix $J$ whose first $n/2$ diagonal
elements have negative real parts and the rest have positive real parts. Next
we construct an $n\times n$ matrix $Z$ such that imaginary and real parts of
elements of $Z$ are uniformly and randomly taken from $(0,1)$. Then we take
$A=ZJZ^{-1}$. Let $\widehat{S}$ be the computed $\operatorname{sign}(A)$. We
measure the accuracy by relative forward error

$\operatorname{Err_{relative}}(\operatorname{sign}(A),\widehat{S})\vcentcolon=\frac{\lVert\operatorname{sign}(A)-\widehat{S}\rVert_{\max}}{\lVert\operatorname{sign}(A)\rVert_{\max}},$

where $\lVert\cdot\rVert_{\max}$ is defined in (47).

In Figure 3, we compare the efficiency and accuracy of Frobenius inversion and
inversion via LU decomposition. From the left plot in Figure 3, it is clear
that the Frobenius inversion is faster than inversion via LU decomposition.
The right plot in Figure 3 indicates that the Frobenius inversion is a bit
less accurate than inversion via LU decomposition. However, in practice, such
a small loss of accuracy can be ignored.

Figure 3. comparison on the matrix sign function

###  7.4. Sylvester equation

In this subsection, we apply Algorithms 1 and 3 to solve the Sylvester
equation:

$AX+XB=C,$ (48)

where $A\in\mathbb{C}^{m\times m},B\in\mathbb{C}^{n\times n},$ and
$C\in\mathbb{C}^{m\times n}$ are given. The goal is to solve for $X$. As noted
in [25, 51], given that $\operatorname{sign}(A)=I_{m}$ and
$\operatorname{sign}(B)=I_{n}$, we have

$\operatorname{sign}\left(\begin{bmatrix}A&-C\\\
0&-B\end{bmatrix}\right)=\begin{bmatrix}I&-2X\\\ 0&-I\end{bmatrix}.$

This implies that the Sylvester equation can be solved by computing the matrix
sign function. Thus the Sylvester equation can be solved by Newton iterations
as well. For simplicity, we denote by $X_{t}$ the variable in Newton
iterations:

$X_{t+1}=\frac{1}{2}(X_{t}+X_{t}^{-1}),\quad t\in\mathbb{N},$

with

$X_{0}=\begin{bmatrix}A&-C\\\ 0&-B\end{bmatrix}.$

It is noticeable that $X_{t}$ does not converge to $X$. Instead, $X_{t}$
converges to $\begin{bmatrix}I&-2X\\\ 0&-I\end{bmatrix}$. In each iteration,
we compute $X_{t}^{-1}$ either by Frobenius inversion or by inversion via LU
decomposition. This gives us two computed values of $X$ for each triple
$(A,B,C)$ of parameters in the Sylvester (48). We measure the progress in
Newton iteration by the relative change in $X_{t}$:

$\delta_{t}\vcentcolon=\frac{\lVert X_{t}-X_{t-1}\rVert_{\max}}{\lVert
X_{t}\rVert_{\max}},$

where $\lVert\cdot\rVert_{\max}$ is defined in (47). We stop the Newton
iteration when either $\delta_{t}\leq 10^{-1}$ or $t\geq 100$ is achieved. Let
$\widehat{X}$ be the computed $X$. We measure the accuracy by relative forward
error

$\operatorname{Err_{relative}}(X,\widehat{X})\vcentcolon=\frac{\lVert
X-\widehat{X}\rVert_{\max}}{\lVert X\rVert_{\max}},$

where $\lVert\cdot\rVert_{\max}$ is defined in (47).

We test and compare Algorithms 1 and 3 on the Sylvester (47) whose matrix
parameters $A,B,C$ have dimension $n$ between $1050$ and $2000$. We first
generate an $n\times n$ matrix $Z$ such that the real and imaginary parts of
its elements are taken uniformly and randomly taken from $(0,1)$. Next we
generate an $n\times n$ diagonal matrix $J$ such that the real and imaginary
parts of its diagonal elements are taken uniformly and randomly from the
interval $(9,10)$. Here we do not use the interval $(0,1)$ since it generates
matrices that are close to singular, which makes the algorithms inaccurate. We
take $A=ZJZ^{-1}$ and we generate $B$ in the same way. Lastly, we generate a
random complex $n\times n$ matrix $X$ such that the real and imaginary parts
of each element of $X$ are chosen uniformly and randomly from $(0,1)$ and take
$C=AX+XB$.

In Figure 4, we compare the efficiency and accuracy of Frobenius inversion and
inversion via LU decomposition on Sylvester equations whose parameters $A,B,C$
are randomly generated by the above procedure. One may easily see from Figure
4 that on Sylvester equations, the Frobenius inversion is faster than
inversion by LU decomposition and the loss of accuracy is ignorable.

Figure 4. comparison on the Sylvester equation

#### Lyapunov equation

As an important special case of the Sylvester equation, the Lyapunov equation
occurs in control and system theory [26]. To obtain Lyapunov equation, we
simply take $B=A^{*}$ in (48):

$AX+XA^{*}=C,$ (49)

where $A,C\in\mathbb{C}^{n\times n}$ are given. We apply Algorithms 1 and 3 to
solve (49) in the same way as what we do for Sylvester equations, except that
we take $B=A^{\ast}$ after we generate $A$.

From Figure 5, we see again that Frobenius inversion is faster than inversion
via LU decomposition with almost no loss of accuracy.

Figure 5. comparison on the Lyapunov equation

###  7.5. polar decomposition

In this subsection, we apply Algorithms 1 and 3 to compute the polar
decomposition. Given a complex matrix $A\in\mathbb{C}^{n\times n}$, its polar
decomposition is a matrix decomposition $A=UH$, where $U$ is unitary and $H$
is hermitian positive semidefinite. The Polar decomposition of $A$ can be
computed by the sign function of a block matrix [24, 25, 31]:

$\operatorname{sign}\left(\begin{bmatrix}0&A\\\
A^{*}&0\end{bmatrix}\right)=\begin{bmatrix}0&U\\\ U^{*}&0\end{bmatrix},$

where $U$ is the unitary factor of $A$ in its polar decomposition. Once we
obtain $U$, we can compute $H$ by $H=U^{\ast}A$. Therefore, we are able to
compute the polar decomposition of a matrix by Newton iterations [24, 25, 31]:

$X_{t+1}=\frac{1}{2}(X_{t}+X_{t}^{-*}),\quad t\in\mathbb{N},$

with $X_{0}=A$. Here $X_{t}^{-*}$ denotes the conjugate transpose of the
inverse of $X_{t}$. In each iteration, we compute $X_{t}^{-*}$ either by
Frobenius inversion or by inversion via LU decomposition. This gives us two
computed polar decompositions of $A$. We measure the progress in Newton
iterations by the relative change in $X_{t}$:

$\delta_{t}\vcentcolon=\frac{\lVert X_{t}-X_{t-1}\rVert_{\max}}{\lVert
X_{t}\rVert_{\max}},$

where $\lVert\cdot\rVert_{\max}$ is defined in (47). We stop the Newton
iterations when either $\delta_{t}\leq 10^{-3}$ or $t\geq 100$ is achieved.

We test Algorithms 1 and 3 on Polar decompositions of matrices whose dimension
$n$ is between $2100$ and $40000$. We generate a random complex $n\times n$
matrix $A$ as follows. We first generate two random complex matrices $B,C$
such that the real and imaginary parts of elements of $B,C$ are uniformly and
randomly taken from $(0,1)$. Then we compute the QR-decomposition of $B$ and
let $U$ be its unitary factor. Lastly, we compute $H=C^{*}C$ and let $A=UH$.

Let $\widehat{U}$ and $\widehat{H}$ be the factors in the computed polar
decomposition of $A$. We measure their accuracy by the relative forward error:

$\operatorname{Err_{relative}}(U,\widehat{U})\vcentcolon=\frac{\lVert
U-\widehat{U}\rVert_{\max}}{\lVert
U\rVert_{\max}},\quad\operatorname{Err_{relative}}(H,\widehat{H})\vcentcolon=\frac{\lVert
H-\widehat{H}\rVert_{\max}}{\lVert H\rVert_{\max}},$

where $\lVert\cdot\rVert_{\max}$ is defined in (47).

From the left plot in Figure 6, we see that Frobenius inversion is faster than
inversion via LU decomposition. Moreover, the right plot in Figure 6 indicates
that the two algorithms on Polar decompositions have the same accuracy.

Figure 6. comparison on the polar decomposition

###  7.6. hermitian positive matrices

In this section, we apply Algorithm 6, 8, 3, and 1 to compute the inverse of a
hermitian positive matrix $X$. Inverting hermitian positive matrices is an
important task in practice. For instance, in some MIMO radios, inversion of
some hermitian positive matrices are needed [58]. Thus, a faster algorithm for
inverting hermitian positive matrices leads to a more efficient procedure in
processing MIMO radios.

Let $3600\leq n\leq 6000$ be a positive integer. We first generate two real
matrices $A,B\in\mathbb{R}^{n\times n}$ where elements of $A$ and $B$ are
generated uniformly and randomly from $(0,1)$. Then we take
$X=(A+iB)(A+iB)^{*}+0.01I$. For each dimension $n$, we generate $10$ random
matrices by the above procedure.

For each randomly generated matrix $X$, we compute the inverse of $X$ by
Algorithm 6, 8, 3, and 1 respectively. For each fixed $n$, we average the
computation time for the $10$ random instances. In Figure 7, we exhibit the
computation time of the four algorithms versus the logarithmic dimension of
matrices. It is clear from Figure 7 that the variant of Frobenius inversion is
the fastest one.

Figure 7. comparison of efficiency for positive matrices

.

Next, we compare the accuracy of the four algorithms. For any hermitian
positive matrix $X$, let $\widehat{X}$ be the computed inverse of $X$. We
assess the accuracy of $\widehat{X}$ in terms of its left and right relative
residual defined in equation (46).

To avoid numerical issues, we only test the four algorithms on well-
conditioned matrices ( i.e., matrices whose condition number is $10$) whose
dimension goes from $2$ to $4096$. We generate a random matrix whose condition
number is $\kappa$ by the procedure that follows. We first generate $(n-2)$
integers from $1$ to $(\kappa-1)$ uniformly and randomly. Then, we consider
the $n$ by $n$ diagonal matrix $D$ whose diagonal elements are the $(n-2)$
integers we just generate together with $1$ and $\kappa$. Then, we generate a
random unitary matrix $U$ as follows. We first generate a random matrix $B$
such that the real and imaginary parts of elements of $B$ are uniformly and
randomly taken from $(0,1)$. Then we compute the QR-decomposition of $B$ and
let $U$ be its unitary factor. We next multiply $U$ to the left of $D$ and the
conjugate transpose of $U$ to the right of $D$ to obtain $UDU^{*}$. Since the
smallest and largest singular values of this matrix are $c$ and $\kappa c$ for
some $c\in\mathbb{R}$, the condition number of this matrix is $\kappa$. In our
experiment, we take $\kappa=10$. We note that there is no specific sparsity
patterns on the matrices we generate so that our experiment is not biased
toward a certain class of matrices with special sparsity patterns.

In Figure 8, we compare the left and right relative residuals of $\widehat{X}$
computed by the four algorithms. It is obvious from the left plot in Figure 8
that the variant of Frobenius inversion is slightly less accurate than the
other algorithms. However, the difference is so small that one can ignore it
for practical purposes.

Figure 8. comparison of relative residuals for positive matrices

## conclusion

In this paper, we analyze the Frobenius inversion. We show that Frobenius
inversion uses the least number of real matrix multiplications and real matrix
inversions among all complex matrix inversion algorithms. Then, we show that
it runs faster than the widely employed inversion algorithm based on the LU
decomposition, both theoretically and empirically. Next, we give three
applications of the Frobenius inversion: matrix sign function, Sylvester
equation, and polar decomposition. Finally, we give a variant of Frobenius
inversion on hermitian positive matrices. We show that this algorithm is
faster than the widely used inversion algorithm based on Cholesky
decomposition, both theoretically and empirically. This leads to a more
efficient procedure to process MIMO radios.

## References

  * [1] Z. Bai, J. Demmel, J. Dongarra, A. Petitet, H. Robinson, and K. Stanley.  The spectral decomposition of nonsymmetric matrices on distributed memory parallel computers.  SIAM J. Sci. Comput., 18(5):1446–1461, 1997.
  * [2] Z. Bai, J. Demmel, J. Dongarra, A. Petitet, H. Robinson, and K. Stanley.  The spectral decomposition of nonsymmetric matrices on distributed memory parallel computers.  SIAM J. Sci. Comput., 18(5):1446–1461, 1997.
  * [3] Z. Bai and J. W. Demmel.  Design of a parallel nonsymmetric eigenroutine toolbox.  University of Kentucky, Department of Mathematics, 1992.
  * [4] Y. Bard.  Nonlinear parameter estimation.  Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1974.
  * [5] A. N. Beavers, Jr. and E. D. Denman.  A computational method for eigenvalues and eigenvectors of a matrix with real eigenvalues.  Numer. Math., 21:389–396, 1973.
  * [6] A. S. Besicovitch.  On the linear independence of fractional powers of integers.  J. London Math. Soc., 15:3–6, 1940.
  * [7] E. Bodewig.  Matrix calculus.  North-Holland Publishing Co., Amsterdam; Interscience Publishers, Inc., New York, enlarged edition, 1959.
  * [8] P. Bürgisser, M. Clausen, and M. A. Shokrollahi.  Algebraic complexity theory, volume 315 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, 1997.  With the collaboration of Thomas Lickteig.
  * [9] R. Byers.  Solving the algebraic Riccati equation with the matrix sign function.  Linear Algebra Appl., 85:267–279, 1987.
  * [10] H. Cramér.  Mathematical Methods of Statistics.  Princeton Mathematical Series, vol. 9. Princeton University Press, Princeton, N. J., 1946.
  * [11] A. Ditkowski, G. Fibich, and N. Gavish.  Efficient solution of $Ax^{(k)}=b^{(k)}$ using $A^{-1}$.  J. Sci. Comput., 32(1):29–44, 2007.
  * [12] A. Druinsky and S. Toledo.  How accurate is $\operatorname{inv}({A})*b$?  arXiv preprint arXiv:1201.6035, 2012.
  * [13] K. Dudeck.  Solving complex systems using spreadsheets: A matrix decomposition approach.  In 2005 ASEE Annual Conference and Exposition, Conference Proceedings, pages 12875–12880, 2005.  2005 ASEE Annual Conference and Exposition: The Changing Landscape of Engineering and Technology Education in a Global World ; Conference date: 12-06-2005 Through 15-06-2005.
  * [14] S. Eberli, D. Cescato, and W. Fichtner.  Divide-and-conquer matrix inversion for linear mmse detection in sdr mimo receivers.  In 2008 NORCHIP, pages 162–167. IEEE, 2008.
  * [15] L. W. Ehrlich.  Complex matrix inversion versus real.  Comm. ACM, 13:561–562, 1970.
  * [16] M. El-Hawary.  Further comments on ”a note on the inversion of complex matrices”.  IEEE Transactions on Automatic Control, 20(2):279–280, 1975.
  * [17] D. K. Faddeev and V. N. Faddeeva.  Computational methods of linear algebra.  W. H. Freeman and Co., San Francisco-London, 1963.  Translated by Robert C. Williams.
  * [18] B. Fischer and J. Modersitzki.  Fast inversion of matrices arising in image processing.  Numer. Algorithms, 22(1):1–11, 1999.
  * [19] D. Y. Grigoriev.  Multiplicative complexity of a pair of bilinear forms and of the polynomial multiplication.  In J. Winkowski, editor, Mathematical Foundations of Computer Science 1978, pages 250–256, Berlin, Heidelberg, 1978. Springer Berlin Heidelberg.
  * [20] L. Guttman.  Enlargement methods for computing the inverse matrix.  Ann. Math. Statistics, 17:336–343, 1946.
  * [21] W. W. Hager.  Updating the inverse of a matrix.  SIAM Rev., 31(2):221–239, 1989.
  * [22] M. T. Heath, G. A. Geist, and J. B. Drake.  Early experience with the intel ipsc/860 at oak ridge national laboratory.  The International Journal of Supercomputing Applications, 5(2):10–26, 1991.
  * [23] H. V. Henderson and S. R. Searle.  On deriving the inverse of a sum of matrices.  SIAM Rev., 23(1):53–60, 1981.
  * [24] N. J. Higham.  Computing the polar decomposition—with applications.  SIAM J. Sci. Statist. Comput., 7(4):1160–1174, 1986.
  * [25] N. J. Higham.  The matrix sign decomposition and its relation to the polar decomposition.  In Proceedings of the 3rd ILAS Conference (Pensacola, FL, 1993), volume 212/213, pages 3–20, 1994.
  * [26] N. J. Higham.  Accuracy and stability of numerical algorithms.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2002.
  * [27] N. J. Higham.  Functions of matrices.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008.  Theory and computation.
  * [28] N. J. Higham and P. Papadimitriou.  A parallel algorithm for computing the polar decomposition.  Parallel Comput., 20(8):1161–1173, 1994.
  * [29] J. L. Howland.  The sign matrix and the separation of matrix eigenvalues.  Linear Algebra Appl., 49:221–232, 1983.
  * [30] A. A. Karatsuba.  The complexity of computations.  Trudy Mat. Inst. Steklov., 211(Optim. Upr. i Differ. Uravn.):186–202, 1995.
  * [31] C. Kenney and A. J. Laub.  On scaling Newton’s method for polar decomposition and the matrix sign function.  SIAM J. Matrix Anal. Appl., 13(3):698–706, 1992.
  * [32] J. Kim and F. L. Teixeira.  Parallel and explicit finite-element time-domain method for Maxwell’s equations.  IEEE Trans. Antennas and Propagation, 59(6, part 2):2350–2356, 2011\.
  * [33] A. Klein and G. Mélard.  Computation of the Fisher information matrix for time series models.  J. Comput. Appl. Math., 64(1-2):57–68, 1995.
  * [34] D. E. Knuth.  The art of computer programming. Vol. 2.  Addison-Wesley, Reading, MA, 1998.  Seminumerical algorithms, Third edition.
  * [35] A. Kochnev and N. Savelov.  Symmetric matrix inversion using modified gaussian elimination.  arXiv preprint arXiv:1504.06734, 2015.
  * [36] C. Lanczos.  Applied analysis.  Prentice-Hall, Inc., Englewood Cliffs, N. J., 1956.
  * [37] J. M. Landsberg.  Tensors: geometry and applications, volume 128 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2012.
  * [38] T. Lickteig.  The computational complexity of division in quadratic extension fields.  SIAM J. Comput., 16(2):278–311, 1987.
  * [39] C.-c. Lin and E. Zmijewski.  A parallel algorithm for computing the eigenvalues of an unsymmetric matrix on an SIMD mesh of processors.  University of California, Santa Barbara, College of Engineering, Department …, 1991.
  * [40] K. Lo.  Several numerical methods for matrix inversion.  International Journal of Electrical Engineering Education, 15(2):131–141, 1978.
  * [41] J. H. Maindonald.  Statistical computation.  Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. John Wiley & Sons, Inc., New York, 1984.
  * [42] P. McCullagh.  Generalized linear models.  European Journal of Operational Research, 16(3):285–292, 1984.
  * [43] P. McCullagh and J. A. Nelder.  Generalized linear models.  Monographs on Statistics and Applied Probability. Chapman & Hall, London, 1989.  Second edition [of MR0727836].
  * [44] P. Mukherjee and L. Satish.  On the inverse of forward adjacency matrix.  arXiv preprint arXiv:1711.09216, 2017.
  * [45] I. Munro.  Some results concerning efficient and optimal algorithms.  In Proceedings of the Third Annual ACM Symposium on Theory of Computing, STOC ’71, page 40–44, New York, NY, USA, 1971. Association for Computing Machinery.
  * [46] S. K. Panda.  Inverses of bicyclic graphs.  Electron. J. Linear Algebra, 32:217–231, 2017.
  * [47] S. Pavlíková.  A note on inverses of labeled graphs.  Australas. J. Combin., 67:222–234, 2017.
  * [48] G. B. Price.  An introduction to multicomplex spaces and functions, volume 140 of Monographs and Textbooks in Pure and Applied Mathematics.  Marcel Dekker, Inc., New York, 1991.  With a foreword by Olga Taussky Todd.
  * [49] C. Radhakrishna Rao.  Information and the accuracy attainable in the estimation of statistical parameters.  Bull. Calcutta Math. Soc., 37:81–91, 1945.
  * [50] K. S. Riedel.  A Sherman-Morrison-Woodbury identity for rank augmenting matrices with application to centering.  SIAM J. Matrix Anal. Appl., 13(2):659–662, 1992.
  * [51] J. D. Roberts.  Linear model reduction and solution of the algebraic Riccati equation by use of the sign function.  Internat. J. Control, 32(4):677–687, 1980.
  * [52] J. Schur.  Über potenzreihen, die im innern des einheitskreises beschränkt sind.  Journal für die reine und angewandte Mathematik (Crelles Journal), 1918:122 – 145, 1918.
  * [53] C. Segre.  Le rappresentazioni reali delle forme complesse e gli enti iperalgebrici.  Math. Ann., 40(3):413–467, 1892.
  * [54] W. W. Smith and J. Smith.  Handbook of real-time fast Fourier transforms.  IEEE New York, 1995.
  * [55] W. W. Smith, Jr. and S. Erdman.  A note on the inversion of complex matrices.  IEEE Trans. Automatic Contol, AC-19:64, 1974.
  * [56] G. W. Stewart.  Matrix algorithms. Vol. I.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 1998\.  Basic decompositions.
  * [57] V. Strassen.  Gaussian elimination is not optimal.  Numer. Math., 13:354–356, 1969.
  * [58] C. Studer, S. Fateh, and D. Seethaler.  Asic implementation of soft-input soft-output mimo detection using mmse parallel interference cancellation.  IEEE Journal of Solid-State Circuits, 46(7):1754–1765, 2011.
  * [59] L. Tornheim.  Inversion of a complex matrix.  Comm. ACM, 4:398, 1961.
  * [60] S. Winograd.  On the number of multiplications necessary to compute certain functions.  Comm. Pure Appl. Math., 23:165–179, 1970.
  * [61] M. A. Woodbury.  The Stability of Out-Input Matrices.  Chicago, Ill., 1949.
  * [62] M. A. Woodbury.  Inverting modified matrices.  Princeton University, Princeton, N. J., 1950.  Statistical Research Group, Memo. Rep. no. 42,.
  * [63] D. Ye, Y. Yang, B. Mandal, and D. J. Klein.  Graph invertibility and median eigenvalues.  Linear Algebra Appl., 513:304–323, 2017.
  * [64] A. Zielinski.  On inversion of complex matrices.  Internat. J. Numer. Methods Engrg., 14(10):1563–1566, 1979.