# Numerical stability and tensor nuclear norm

Zhen Dai  and  Lek-Heng Lim  Computational and Applied Mathematics Initiative,
Department of Statistics, University of Chicago, Chicago, IL 60637-1514.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We present a notion of bilinear stability, which is to numerical stability
what bilinear complexity is to time complexity. In bilinear complexity, an
algorithm for evaluating a bilinear operator $\beta\mathrel{\mathop{\mathchar
58\relax}}\mathbb{U}\times\mathbb{V}\to\mathbb{W}$ is a decomposition
$\beta=\varphi_{1}\otimes\psi_{1}\otimes
w_{1}+\dots+\varphi_{r}\otimes\psi_{r}\otimes w_{r}$; the number of terms $r$
captures the speed of the algorithm; and its smallest possible value, i.e.,
the tensor rank of $\beta$, quantifies the speed of a fastest algorithm.
Bilinear stability introduces norms to the mix: The growth factor of the
algorithm $\lVert\varphi_{1}\rVert_{*}\lVert\psi_{1}\rVert_{*}\lVert
w_{1}\rVert+\dots+\lVert\varphi_{r}\rVert_{*}\lVert\psi_{r}\rVert_{*}\lVert
w_{r}\rVert$ captures the accuracy of the algorithm; and its smallest possible
value, i.e., the tensor nuclear norm of $\beta$, quantifies the accuracy of a
stablest algorithm. To substantiate this notion, we establish a bound for the
forward error in terms of the growth factor and present numerical evidence
comparing various fast algorithms for matrix and complex multiplications,
showing that larger growth factors correlate with less accurate results.
Compared to similar studies of numerical stability, bilinear stability is more
general, applying to any bilinear operators and not just matrix or complex
multiplications; is more simplistic, bounding forward error in terms of a
single (growth) factor; and is truly tensorial like bilinear complexity,
invariant under any orthogonal change of coordinates. As an aside, we study a
new algorithm for computing complex multiplication in terms of real, much like
Gauss’s, but is optimally fast and stable in that it attains both tensor rank
and nuclear norm.

##  1\. Introduction

More than fifty years ago, Volker Strassen announced an astounding result: A
pair of $2\times 2$ matrices may be multiplied with seven multiplications
[30]. A consequence is that linear systems can be solved in $O(n^{\log_{2}7})$
time complexity, a surprise at that time as existing works such as [21]
purportedly showed that $O(n^{3})$ was the lowest possible.

Strassen’s algorithm is in the spirit of the well-known algorithm, often
attributed to Gauss,111See [26, p. 37], [27, p. 8] for example. for
multiplying a pair of complex numbers with three real multiplications [15],

$(a+bi)(c+di)=(ac-bd)+i[(a+b)(c+d)-ac-bd],$ (1.1)

but is notable in that Strassen’s applies to a noncommutative product (matrix
multiplication) as opposed to a commutative one (complex scalar
multiplication). It led to a plethora of followed-up works and ultimately to
the realization that there is a unified framework underlying the algorithms of
Gauss and Strassen, namely, in evaluating a bilinear operator
$\beta\colon\mathbb{U}\times\mathbb{V}\to\mathbb{W}$, viewed as a $3$-tensor
in $\mathbb{U}^{*}\otimes\mathbb{V}^{*}\otimes\mathbb{W}$, any decomposition

$\beta=\varphi_{1}\otimes\psi_{1}\otimes
w_{1}+\dots+\varphi_{r}\otimes\psi_{r}\otimes w_{r}$ (1.2)

into linear functionals $\varphi_{i}\colon\mathbb{U}\to\mathbb{R}$,
$\psi_{i}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{V}\to\mathbb{R}$, and
vectors $w_{i}\in\mathbb{W}$, $i=1,\dots,r$, gives us an algorithm for
computing $\beta$. Furthermore, the number of terms $r$ in such a
decomposition counts precisely the number of multiplications, and thus the
minimal value of $r$, i.e., the _tensor rank_ of $\beta$, gives the optimal
complexity for evaluating $\beta$ in an appropriate sense [31] (see Section
2). Both Gauss’s and Strassen’s algorithms are the fastest possible according
to this measure, that is, they attain the tensor ranks of complex
multiplication (three) and $2\times 2$ matrix product (seven) respectively
[35].

As is well-known to practitioners, speed is not all that matters in an
algorithm, numerical stability is arguably more important in finite-precision
computations as rounding errors may result in an unstable algorithm producing
no correct digits. While the stability of algorithms for evaluating bilinear
operators has been studied for specific algorithms or operators in isolation,
e.g., for Gauss’s algorithm in [15], Strassen’s algorithm in [7], and other
fast matrix multiplication algorithms in [2, 4], there has been no unfied
treatment that applies to all bilinear operators $\beta$ as in the case of
speed. There is no analysis that quantifies stability in terms of some
tensorial property of $\beta$ analogous to how speed is quantified in terms of
its tensor rank. The goal of the present article is to fill this gap. We will
show that just as the number of terms $r$ in the decomposition (1.2) controls
the speed of the algorithm, the _growth factor_ , defined as

$\lVert\varphi_{1}\rVert_{*}\lVert\psi_{1}\rVert_{*}\lVert
w_{1}\rVert+\dots+\lVert\varphi_{r}\rVert_{*}\lVert\psi_{r}\rVert_{*}\lVert
w_{r}\rVert,$ (1.3)

controls the stability of the algorithm; and just as the tensor rank of
$\beta$ measures the optimal speed, the _tensor nuclear norm_ of $\beta$,
defined as

$\lVert\beta\rVert_{\nu}\coloneqq\inf\biggl{\\{}\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert\mathrel{\mathop{\mathchar
58\relax}}\beta=\sum_{i=1}^{r}\varphi_{i}\otimes\psi_{i}\otimes
w_{i}\biggr{\\}},$ (1.4)

measures the optimal stability, the precise meaning of which we will state in
due course.

Although we have alluded to the relation between tensor nuclear norm and
numerical stability in earlier works [25, 36, 13], we have never stated a
precise relation nor carried out numerical experiments to demonstrate the
relation. This article provides both. Theorem 3.3 gives a general relation
between the growth factor of a bilinear algorithm and its forward error, from
which a relation between tensor nuclear norm and forward error may be deduced
as in Corollary 3.4. We then perform a range of numerical experiments
involving Gauss’s and Strassen’s algorithms to substantiate our theoretical
findings:

Matrix multiplication:

We compare Strassen’s algorithm with a well-known variant due to Winograd [17,
22]. While both attain the optimal seven multiplications, Winograd’s variant
is often favored because it requires only fifteen additions, compared to
Strassen’s eighteen. Nevertheless we will show that Strassen’s algorithm has a
growth factor of $12+2\sqrt{2}\approx 14.83$ whereas Winograd’s variant has a
growth factor of $7+4\sqrt{2}+3\sqrt{3}\approx 17.85$. For comparison, the
conventional algorithm for $2\times 2$ matrix product has eight
multiplications and a growth factor of $8$. Our numerical experiments confirm
that in terms of accuracy Winograd’s is indeed worse than Strassen’s, which is
in turn worse than the conventional algorithm, as Theorem 3.3 indicates.

Complex multiplication:

We compare the regular algorithm for complex multiplication, which requires
four real multiplications and has a growth factor of $4$; Gauss’s algorithm,
which requires three real multiplications but has a larger growth factor of
$2(1+\sqrt{2})\approx 4.83$; and a new algorithm:

$\displaystyle(a+bi)(c+di)$
$\displaystyle=\frac{1}{2}\biggl{[}\biggl{(}a+\frac{1}{\sqrt{3}}b\biggr{)}\biggl{(}c+\frac{1}{\sqrt{3}}d\biggr{)}+\biggl{(}a-\frac{1}{\sqrt{3}}b\biggr{)}\biggl{(}c-\frac{1}{\sqrt{3}}d\biggr{)}-\frac{8}{3}bd\biggr{]}$
(1.5)
$\displaystyle\qquad+\frac{i\sqrt{3}}{2}\biggl{[}\biggl{(}a+\frac{1}{\sqrt{3}}b\biggr{)}\biggl{(}c+\frac{1}{\sqrt{3}}d\biggr{)}-\biggl{(}a-\frac{1}{\sqrt{3}}b\biggr{)}\biggl{(}c-\frac{1}{\sqrt{3}}d\biggr{)}\biggr{]}.$

This new algorithm has the best features of both the regular and Gauss’s
algorithms, requiring three real multiplications and yet has the smaller (in
fact, smallest, as we will see) growth factor of $4$. Again the results are
consistent with the prediction of Theorem 3.3.

For the uniniated, we would like to stress that the aforementioned algorithms
only begin to make a difference when they are applied recursively, or applied
to matrices, or both. For instance, Gauss’s algorithm (1.1) is really quite
useless for multiplying a pair of complex numbers, whether ‘by hand’ or on a
computer. It only becomes useful when applied recursively in the form of
Karatsuba’s algorithm [19] for integer multiplication, with $i$ replaced by
the number base; or when applied to complex matrices [12]:

$(A+iB)(C+iD)=(AC-BD)+i[(A+B)(C+D)-AC-BD],$ (1.6)

with $A+iB,C+iD\in\mathbb{C}^{n\times n}$, $A,B,C,D\in\mathbb{R}^{n\times n}$.
As multiplication of matrices is much more expensive than addition of
matrices, so (1.6) really does represent an enormous savings in speed over the
regular algorithm:

$(A+iB)(C+iD)=(AC-BD)+i(BC+AD).$ (1.7)

Likewise, our new algorithm (1.5) only begins to make a difference when the
quantities involved are matrices. For the same reasons, the algorithms of
Strassen and Winograd are only worth the trouble when applied recursively to a
product of $n\times n$ matrices partitioned recursively into $2\times 2$
blocks.

To address another related point early on, a surprisingly common complaint
among early feedbacks is that there are a lot of $\sqrt{3}$’s in our algorithm
(1.5). Certainly, if one computes these products ‘by hand,’ it would be easier
to use the regular or Gauss’s algorithm since they do not involve irrational
coefficients. But when performed by a computer this is completely immaterial.
In case it is not clear, it does not matter whether we multiply by $3$ or by
$\sqrt{3}$; to a computer (or any IEEE 754-compliant equipment) both are
binary strings of $0$’s and $1$’s and arithmetic takes one flop regardless.
Maybe there would be some minor savings when a constant happens to be a power
of $2$, because of binary arithmetic; but aside from that, it makes no
difference what constants we have in our algorithm.

For the matrix multiplication experiments, our goal is simply to illustrate
Theorem 3.3 by comparing the known algorithms of Strassen and Winograd. But
for the complex multiplication experiments, we also have the additional goal
of testing, for the first time, the new algorithm (1.5) applied to multiply
complex matrices, which we will see is

  * •

nearly as fast as Gauss’s algorithm (1.6), and

  * •

nearly as accurate as the regular algorithm (1.7).

To substantiate these claims, we perform more extensive experiments to compare
(1.5), (1.6), and (1.7), including three practical applications: evaluation of
matrix polynomials via Horner’s method [18], unitary transform, and complex-
valued neural networks [1, 3, 9, 29, 34, 37]. All our codes are available from
https://github.com/zhen06/Complex-Matrix-Multiplication.

### Notations

To reduce notational clutter, we denote norms on different vector spaces
$\mathbb{U},\mathbb{V},\mathbb{W}$ by the same $\lVert\,\cdot\,\rVert$. There
is no cause for confusion since we always use it in a form like $\lVert
v\rVert$ for some $v\in\mathbb{V}$, where it is clear from context that
$\lVert\,\cdot\,\rVert$ refers to a norm on $\mathbb{V}$. Likewise the
corresponding dual norms on $\mathbb{U}^{*},\mathbb{V}^{*},\mathbb{W}^{*}$
will be denoted by the same $\lVert\,\cdot\,\rVert_{*}$. Recall that for
$\varphi\in\mathbb{V}^{*}$, i.e., $\varphi\mathrel{\mathop{\mathchar
58\relax}}\mathbb{V}\to\mathbb{R}$ is a linear functional, this is defined by

$\lVert\varphi\rVert_{*}\coloneqq\sup\\{\lvert\varphi(v)\rvert\mathrel{\mathop{\mathchar
58\relax}}\lVert v\rVert\leq 1\\}.$

##  2\. Bilinear complexity

We provide a brief review of bilinear complexity, usually studied in Algebraic
Computational Complexity [6, 8, 24, 33], for numerical analysts. Our goals
here are to (i) highlight certain departures from typical practice in
numerical linear algebra; and (ii) show a parallel with our notion of bilinear
stability in the next section.

Let $\mathbb{U}$, $\mathbb{V}$, $\mathbb{W}$ be finite-dimensional vector
spaces, assume to be over $\mathbb{R}$ for simplicity. Let
$\beta\mathrel{\mathop{\mathchar
58\relax}}\mathbb{U}\times\mathbb{V}\to\mathbb{W}$ be a bilinear operator.
Depending on one’s definition of a tensor, we have
$\beta\in\mathbb{U}^{*}\otimes\mathbb{V}^{*}\otimes\mathbb{W}$ either through
definition [25, Definition 3.3] or by the universal mapping property [25,
Equation 4.88]. A _bilinear algorithm_ for evaluating $\beta$ is a
decomposition of the form (1.2). In other words, for any $u\in\mathbb{U}$ and
$v\in\mathbb{V}$, we evaluate $\beta(u,v)$ by performing the algorithm given
by the decomposition on the right:

$\beta(u,v)=\sum_{i=1}^{r}\varphi_{i}(u)\psi_{i}(v)w_{i}.$ (2.1)

In practice, the vector spaces involved are usually Euclidean spaces of
vectors $\mathbb{R}^{n}$ or matrices $\mathbb{R}^{m\times n}$. Riesz
representation theorem guarantees that any linear functional
$\varphi\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\to\mathbb{R}$ must
take the form $\varphi(x)=a^{\scriptscriptstyle\mathsf{T}}x$ for some
$a\in\mathbb{R}^{n}$ and likewise any functional
$\varphi\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{m\times
n}\to\mathbb{R}$ must take the form
$\varphi(X)=\operatorname{tr}(A^{\scriptscriptstyle\mathsf{T}}X)$ for some
$A\in\mathbb{R}^{m\times n}$.

Each rank-one term $\varphi_{i}(u)\psi_{i}(v)w_{i}$ in (2.1) accounts for one
multiplication but herein lies a pitfall — the ‘multiplication’ refers to the
product of $\varphi_{i}(u)$ and $\psi_{i}(v)$; note that this a variable
product, i.e., the value depends on variables $u$ and $v$, as opposed to a
scalar product. Take a randomly made-up example222Genuine examples to follow
in Sections 4 and 5. with $\mathbb{U}=\mathbb{R}^{2\times 2}$,
$\mathbb{V}=\mathbb{R}^{2}$, $\mathbb{W}=\mathbb{R}^{3}$, and

$\varphi_{i}\bigl{(}\begin{bmatrix}a&b\\\
c&d\end{bmatrix}\bigr{)}=\operatorname{tr}\bigl{(}\begin{bmatrix}-1&0\\\
1&2\end{bmatrix}^{\scriptscriptstyle\mathsf{T}}\begin{bmatrix}a&b\\\
c&d\end{bmatrix}\bigr{)}=-a+c+2d,\quad\psi_{i}\bigl{(}\begin{bmatrix}x\\\
y\end{bmatrix}\bigr{)}=\bigl{[}\begin{smallmatrix}3\\\
-1/2\end{smallmatrix}\bigr{]}^{\scriptscriptstyle\mathsf{T}}\begin{bmatrix}x\\\
y\end{bmatrix}=3x-y/2,\quad w_{i}=\begin{bmatrix}-3\\\ 4\\\
\sqrt{5}\end{bmatrix},$

then there is exactly one multiplication in

$\varphi_{i}(u)\psi_{i}(v)w_{i}=\begin{bmatrix}-3(-a+c+2d)(3x-y/2)\\\
4(-a+c+2d)(3x-y/2)\\\ \sqrt{5}(-a+c+2d)(3x-y/2)\end{bmatrix}.$

The scalar products like $2d$ or $-y/2$ or $\sqrt{5}t$ are discounted in
Strassen’s model of bilinear complexity [31, 32] and for good reasons — these
constants coefficients are fixed in the algorithm and can be hardcoded or
hardwired, unlike the product between $-a+c+2d$ and $3x-y/2$, which depends on
the variable inputs $u=\begin{bmatrix}a&b\\\ c&d\end{bmatrix}$ and
$v=\begin{bmatrix}x\\\ y\end{bmatrix}$. In particular, Strassen’s measure of
speed, called _bilinear complexity_ , is independent of the values of these
constant coefficients, but we will show in the next section that these will
affect numerical stability of the algorithm.

To emphasize its distinction from scalar products, Strassen calls a variable
product in the above sense a _nonscalar product_ [32]. In other words,
bilinear complexity measures speed purely in terms of the number of nonscalar
products. The bilinear complexity of the algorithm in (2.1) is given by the
number terms in the decomposition $r$ and the optimal speed of evaluating
$\beta$ is therefore given by the tensor rank [31]

$\operatorname{rank}(\beta)\coloneqq\min\biggl{\\{}r\mathrel{\mathop{\mathchar
58\relax}}\beta=\sum_{i=1}^{r}\varphi_{i}\otimes\psi_{i}\otimes
w_{i}\biggr{\\}}.$ (2.2)

A _tensor rank decomposition_ of $\beta$, i.e., one that attains its tensor
rank, is then a fastest algorithm in the context of bilinear complexity.

In realistic scenarios, storage and computations both have finite-precision.
Given $u$ and $v$, we do not need to know $\beta(u,v)$ exactly; in fact
computing anything beyond $16$ decimal digits of accuracy is wasted effort
since we do not store more than $16$ digits in IEEE double precision. So the
tensor rank of $\beta$ is less relevant than the _border rank_ [5] of $\beta$,
which is the smallest $r$ so that

$\lVert\beta-\varphi_{1}^{\varepsilon}\otimes\psi_{1}^{\varepsilon}\otimes
w_{1}^{\varepsilon}-\varphi_{2}^{\varepsilon}\otimes\psi_{2}^{\varepsilon}\otimes
w_{2}^{\varepsilon}-\dots-\varphi_{r}^{\varepsilon}\otimes\psi_{r}^{\varepsilon}\otimes
w_{r}^{\varepsilon}\rVert<\varepsilon$

for all $\varepsilon>0$, or, formally,

$\operatorname{\overline{\operatorname{rank}}}(\beta)\coloneqq\min\biggl{\\{}r\mathrel{\mathop{\mathchar
58\relax}}\beta=\lim_{\varepsilon\to
0^{+}}\sum_{i=1}^{r}\varphi_{i}^{\varepsilon}\otimes\psi_{i}^{\varepsilon}\otimes
w_{i}^{\varepsilon}\biggr{\\}}.$ (2.3)

For the two problems studied in our article, namely, matrix multiplication,

$\beta_{m,n,p}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{m\times
n}\times\mathbb{R}^{n\times p}\to\mathbb{R}^{m\times p},\quad(A,B)\mapsto AB,$

and complex multiplication,

$\beta_{\mathbb{C}}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{C}\times\mathbb{C}\to\mathbb{C},\quad(w,z)\mapsto wz,$

(noting that $\mathbb{C}$ is a two-dimensional real vector space), we have
[23, 35]

$\operatorname{rank}(\beta_{2,2,2})=\operatorname{\overline{\operatorname{rank}}}(\beta_{2,2,2})=7,\qquad\operatorname{rank}(\beta_{\mathbb{C}})=\operatorname{\overline{\operatorname{rank}}}(\beta_{\mathbb{C}})=3.$

It is in general difficult to find such exact values. For instance, the values
of $\operatorname{rank}(\beta_{3,3,3})$ and
$\operatorname{\overline{\operatorname{rank}}}(\beta_{3,3,3})$ are still
unknown. Most of the efforts in studying matrix multiplication go towards
determining the asymptotic value
$\omega\coloneqq\inf\\{p\in\mathbb{R}\colon\operatorname{rank}(\beta_{n,n,n})=O(n^{p})\\}$,
called the _exponent of matrix multiplication_. An advantage is that
asymptotically, the full arithmetic complexity, i.e., counting all operations
and not just nonscalar multiplications, is also $O(n^{\omega})$. More
importantly, the role of $\omega$ stretches far beyond matrix multiplication,
governing the full arithmetic complexity of computing inverse, determinant,
null basis, linear systems, LU/QR/eigenvalue/Hessenberg decompositions,
characteristic polynomials, sparsification, and even linear programming — note
in particular that none of these are bilinear operations [32] (see also [8,
Chapter 16] and [25, Examples 3.10 and 4.40].

##  3\. Bilinear stability

We would like to state at the outset that numerical stability is a moderately
complicated issue that depends on many factors and cannot be completely
represented by any single number. Designing numerically stable algorithms is
as much an art as it is a science. However the six _Higham guidelines_ for
numerical stability [17, Section 1.18] capture the most salient aspects. Among
them, the second guideline to “minimize the size of intermediate quantities
relative to the final solution” is one of the most unequivocal, lends itself
to precise quantification, and is what we will focus on in this section.
Consideration of Higham’s second guideline for bilinear algorithms leads us
naturally to the notion of _bilinear stability_ , which relates to accuracy
the way bilinear complexity relates to speed. More precisely, the growth
factor (1.3) and tensor nuclear norm (1.4) are to accuracy in bilinear
stability what the number of rank-$1$ terms in (2.1) and the tensor rank (2.2)
are to speed in bilinear complexity. Here accuracy refers to the size of
relative forward error.

Bilinear stability differs from existing studies of numerical stability of
bilinear algorithms such as those in [2, 4, 7, 15] in three ways: (i) it is
more general, applying to any bilinear operators as opposed to specific ones
like matrix multiplication; (ii) it is more simplistic, relating forward error
to just growth factor as opposed to two or three different factors in the
approaches of [2, 4]; (iii) it is truly tensorial, as growth factor and tensor
nuclear norm are invariant under any orthogonal change-of-coordinates, just as
tensor rank is invariant under any invertible change-of-coordinates. The
factors (i) and (ii), i.e., generality and simplicity, may often be sacrificed
for better bounds: Given any specific bilinear operator, we may often obtain
smaller forward error bounds by performing a more precise analysis tailored to
that given operator. We will do see this in Section 5.2.

One difference between bilinear complexity and bilinear stability is that the
latter requires a norm. While there are many excellent treatises on tensor
norms [10, 11, 28], they are excessive for our purpose. All the reader needs
to know is that for a vector space $\mathbb{V}_{i}$ with norm
$\lVert\,\cdot\,\rVert_{i}$, $i=1,\dots,d$, a tensor norm
$\lVert\,\cdot\,\rVert$ on
$\mathbb{V}_{1}\otimes\mathbb{V}_{2}\otimes\dots\otimes\mathbb{V}_{d}$
satisfies the multiplicativity property for rank-$1$ tensors:

$\lVert v_{1}\otimes v_{2}\otimes\dots\otimes v_{d}\rVert=\lVert
v_{1}\rVert_{1}\lVert v_{2}\rVert_{2}\cdots\lVert v_{d}\rVert_{d},$

where $v_{i}\in\mathbb{V}_{i}$. In particular, the spectral, Frobenius (also
called Hilbert–Schmidt), nuclear norms [25, p. 561 and Example 4.17] are all
equal on rank-$1$ tensors in
$\mathbb{U}^{*}\otimes\mathbb{V}^{*}\otimes\mathbb{W}$, i.e.,

$\lVert\varphi\otimes\psi\otimes
w\rVert_{\sigma}=\lVert\varphi\otimes\psi\otimes
w\rVert_{\scriptscriptstyle\mathsf{F}}=\lVert\varphi\otimes\psi\otimes
w\rVert_{\nu}=\lVert\varphi\rVert_{*}\lVert\psi\rVert_{*}\lVert w\rVert$

for all $\varphi\in\mathbb{U}^{*}$, $\psi_{i}\in\mathbb{V}^{*}$,
$w\in\mathbb{W}$. Consequently, when we speak of the norm of a rank-$1$ tensor
$\varphi\otimes\psi\otimes w$, it does not matter which of these three norms
we choose, and we will simply write

$\lVert\varphi\otimes\psi\otimes
w\rVert\coloneqq\lVert\varphi\rVert_{*}\lVert\psi\rVert_{*}\lVert w\rVert.$

We first present a straightforward heurstic that motivates our definition of
the growth factor, deferring the more formal forward error analysis to Theorem
3.3. If we apply the rank-one bilinear operator
$\varphi_{i}\otimes\psi_{i}\otimes w_{i}$ to $u$ and $v$,

$\displaystyle\lVert(\varphi_{i}\otimes\psi_{i}\otimes
w_{i})(u,v)\rVert=\lVert\varphi_{i}(u)\psi_{i}(v)w_{i}\rVert$
$\displaystyle=\lvert\varphi_{i}(u)\rvert\lvert\psi_{i}(v)\rvert\lVert
w_{i}\rVert$ $\displaystyle\leq\lVert\varphi_{i}\rVert_{*}\lVert
u\rVert\lVert\psi_{i}\rVert_{*}\lVert v\rVert\lVert
w_{i}\rVert=\lVert\varphi_{i}\otimes\psi_{i}\otimes w_{i}\rVert\lVert
u\rVert\lVert v\rVert.$

So $\varphi_{i}\otimes\psi_{i}\otimes w_{i}$ magnifies the errors in $u$ and
$v$ by an amount bounded by its tensor norm
$\lVert\varphi_{i}\otimes\psi_{i}\otimes w_{i}\rVert$. Therefore, in a
bilinear algorithm given by the right side of (2.1) for evaluating $\beta$,
triangle inequality gives

$\lVert\beta(u,v)\rVert=\biggl{\lVert}\sum_{i=1}^{r}(\varphi_{i}\otimes\psi_{i}\otimes
w_{i})(u,v)\biggr{\rVert}\leq\biggl{[}\sum_{i=1}^{r}\lVert\varphi_{i}\otimes\psi_{i}\otimes
w_{i}\rVert\biggr{]}\lVert u\rVert\lVert v\rVert.$

The $i$th step of the algorithm magnifies the error in the inputs $(u,v)$ by
an amount bounded by $\lVert\varphi_{i}\otimes\psi_{i}\otimes w_{i}\rVert$ and
over the course of $r$ steps in the algorithm, the accumulated error is
bounded by a factor of

$\sum_{i=1}^{r}\lVert\varphi_{i}\otimes\psi_{i}\otimes
w_{i}\rVert=\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert,$ (3.1)

which we will define as the _growth factor_ of the algorithm or decomposition
(2.1). Its minimum value over all possible bilinear algorithms for evaluating
$\beta$ or, equivalently, over all decomposition of $\beta$ as a $3$-tensor is
therefore given by the nuclear norm (1.4). This idea was first floated in [36,
Section 3.2]. Note that the growth factor depends on the
algorithm/decomposition for $\beta$ but the nulcear norm depends only on
$\beta$.

We now state a formal definition to make precise the terms used in the
preceding discussions.

######  Definition 3.1.

Let $\mathbb{U},\mathbb{V},\mathbb{W}$ be three finite-dimensional real vector
spaces. A _decomposition_ of a bilinear operator
$\beta\mathrel{\mathop{\mathchar
58\relax}}\mathbb{U}\times\mathbb{V}\to\mathbb{W}$ is a sequence
$\mathcal{D}=(\varphi_{i},\psi_{i},w_{i})_{i=1}^{r}$ with

$\beta=\sum_{i=1}^{r}\varphi_{i}\otimes\psi_{i}\otimes w_{i},$ (3.2)

where $\varphi_{i}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{U}\to\mathbb{R}$ and $\psi_{i}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{V}\to\mathbb{R}$ are linear functionals and
$w_{i}\in\mathbb{W}$, $i=1,\dots,r$. An _algorithm_
$\widehat{\beta}_{\mathcal{D}}$ given by the decomposition $\mathcal{D}$ takes
$(u,v)\in\mathbb{U}\times\mathbb{V}$ as inputs and computes the output
$\beta(u,v)$ in three steps:

  1. (i)

computes $\varphi_{i}(u)$ and $\psi_{i}(v)$, $i=1,\dots,r$;

  2. (ii)

computes $\varphi_{i}(u)\psi_{i}(v)w_{i}$, $i=1,\dots,r$;

  3. (iii)

computes $\sum_{i=1}^{r}\varphi_{i}(u)\psi_{i}(v)w_{i}$.

The _growth factor_ of the algorithm $\widehat{\beta}_{\mathcal{D}}$ is
defined as

$\gamma(\widehat{\beta}_{\mathcal{D}})\coloneqq\sum_{i=1}^{r}\lVert\varphi_{i}\otimes\psi_{i}\otimes
w_{i}\rVert=\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert.$

As noted in Section 2, only the variable multiplication in step (ii) counts in
bilinear complexity; the other two steps comprising scalar multiplications and
additions are discounted. In bilinear stability all three steps contribute to
the growth factor.

######  Proposition 3.2.

The minimal growth factor is given by nucler norm of the $\beta$, i.e.,

$\min_{\mathcal{D}}\gamma(\widehat{\beta}_{\mathcal{D}})=\lVert\beta\rVert_{\nu},$

with $\mathcal{D}$ running over all decomposition. Furthermore, there is
always an algorithm that attains the minimal growth factor.

The above equality is just stating (1.4) in terms of the growth factor. That
there is always an algorithm attaining the minimal growth factor, justifying
our writing $\min$ instead of $\inf$, follows from the existence of a _nuclear
decomposition_ [13, Proposition 3.1], i.e., a decomposition that attains the
nuclear norm. Just as a rank decomposition of $\beta$ represents a fastest
algorithm in bilinear complexity, a nuclear decomposition of $\beta$
represents a stablest algorithm in bilinear stability.

We next establish a rigorous relationship between growth factor and numerical
stability by proving a forward error bound in terms of the growth factor of a
bilinear algorithm. We assume a system of floating point arithmetic obeying
the standard model as in [17]: For $x,y\in\mathbb{R}$

$\operatorname{\mathsf{fl}}(x\mathbin{\mathsf{op}}y)=(x\mathbin{\mathsf{op}}y)(1+\delta),\quad\quad|\delta|\leq\mathsf{u},\quad\mathbin{\mathsf{op}}=+,-,*,/$
(3.3)

with $\mathsf{u}$ the unit roundoff. We assume that
$\mathbb{U},\mathbb{V},\mathbb{W}$ are vector spaces of dimensions $m,n,p$ and
that appropriate computational bases have been chosen on them so that we may
identify $\mathbb{U}\cong\mathbb{R}^{m}$, $\mathbb{V}\cong\mathbb{R}^{n}$,
$\mathbb{W}\cong\mathbb{R}^{p}$. The computational bases do not need to be the
standard bases and may instead be Fourier, Krylov, Haar, wavelet bases, etc.
This is another reason why we cast our discussions in terms of abstract vector
spaces and do not choose bases until absolutely necessary. However, once a
choice of bases has been made, the result below depends only on the dimensions
of $\mathbb{U},\mathbb{V},\mathbb{W}$; if say, $\mathbb{U}=\mathbb{R}^{m\times
n}$, then only the fact that it has dimension $mn$ matters, i.e.,
$\mathbb{U}\cong\mathbb{R}^{mn}$.

######  Theorem 3.3 (Growth factor and forward error).

Let $\beta\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}^{p}$ be a bilinear
operator, $\mathcal{D}=(\varphi_{i},\psi_{i},w_{i})_{i=1}^{r}$ a
decomposition, and $\widehat{\beta}_{\mathcal{D}}$ the corresponding
algorithm. If $\widehat{\beta}_{\mathcal{D}}(u,v)$ is the output of
$\widehat{\beta}_{\mathcal{D}}$ with $u\in\mathbb{R}^{m}$ and
$v\in\mathbb{R}^{n}$ as inputs, then

$\lVert\beta(u,v)-\widehat{\beta}_{\mathcal{D}}(u,v)\rVert_{\infty}\leq(m+n+r+1)\gamma(\widehat{\beta}_{\mathcal{D}})\lVert
u\rVert\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}).$

###### Proof.

We first show that the result reduces to the case $p=1$. It suffices to show
that

$\lvert\beta(u,v)_{k}-\widehat{\beta}_{\mathcal{D}}(u,v)_{k}\rvert\leq(m+n+r+1)\gamma(\widehat{\beta}_{\mathcal{D}})\lVert
u\rVert\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2})$ (3.4)

for all $k=1,\dots,p$, where the subscript $k$ refers to the $k$th coordinate
of a vector in $\mathbb{R}^{p}$. Since

$\gamma(\widehat{\beta}_{\mathcal{D}})=\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert\geq\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}|w_{ik}|,$

with $w_{ik}$ the $k$th coordinate of $w_{k}\in\mathbb{R}^{p}$, to show (3.4),
it suffices to show

$|\beta(u,v)_{k}-\widehat{\beta}_{\mathcal{D}}(u,v)_{k}|\leq(m+n+r+1)\Biggl{[}\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}|w_{ik}|\Biggr{]}\lVert
u\rVert\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}),$

which is equivalent to the case $p=1$. In the following, we will assume that
$p=1$.

Since $\varphi_{i}$ and $\psi_{i}$ are linear functionals on $\mathbb{R}^{m}$
and $\mathbb{R}^{n}$, there exist $u_{i}\in\mathbb{R}^{m}$ and
$v_{i}\in\mathbb{R}^{n}$ such that

$\varphi_{i}(u)=u_{i}^{\scriptscriptstyle\mathsf{T}}u\quad\textrm{and}\quad\psi_{i}(v)=v_{i}^{\scriptscriptstyle\mathsf{T}}v,$

for all $u\in\mathbb{R}^{m}$ and $v\in\mathbb{R}^{n}$. By [17, equation 3.7],

$|x^{\scriptscriptstyle\mathsf{T}}y-\operatorname{\mathsf{fl}}(x^{\scriptscriptstyle\mathsf{T}}y)|\leq
n|x|^{\scriptscriptstyle\mathsf{T}}|y|\mathsf{u}+O(\mathsf{u}^{2}),$

for any $x,y\in\mathbb{R}^{n}$ where $|\cdot|$ and $\leq$ apply
coordinatewise. So for each $i=1,\dots,r$,

$\begin{split}|\varphi_{i}(u)-\operatorname{\mathsf{fl}}(\varphi_{i}(u))|&=|u_{i}^{\scriptscriptstyle\mathsf{T}}u-\operatorname{\mathsf{fl}}(u_{i}^{\scriptscriptstyle\mathsf{T}}u)|\leq
m|u_{i}|^{\scriptscriptstyle\mathsf{T}}|u|\mathsf{u}+O(\mathsf{u}^{2})\\\
&\leq m\lVert u_{i}\rVert\lVert
u\rVert\mathsf{u}+O(\mathsf{u}^{2})=m\lVert\varphi_{i}\rVert_{*}\lVert
u\rVert\mathsf{u}+O(\mathsf{u}^{2}).\end{split}$ (3.5)

Likewise, for each $i=1,\dots,r$,

$|\psi_{i}(v)-\operatorname{\mathsf{fl}}(\psi_{i}(v))|\leq
n\lVert\psi_{i}\rVert_{*}\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}).$ (3.6)

Let $\Delta_{1,i}=\operatorname{\mathsf{fl}}(\varphi_{i}(u))-\varphi_{i}(u)$
and $\Delta_{2,i}=\operatorname{\mathsf{fl}}(\psi_{i}(v))-\psi_{i}(v)$. By
(3.5) and (3.6),

$|\Delta_{1,i}|\leq m\lVert\varphi_{i}\rVert_{*}\lVert
u\rVert\mathsf{u}+O(\mathsf{u}^{2}),\quad|\Delta_{2,i}|\leq
n\lVert\psi_{i}\rVert_{*}\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}).$ (3.7)

Let $c_{i}=\varphi_{i}(u)\psi_{i}(v)$ and $\widehat{c}_{i}$ be its computed
value. By (3.7), there exists $\delta_{i}$ with $|\delta_{i}|\leq\mathsf{u}$
such that

$\begin{split}\widehat{c}_{i}&=(\varphi_{i}(u)+\Delta_{1,i})(\psi_{i}(v)+\Delta_{2,i})(1+\delta_{i})\\\
&=\varphi_{i}(u)\psi_{i}(v)+\Delta_{1,i}\psi_{i}(v)+\varphi_{i}(u)\Delta_{2,i}+\delta_{i}\varphi_{i}(u)\psi_{i}(v)+O(\mathsf{u}^{2}).\end{split}$
(3.8)

By (3.7) and (3.8),

$\begin{split}\quad|c_{i}-\widehat{c}_{i}|&\leq
m\lVert\varphi_{i}\rVert_{*}\lVert
u\rVert|\psi_{i}(v)|\mathsf{u}+|\varphi_{i}(u)|n\lVert\psi_{i}\rVert_{*}\lVert
v\rVert\mathsf{u}+|\varphi_{i}(u)\psi_{i}(v)|\mathsf{u}+O(\mathsf{u}^{2})\\\
&\leq(m+n+1)\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
u\rVert\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}).\end{split}$ (3.9)

Let $\Delta_{i}=\widehat{c}_{i}-c_{i}$. By (3.9),

$|\Delta_{i}|\leq(m+n+1)\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
u\rVert\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}).$ (3.10)

Let $d_{i}=c_{i}w_{i}$ and $\widehat{d}_{i}$ be the computed value of $d_{i}$.
By (3.10), there exists $\delta_{i}^{\prime}$ with
$|\delta_{i}^{\prime}|\leq\mathsf{u}$ such that

$\widehat{d}_{i}=(c_{i}+\Delta_{i})w_{i}(1+\delta_{i}^{\prime})=c_{i}w_{i}+\Delta_{i}w_{i}+\delta_{i}^{\prime}c_{i}w_{i}+O(\mathsf{u}^{2}).$
(3.11)

Let $\Delta_{i}^{\prime}=\widehat{d}_{i}-d_{i}$. By (3.10) and (3.11),

$\begin{split}\quad|\Delta_{i}^{\prime}|&\leq(m+n+1)\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
u\rVert\lVert v\rVert\lVert
w_{i}\rVert\mathsf{u}+|\varphi_{i}(u)\psi_{i}(v)|\lVert
w_{i}\rVert\mathsf{u}+O(\mathsf{u}^{2})\\\
&\leq(m+n+2)\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert\lVert u\rVert\lVert
v\rVert\mathsf{u}+O(\mathsf{u}^{2}).\end{split}$ (3.12)

Finally, let $a=\sum_{i=1}^{r}\varphi_{i}(u)\psi_{i}(v)w_{i}$ and
$\widehat{a}$ be the computed value of $a$. By (3.12), there exists $\delta$
with $|\delta|\leq\mathsf{u}$ such that

$\begin{split}\widehat{a}&=\widehat{d}_{1}(1+\delta)^{r-1}+\widehat{d}_{2}(1+\delta)^{r-1}+\widehat{d}_{3}(1+\delta)^{r-2}+\dots+\widehat{d}_{r}(1+\delta),\end{split}$

where we compute the sum
$\widehat{d}_{1}+\widehat{d}_{2}+\dots\widehat{d}_{r}$ from left to right.
Hence we obtain

$\begin{split}|a-\widehat{a}|&\leq(m+n+2)\lVert u\rVert\lVert
v\rVert\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert\mathsf{u}+(r-1)\Bigg{\lVert}\sum_{i=1}^{r}c_{i}w_{i}\Bigg{\rVert}\mathsf{u}+O(\mathsf{u}^{2})\\\
&\leq(m+n+r+1)\lVert u\rVert\lVert
v\rVert\sum_{i=1}^{r}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
w_{i}\rVert\mathsf{u}+O(\mathsf{u}^{2})\\\
&=(m+n+r+1)\gamma(\widehat{\beta}_{\mathcal{D}})\lVert u\rVert\lVert
v\rVert\mathsf{u}+O(\mathsf{u}^{2}).\qed\end{split}$

Theorem 3.3 essentially says that that algorithms with small growth factors
have small forward errors. Combined with Proposition 3.2, we see that the
optimally stable algorithm in this context is the one corresponding to a
nuclear decomposition of $\beta$.

######  Corollary 3.4 (Tensor nuclear norm and forward error).

Let $\beta\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{p}$ be a bilinear
operator, $\mathcal{D}=(\varphi_{i},\psi_{i},w_{i})_{i=1}^{r}$ a nuclear
decomposition, and $\widehat{\beta}_{\mathcal{D}}$ the corresponding
algorithm. Then

$\lVert\beta(u,v)-\widehat{\beta}_{\mathcal{D}}(u,v)\rVert_{\infty}\leq(m+n+r+1)\lVert\beta\rVert_{\nu}\lVert
u\rVert\lVert v\rVert\mathsf{u}+O(\mathsf{u}^{2}).$

In principle, there is no reason to expect there to be an algorithm that is
both fastest in the sense of Section 2 and stablest in the sense of this
section, i.e., having a decomposition that attains both tensor rank and
nuclear norm. In Section 5, we will see that such an algorithm exists for
complex multiplication and we will study its properties when applied to
complex matrix multiplication.

##  4\. Fast matrix multiplications

As an illustration of bilinear stability in the last section, we will
calculate the growth factors of Strassen’s algorithm [30] and Winograd’s
variant [17, 22] for fast matrix multiplication and compare their stability
empirically. We will see that the growth factor of Strassen’s algorithm is
smaller than that of Winograd’s variant, and, consistent with the prediction
of Theorem 3.3, numerical experiments indeed show that the former gives more
accurate results.

###  4.1. Bilinear stability of Strassen multiplication

Given two block matrices

$A=\begin{bmatrix}A_{11}&A_{12}\\\ A_{21}&A_{22}\end{bmatrix},\quad
B=\begin{bmatrix}B_{11}&B_{12}\\\ B_{21}&B_{22}\end{bmatrix},$

Strassen’s algorithm [30] first computes

$\displaystyle M_{1}$ $\displaystyle=(A_{11}+A_{22})(B_{11}+B_{22}),\qquad$
$\displaystyle M_{5}$ $\displaystyle=(A_{11}+A_{12})B_{22},$ $\displaystyle
M_{2}$ $\displaystyle=(A_{21}+A_{22})B_{11},$ $\displaystyle M_{6}$
$\displaystyle=(A_{21}-A_{11})(B_{11}+B_{12}),$ $\displaystyle M_{3}$
$\displaystyle=A_{11}(B_{12}-B_{22}),$ $\displaystyle M_{7}$
$\displaystyle=(A_{12}-A_{22})(B_{21}+B_{22}),$ $\displaystyle M_{4}$
$\displaystyle=A_{22}(B_{21}-B_{11}),$

and then computes the product via

$AB=\begin{bmatrix}M_{1}+M_{4}-M_{5}+M_{7}&M_{3}+M_{5}\\\
M_{2}+M_{4}&M_{1}-M_{2}+M_{3}+M_{6}\end{bmatrix}.$

Note that this may be applied recursively. Let
$\widehat{\beta}_{\scriptscriptstyle\mathsf{S}}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{2\times 2}\times\mathbb{R}^{2\times
2}\to\mathbb{R}^{2\times 2}$ denote the Strassen’s algorithm for $2\times 2$
matrices. It is routine to check that for $A,B\in\mathbb{R}^{2\times 2}$,

$\widehat{\beta}_{\scriptscriptstyle\mathsf{S}}(A,B)=\sum_{i=1}^{7}\varphi_{i}(A)\psi_{i}(B)W_{i},$

where
$\varphi_{i}(A)=\operatorname{tr}(U_{i}^{\scriptscriptstyle\mathsf{T}}A)$ and
$\psi_{i}(B)=\operatorname{tr}(V_{i}^{\scriptscriptstyle\mathsf{T}}B)$ with

$\displaystyle U_{1}$ $\displaystyle=\begin{bmatrix}1&0\\\
0&1\end{bmatrix},\quad$ $\displaystyle V_{1}$
$\displaystyle=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix},\quad$ $\displaystyle
W_{1}$ $\displaystyle=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix};$ $\displaystyle
U_{2}$ $\displaystyle=\begin{bmatrix}0&0\\\ 1&1\end{bmatrix},\quad$
$\displaystyle V_{2}$ $\displaystyle=\begin{bmatrix}1&0\\\
0&0\end{bmatrix},\quad$ $\displaystyle W_{2}$
$\displaystyle=\begin{bmatrix}0&0\\\ 1&-1\end{bmatrix};$ $\displaystyle U_{3}$
$\displaystyle=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix},\quad$ $\displaystyle
V_{3}$ $\displaystyle=\begin{bmatrix}0&1\\\ 0&-1\end{bmatrix},\quad$
$\displaystyle W_{3}$ $\displaystyle=\begin{bmatrix}0&1\\\ 0&1\end{bmatrix};$
$\displaystyle U_{4}$ $\displaystyle=\begin{bmatrix}0&0\\\
0&1\end{bmatrix},\quad$ $\displaystyle V_{4}$
$\displaystyle=\begin{bmatrix}-1&0\\\ 1&0\end{bmatrix},\quad$ $\displaystyle
W_{4}$ $\displaystyle=\begin{bmatrix}1&0\\\ 1&0\end{bmatrix};$ $\displaystyle
U_{5}$ $\displaystyle=\begin{bmatrix}1&1\\\ 0&0\end{bmatrix},\quad$
$\displaystyle V_{5}$ $\displaystyle=\begin{bmatrix}0&0\\\
0&1\end{bmatrix},\quad$ $\displaystyle W_{5}$
$\displaystyle=\begin{bmatrix}-1&1\\\ 0&0\end{bmatrix};$ $\displaystyle U_{6}$
$\displaystyle=\begin{bmatrix}-1&0\\\ 1&0\end{bmatrix},\quad$ $\displaystyle
V_{6}$ $\displaystyle=\begin{bmatrix}1&1\\\ 0&0\end{bmatrix},\quad$
$\displaystyle W_{6}$ $\displaystyle=\begin{bmatrix}0&0\\\ 0&1\end{bmatrix};$
$\displaystyle U_{7}$ $\displaystyle=\begin{bmatrix}0&1\\\
0&-1\end{bmatrix},\quad$ $\displaystyle V_{7}$
$\displaystyle=\begin{bmatrix}0&0\\\ 1&1\end{bmatrix},\quad$ $\displaystyle
W_{7}$ $\displaystyle=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix}.$

For simplicity we will use the Frobenius norm on $\mathbb{R}^{2\times 2}$
since it is self dual. The growth factor of Strassen’s algorithm is then given
by

$\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{S}})=\sum_{i=1}^{7}\lVert\varphi_{i}\rVert_{*}\lVert\psi_{i}\rVert_{*}\lVert
W_{i}\rVert=\sum_{i=1}^{7}\lVert
U_{i}\rVert_{\scriptscriptstyle\mathsf{F}}\lVert
V_{i}\rVert_{\scriptscriptstyle\mathsf{F}}\lVert
W_{i}\rVert_{\scriptscriptstyle\mathsf{F}}\\\ =12+2\sqrt{2}\approx 14.83.$
(4.1)

###  4.2. Bilinear stability of Winograd multiplication

Winograd’s algorithm [17, 22] computes a different set of intermediate
quantities

$\displaystyle M^{\prime}_{1}$
$\displaystyle=(A_{21}+A_{22}-A_{11})(B_{11}+B_{22}-B_{12}),\qquad$
$\displaystyle M^{\prime}_{5}$ $\displaystyle=(A_{21}+A_{22})(B_{12}-B_{11}),$
$\displaystyle M^{\prime}_{2}$ $\displaystyle=A_{11}B_{11},$ $\displaystyle
M^{\prime}_{6}$ $\displaystyle=(A_{11}+A_{12}-A_{21}-A_{22})B_{22},$
$\displaystyle M^{\prime}_{3}$ $\displaystyle=A_{12}B_{21},$ $\displaystyle
M^{\prime}_{7}$ $\displaystyle=A_{22}(B_{11}+B_{22}-B_{12}-B_{21}),$
$\displaystyle M^{\prime}_{4}$ $\displaystyle=(A_{11}-A_{21})(B_{22}-B_{12}),$

and then compute the product via

$AB=\begin{bmatrix}M^{\prime}_{2}+M^{\prime}_{3}&M^{\prime}_{1}+M^{\prime}_{2}+M^{\prime}_{5}+M^{\prime}_{6}\\\
M^{\prime}_{1}+M^{\prime}_{2}+M^{\prime}_{4}-M^{\prime}_{7}&M^{\prime}_{1}+M^{\prime}_{2}+M^{\prime}_{4}+M^{\prime}_{5}\end{bmatrix}.$

Again this can be applied recursively. Let
$\widehat{\beta}_{\scriptscriptstyle\mathsf{W}}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{2\times 2}\times\mathbb{R}^{2\times
2}\to\mathbb{R}^{2\times 2}$ denote the Winograd’s algorithm for $2\times 2$
matrices. It is again routine to check that for $A,B\in\mathbb{R}^{2\times
2}$,

$\widehat{\beta}_{\scriptscriptstyle\mathsf{W}}(A,B)=\sum_{i=1}^{7}\varphi_{i}^{\prime}(A)\psi_{i}^{\prime}(B)W_{i}^{\prime},$

where
$\varphi_{i}^{\prime}(A)=\operatorname{tr}(U_{i}^{\prime{\scriptscriptstyle\mathsf{T}}}A)$
and
$\psi_{i}^{\prime}(B)=\operatorname{tr}(V_{i}^{\prime{\scriptscriptstyle\mathsf{T}}}B)$
with

$\displaystyle U_{1}^{\prime}$ $\displaystyle=\begin{bmatrix}-1&0\\\
1&1\end{bmatrix},\quad$ $\displaystyle V_{1}^{\prime}$
$\displaystyle=\begin{bmatrix}1&-1\\\ 0&1\end{bmatrix},\quad$ $\displaystyle
W_{1}^{\prime}$ $\displaystyle=\begin{bmatrix}0&1\\\ 1&1\end{bmatrix};$
$\displaystyle U_{2}^{\prime}$ $\displaystyle=\begin{bmatrix}1&0\\\
0&0\end{bmatrix},\quad$ $\displaystyle V_{2}^{\prime}$
$\displaystyle=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix},\quad$ $\displaystyle
W_{2}^{\prime}$ $\displaystyle=\begin{bmatrix}1&1\\\ 1&1\end{bmatrix};$
$\displaystyle U_{3}^{\prime}$ $\displaystyle=\begin{bmatrix}0&1\\\
0&0\end{bmatrix},\quad$ $\displaystyle V_{3}^{\prime}$
$\displaystyle=\begin{bmatrix}0&0\\\ 1&0\end{bmatrix},\quad$ $\displaystyle
W_{3}^{\prime}$ $\displaystyle=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix};$
$\displaystyle U_{4}^{\prime}$ $\displaystyle=\begin{bmatrix}1&0\\\
-1&0\end{bmatrix},\quad$ $\displaystyle V_{4}^{\prime}$
$\displaystyle=\begin{bmatrix}0&-1\\\ 0&1\end{bmatrix},\quad$ $\displaystyle
W_{4}^{\prime}$ $\displaystyle=\begin{bmatrix}0&0\\\ 1&1\end{bmatrix};$
$\displaystyle U_{5}^{\prime}$ $\displaystyle=\begin{bmatrix}0&0\\\
1&1\end{bmatrix},\quad$ $\displaystyle V_{5}^{\prime}$
$\displaystyle=\begin{bmatrix}-1&1\\\ 0&0\end{bmatrix},\quad$ $\displaystyle
W_{5}^{\prime}$ $\displaystyle=\begin{bmatrix}0&1\\\ 0&1\end{bmatrix};$
$\displaystyle U_{6}^{\prime}$ $\displaystyle=\begin{bmatrix}1&1\\\
-1&-1\end{bmatrix},\quad$ $\displaystyle V_{6}^{\prime}$
$\displaystyle=\begin{bmatrix}0&0\\\ 0&1\end{bmatrix},\quad$ $\displaystyle
W_{6}^{\prime}$ $\displaystyle=\begin{bmatrix}0&1\\\ 0&0\end{bmatrix};$
$\displaystyle U_{7}^{\prime}$ $\displaystyle=\begin{bmatrix}0&0\\\
0&1\end{bmatrix},\quad$ $\displaystyle V_{7}^{\prime}$
$\displaystyle=\begin{bmatrix}1&-1\\\ -1&1\end{bmatrix},\quad$ $\displaystyle
W_{7}^{\prime}$ $\displaystyle=\begin{bmatrix}0&0\\\ -1&0\end{bmatrix}.$

With respect to the Frobenius norm, the growth factor of Winograd’s algorithm
is

$\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{W}})=\sum_{i=1}^{7}\lVert\varphi_{i}^{\prime}\rVert_{*}\lVert\psi_{i}^{\prime}\rVert_{*}\lVert
W_{i}^{\prime}\rVert_{\scriptscriptstyle\mathsf{F}}=\sum_{i=1}^{7}\lVert
U_{i}^{\prime}\rVert_{\scriptscriptstyle\mathsf{F}}\lVert
V_{i}^{\prime}\rVert_{\scriptscriptstyle\mathsf{F}}\lVert
W_{i}^{\prime}\rVert_{\scriptscriptstyle\mathsf{F}}\\\
=7+4\sqrt{2}+3\sqrt{3}\approx 17.85.$ (4.2)

###  4.3. Bilinear stability of conventional matrix multiplication

For completeness we state the growth factor of the conventional algorithm for
matrix multiplication
$\widehat{\beta}_{\scriptscriptstyle\mathsf{C}}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{2\times 2}\times\mathbb{R}^{2\times
2}\to\mathbb{R}^{2\times 2}$,

$\widehat{\beta}_{\scriptscriptstyle\mathsf{C}}(A,B)=\sum_{i,j,k=1}^{2}\operatorname{tr}(E_{ij}^{\scriptscriptstyle\mathsf{T}}A)\operatorname{tr}(E_{jk}^{\scriptscriptstyle\mathsf{T}}B)E_{ik},$

where $E_{ij}\in\mathbb{R}^{2\times 2}$ denotes the standard basis matrix. Its
growth factor is easily seen to be

$\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{C}})=\sum_{i,j,k=1}^{2}\lVert
E_{ij}\rVert_{\scriptscriptstyle\mathsf{F}}\lVert
E_{jk}\rVert_{\scriptscriptstyle\mathsf{F}}\lVert
E_{ik}\rVert_{\scriptscriptstyle\mathsf{F}}=8.$

From (4.1) and (4.2), we see that

$\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{W}})>\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{S}})>\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{C}}).$
(4.3)

The first inequality will be verified in the numerical experiments below; the
second is consistent with the well-known fact [17] that Strassen’s algorithm
is less accurate than conventional multiplication.

###  4.4. Numerical experiments for fast matrix multiplications

By Theorem 3.3 and the sizes of the growth factors in (4.3), we expect
Strassen’s algorithm to give more accurate results than Winograd’s variant
since it has a smaller growth factor. We test this statement with random
matrices generated in three different ways: with (a) real entries drawn from
the uniform distribution on $[-1,1]$, (b) real entries drawn from the standard
normal distribution, (c) complex entries whose real and imaginary parts are
drawn from the uniform distribution on $[-1,1]$. In the last case, note that
our earlier discussions over $\mathbb{R}$ apply verbatim over $\mathbb{C}$
with the same growth factors.

(a) Real random matrices $\mathcal{U}[-1,1]$.

(b) Real random matrices $\mathcal{N}(0,1)$.

(c) Complex random matrices $\mathcal{U}[-1,1]+\mathcal{U}[-1,1]i$.

Figure 1. Accuracy of Strassen’s algorithm and Winograd’s variant.

In all cases, we compute $\widehat{\beta}_{\scriptscriptstyle\mathsf{S}}(A,B)$
and $\widehat{\beta}_{\scriptscriptstyle\mathsf{W}}(A,B)$ using Strassen’s
algorithm and Winograd’s variant respectively and compare the results against
the exact value $\beta(A,B)=AB$ computed using the Matlab symbolic toolbox.
From Figure 1, we see that Strassen’s algorithm is indeed more accurate than
Winograd’s variant, substantiating Theorem 3.3. Even though the $14.83$ growth
factor of Strassen’s algorithm appears to differ only moderately from the
$17.85$ growth factor of Winograd’s variant, the effect is magnified multifold
as a result of recursion — these algorithms are applied recursively to an
$n\times n$ matrix as a block $2\times 2$ matrix $\lfloor\log_{2}n\rfloor$
times. The conventional algorithm, which has a growth factor of $8$, is
included in these plots for comparison.

##  5\. Complex multiplication

As described towards the end of Section 2, complex multiplication is an
$\mathbb{R}$-bilinear operator
$\beta_{\mathbb{C}}\in\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb{R}^{2}$
when we identify $\mathbb{C}\cong\mathbb{R}^{2}$, with the standard basis
vectors in $\mathbb{R}^{2}$

$e_{1}=\begin{bmatrix}1\\\ 0\end{bmatrix},\qquad e_{2}=\begin{bmatrix}0\\\
1\end{bmatrix}$

corresponding to $1,i\in\mathbb{C}$. We write
$e_{1}^{*},e_{2}^{*}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{2}\to\mathbb{R}$ for the dual basis, i.e., linear
functionals with

$e_{1}^{*}\biggl{(}\begin{bmatrix}a\\\ b\end{bmatrix}\biggr{)}=a,\qquad
e_{2}^{*}\biggl{(}\begin{bmatrix}a\\\ b\end{bmatrix}\biggr{)}=b.$

We will denote the regular algorithm $(a+bi)(c+di)=(ac-bd)+i(bc+ad)$, Gauss’s
algorithm (1.1), and our new algorithm (1.5) by
$\widehat{\beta}_{\scriptscriptstyle\mathsf{R}},\widehat{\beta}_{\scriptscriptstyle\mathsf{G}},\widehat{\beta}_{\scriptscriptstyle\mathsf{N}}$
respectively. For easy reference,

$\displaystyle\widehat{\beta}_{\scriptscriptstyle\mathsf{R}}\biggl{(}\begin{bmatrix}a\\\
b\end{bmatrix},\begin{bmatrix}c\\\ d\end{bmatrix}\biggr{)}$
$\displaystyle=\begin{bmatrix}ac-bd\\\
bc+ad\end{bmatrix},\qquad\widehat{\beta}_{\scriptscriptstyle\mathsf{G}}\biggl{(}\begin{bmatrix}a\\\
b\end{bmatrix},\begin{bmatrix}c\\\ d\end{bmatrix}\biggr{)}=\begin{bmatrix}ac-
bd\\\ (a+b)(c+d)-ac-bd\end{bmatrix},$
$\displaystyle\widehat{\beta}_{\scriptscriptstyle\mathsf{N}}\biggl{(}\begin{bmatrix}a\\\
b\end{bmatrix},\begin{bmatrix}c\\\ d\end{bmatrix}\biggr{)}$
$\displaystyle=\begin{bmatrix}\frac{1}{2}\Bigl{[}\Bigl{(}a+\frac{1}{\sqrt{3}}b\Bigr{)}\Bigl{(}c+\frac{1}{\sqrt{3}}d\Bigr{)}+\Bigl{(}a-\frac{1}{\sqrt{3}}b\Bigr{)}\Bigl{(}c-\frac{1}{\sqrt{3}}d\Bigr{)}-\frac{8}{3}bd\Bigr{]}\\\\[8.61108pt]
\frac{i\sqrt{3}}{2}\Bigl{[}\Bigl{(}a+\frac{1}{\sqrt{3}}b\Bigr{)}\Bigl{(}c+\frac{1}{\sqrt{3}}d\Bigr{)}-\Bigl{(}a-\frac{1}{\sqrt{3}}b\Bigr{)}\Bigl{(}c-\frac{1}{\sqrt{3}}d\Bigr{)}\Bigr{]}\end{bmatrix}.$

They correspond to the decompositions

$\displaystyle\widehat{\beta}_{\scriptscriptstyle\mathsf{R}}$
$\displaystyle=(e_{1}^{*}\otimes e_{1}^{*}-e_{2}^{*}\otimes e_{2}^{*})\otimes
e_{1}+(e_{1}^{*}\otimes e_{2}^{*}+e_{2}^{*}\otimes e_{1}^{*})\otimes e_{2},$
(5.1) $\displaystyle\widehat{\beta}_{\scriptscriptstyle\mathsf{G}}$
$\displaystyle=(e_{1}^{*}+e_{2}^{*})\otimes(e_{1}^{*}+e_{2}^{*})\otimes
e_{2}+e_{1}^{*}\otimes e_{1}^{*}\otimes(e_{1}-e_{2})-e_{2}^{*}\otimes
e_{2}^{*}\otimes(e_{1}+e_{2}),$ (5.2)
$\displaystyle\widehat{\beta}_{\scriptscriptstyle\mathsf{N}}$
$\displaystyle=\frac{4}{3}\biggl{(}\biggl{[}\frac{\sqrt{3}}{2}e_{1}^{*}+\frac{1}{2}e_{2}^{*}\biggr{]}\otimes\biggl{[}\frac{\sqrt{3}}{2}e_{1}^{*}+\frac{1}{2}e_{2}^{*}\biggr{]}\otimes\biggl{[}\frac{1}{2}e_{1}+\frac{\sqrt{3}}{2}e_{2}\biggr{]}$
$\displaystyle\qquad+\biggl{[}\frac{\sqrt{3}}{2}e_{1}^{*}-\frac{1}{2}e_{2}^{*}\biggr{]}\otimes\biggl{[}\frac{\sqrt{3}}{2}e_{1}^{*}-\frac{1}{2}e_{2}^{*}\biggr{]}\otimes\biggl{[}\frac{1}{2}e_{1}-\frac{\sqrt{3}}{2}e_{2}\biggr{]}-e_{2}^{*}\otimes
e_{2}^{*}\otimes e_{1}\biggr{)}.$ (5.3)

###  5.1. Bilinear stability of complex multiplication algorithms

Recall from Section 2 that
$\operatorname{rank}(\beta_{\mathbb{C}})=3=\operatorname{\overline{\operatorname{rank}}}(\beta_{\mathbb{C}})$,
i.e., both Gauss’s algorithm and our new algorithm have optimal bilinear
complexity whether in the exact or approximate sense. One may also show that
$\beta_{\mathbb{C}}$ has nuclear norm [13, Lemma 6.1] is given by

$\lVert\beta_{\mathbb{C}}\rVert_{\nu}=4.$

The growth factor of the regular algorithm (5.1) attains this minimum value,

$\displaystyle\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{R}})$
$\displaystyle=\lVert e_{1}^{*}\rVert_{*}\lVert e_{1}^{*}\rVert_{*}\lVert
e_{1}\rVert+\lVert-e_{2}^{*}\rVert_{*}\lVert e_{2}^{*}\rVert_{*}\lVert
e_{1}\rVert+\lVert e_{1}^{*}\rVert_{*}\lVert e_{2}^{*}\rVert_{*}\lVert
e_{2}\rVert+\lVert e_{2}^{*}\rVert_{*}\lVert e_{1}^{*}\rVert_{*}\lVert
e_{2}\rVert$ $\displaystyle=4=\lVert\beta_{\mathbb{C}}\rVert_{\nu},$

as does our new algorithm (5.3),

$\displaystyle\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{N}})$
$\displaystyle=\frac{4}{3}\biggl{(}\biggl{\lVert}\frac{\sqrt{3}}{2}e_{1}^{*}+\frac{1}{2}e_{2}^{*}\biggr{\rVert}_{*}\biggl{\lVert}\frac{\sqrt{3}}{2}e_{1}^{*}+\frac{1}{2}e_{2}^{*}\biggr{\rVert}_{*}\biggl{\lVert}\frac{1}{2}e_{1}+\frac{\sqrt{3}}{2}e_{2}\biggr{\rVert}$
$\displaystyle\qquad+\biggl{\lVert}\frac{\sqrt{3}}{2}e_{1}^{*}-\frac{1}{2}e_{2}^{*}\biggr{\rVert}_{*}\biggl{\lVert}\frac{\sqrt{3}}{2}e_{1}^{*}-\frac{1}{2}e_{2}^{*}\biggr{\rVert}_{*}\biggl{\lVert}\frac{1}{2}e_{1}-\frac{\sqrt{3}}{2}e_{2}\biggr{\rVert}+\lVert
e_{2}^{*}\rVert_{*}\lVert e_{2}^{*}\rVert_{*}\lVert e_{1}\rVert\biggr{)}$
$\displaystyle=4=\lVert\beta_{\mathbb{C}}\rVert_{\nu},$

but not Gauss’s algorithm (5.2),

$\displaystyle\gamma(\widehat{\beta}_{\scriptscriptstyle\mathsf{G}})$
$\displaystyle=\lVert e_{1}^{*}+e_{2}^{*}\rVert_{*}\lVert
e_{1}^{*}+e_{2}^{*}\rVert_{*}\lVert e_{2}\rVert+\lVert
e_{1}^{*}\rVert_{*}\lVert e_{1}^{*}\rVert_{*}\lVert e_{1}-e_{2}\rVert+\lVert-
e_{2}^{*}\rVert_{*}\lVert e_{2}^{*}\rVert_{*}\lVert e_{1}+e_{2}\rVert$
$\displaystyle=2(1+\sqrt{2})>\lVert\beta_{\mathbb{C}}\rVert_{\nu}.$

So Gauss’s algorithm $\widehat{\beta}_{\scriptscriptstyle\mathsf{G}}$ is
faster (by bilinear complexity) but less stable (by bilinear stability) than
the regular algorithm. Our new algorithm
$\widehat{\beta}_{\scriptscriptstyle\mathsf{N}}$ on the other hand is optimal
in both measures, attaining both $\operatorname{rank}(\beta_{\mathbb{C}})$ and
$\lVert\beta_{\mathbb{C}}\rVert_{\nu}$.

We stress that numerical stability is too complicated an issue to be
completely covered by the simple framework of bilinear stability. For
instance, from the perspective of cancellation errors, our new algorithm also
suffers from the issue pointed out in [17, Section 23.2.4] for Gauss’s
algorithm. By choosing $z=w$ and $b=\sqrt{3}/a$, our algorithm (5.3) computes

$\frac{1}{2}\Bigl{[}\Bigl{(}a+\frac{1}{a}\Bigr{)}^{2}+\Bigl{(}a-\frac{1}{a}\Bigr{)}^{2}-\frac{8}{a^{2}}\Bigr{]}+\frac{i\sqrt{3}}{2}\Bigl{[}\Bigl{(}a+\frac{1}{a}\Bigr{)}^{2}-\Bigl{(}a-\frac{1}{a}\Bigr{)}^{2}\Bigr{]}\eqqcolon
x+iy.$

There will be cancellation error in the computed real part $\widehat{x}$ when
$\lvert a\rvert$ is small and likewise in the computed imaginary part
$\widehat{y}$ when $\lvert a\rvert$ is large. Nevertheless, as discussed in
[17, Section 23.2.4], the new algorithm (5.3) is still stable in the weaker
sense of having acceptably small $\lvert x-\widehat{x}\rvert/\lvert z\rvert$
and $\lvert y-\widehat{y}\rvert/\lvert z\rvert$ even if $\lvert
x-\widehat{x}\rvert/\lvert x\rvert$ or $\lvert y-\widehat{y}\rvert/\lvert
y\rvert$ might be large.

###  5.2. Error analysis of new algorithm applied to matrices

While using Gauss’s algorithm or our new algorithm for multiplying of complex
_numbers_ is a pointless overkill, they become useful when applied to the
multiplication of complex _matrices_. Note that any complex matrices
$A+iB,C+iD\in\mathbb{C}^{n\times n}$ may be multiplied via their real and
imaginary parts $A,B,C,D\in\mathbb{R}^{n\times n}$:

$(A+iB)(C+iD)=(AC-BD)+i[AD+BC],$ (5.4)

allowing us to focus our attention on designing algorithms for real matrix
products. In this regard, Gauss’s algorithm applied in the form

$(A+iB)(C+iD)=(AC-BD)+i[(A+B)(C+D)-AC-BD]$ (5.5)

reduces the number of real matrix products from four to three at the expense
of more matrix additions. This represents an enormous saving as matrix
products are invariably orders of magnitude more costly than matrix additions
— note that this statement remains true even if the exponent of matrix
mulitplication $\omega$ is $2$. Our new algorithm (1.5) likewise applies in
the form

$\displaystyle(A+iB)(C+iD)$
$\displaystyle=\frac{1}{2}\biggl{[}\biggl{(}A+\frac{1}{\sqrt{3}}B\biggr{)}\biggl{(}C+\frac{1}{\sqrt{3}}D\biggr{)}+\biggl{(}A-\frac{1}{\sqrt{3}}B\biggr{)}\biggl{(}C-\frac{1}{\sqrt{3}}D\biggr{)}-\frac{8}{3}BD\biggr{]}$
$\displaystyle\qquad+\frac{i\sqrt{3}}{2}\biggl{[}\biggl{(}A+\frac{1}{\sqrt{3}}B\biggr{)}\biggl{(}C+\frac{1}{\sqrt{3}}D\biggr{)}-\biggl{(}A-\frac{1}{\sqrt{3}}B\biggr{)}\biggl{(}C-\frac{1}{\sqrt{3}}D\biggr{)}\biggr{]},$
(5.6)

trading expensive matrix products for inexpensive scalar multiplications and
additions.

The following is an error analysis of (5.6), i.e., our new algorithm applied
to complex matrix multiplication. We emulate a similar analysis for Gauss’s
algorithm in [15, 17], assuming in particular that the real matrix
multiplications involved are performed using the conventional algorithm (as
opposed to Strassen’s or Winograd’s). We remind the reader that conventional
matrix multiplication has the simple error bound

$|AB-\operatorname{\mathsf{fl}}(AB)|\leq n|A||B|\mathsf{u}+O(\mathsf{u}^{2})$
(5.7)

for $A,B\in\mathbb{R}^{n\times n}$.

######  Theorem 5.1 (Error analysis for our new algorithm).

Let $(A+iB)(C+iD)=F+iG$ with $F,G\in\mathbb{R}^{n\times n}$ and let
$\widehat{F}_{\scriptscriptstyle\mathsf{N}},\widehat{G}_{\scriptscriptstyle\mathsf{N}}$
be computed via (5.6) in floating point arithmetic satisfying (3.3). Then

$\displaystyle|F-\widehat{F}_{\scriptscriptstyle\mathsf{N}}|$
$\displaystyle\leq(n+7)\biggl{(}|A|+\frac{1}{\sqrt{3}}|B|\biggr{)}\biggl{(}|C|+\frac{1}{\sqrt{3}}|D|\biggr{)}\mathsf{u}+\biggl{(}\frac{4}{3}n+4\biggr{)}|B||D|\mathsf{u}+O(\mathsf{u}^{2}),$
(5.8) $\displaystyle|G-\widehat{G}_{\scriptscriptstyle\mathsf{N}}|$
$\displaystyle\leq\sqrt{3}(n+6)\biggl{(}|A|+\frac{1}{\sqrt{3}}|B|\biggr{)}\biggl{(}|C|+\frac{1}{\sqrt{3}}|D|\biggr{)}\mathsf{u}+O(\mathsf{u}^{2}),$
(5.9)

where the inequality $\leq$ and absolute value $\lvert\,\cdot\,\rvert$ both
apply in a coordinatewise sense.

###### Proof.

Following [17], we use the same letter $\delta$ to denote the error incurred
in each step of our algorithm. So, for example,

$\operatorname{\mathsf{fl}}(B/\sqrt{3})=B/\sqrt{3}+\delta B/\sqrt{3}.$

In the following we will define matrices $H_{i}$ and let $\widehat{H}_{i}$ be
its computed value, $i=1,\dots,8$.

Let $H_{1}\coloneqq A+B/\sqrt{3}$. Then

$\displaystyle\widehat{H}_{1}$
$\displaystyle=\operatorname{\mathsf{fl}}(A+B/\sqrt{3}+\delta
B/\sqrt{3})=(A+B/\sqrt{3}+\delta B/\sqrt{3})(1+\delta)$
$\displaystyle=A+B/\sqrt{3}+\delta(A+2B/\sqrt{3})+O(\mathsf{u}^{2})$
$\displaystyle=H_{1}+2\Delta_{1}+O(\mathsf{u}^{2}),\qquad|\Delta_{1}|\leq(|A|+|B|/\sqrt{3})\mathsf{u}.$

Similarly $H_{2}\coloneqq C+D/\sqrt{3}$ satisfies

$\widehat{H}_{2}=H_{2}+2\Delta_{2}+O(\mathsf{u}^{2}),\qquad|\Delta_{2}|\leq(|C|+|D|/\sqrt{3})\mathsf{u}.$

Let $H_{3}\coloneqq(A+B/\sqrt{3})(C+D/\sqrt{3})$. By (5.7),

$\widehat{H}_{3}=(A+B/\sqrt{3}+2\Delta_{1})(C+D/\sqrt{3}+2\Delta_{2})+n\Delta_{3}+O(\mathsf{u}^{2})$
(5.10)

where

$\displaystyle|\Delta_{3}|$
$\displaystyle\leq|(A+B/\sqrt{3}+2\Delta_{1})||(C+D/\sqrt{3}+2\Delta_{2})|\mathsf{u}$
(5.11)
$\displaystyle\leq(|A|+|B|/\sqrt{3}+2|\Delta_{1}|)(|C|+|D|/\sqrt{3}+2|\Delta_{2}|)\mathsf{u}$
$\displaystyle\leq(|A|+|B|/\sqrt{3}+2\mathsf{u}(|A|+|B|/\sqrt{3}))(|C|+|D|/\sqrt{3}+2\mathsf{u}(|C|+|D|/\sqrt{3}))\mathsf{u}$
$\displaystyle\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u}+O(\mathsf{u}^{2}).$

By (5.10) and (5.11),

$\displaystyle\widehat{H}_{3}$
$\displaystyle=(A+B/\sqrt{3})(C+D/\sqrt{3})+2\Delta_{1}(C+D/\sqrt{3})$ (5.12)
$\displaystyle\qquad+2(A+B/\sqrt{3})\Delta_{2}+n\Delta_{3}+O(\mathsf{u}^{2})$
$\displaystyle=H_{3}+(n+4)\Delta_{4}+O(\mathsf{u}^{2})$

where

$|\Delta_{4}|\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u}.$

Similarly $H_{4}\coloneqq(A-B/\sqrt{3})(C-D/\sqrt{3})$ satisfies

$\widehat{H}_{4}=H_{4}+(n+4)\Delta_{5}+O(\mathsf{u}^{2})$ (5.13)

where

$|\Delta_{5}|\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u}.$

Let $H_{5}\coloneqq(A+B/\sqrt{3})(C+D/\sqrt{3})+(A-B/\sqrt{3})(C-D/\sqrt{3})$.
By (5.12) and (5.13),

$\displaystyle\widehat{H}_{5}$
$\displaystyle=[H_{3}+(n+4)\Delta_{4}+H_{4}+(n+4)\Delta_{5}](1+\delta)+O(u^{2})$
(5.14) $\displaystyle=H_{5}+(2n+10)\Delta_{6}+O(\mathsf{u}^{2})$

where

$|\Delta_{6}|\leq\mathsf{u}(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3}).$

Let $H_{6}\coloneqq 8/3BD$. Then

$\displaystyle\widehat{H}_{6}$
$\displaystyle=\operatorname{\mathsf{fl}}(8/3(BD+n\Delta_{7}))+O(\mathsf{u}^{2})$
(5.15) $\displaystyle=8/3(BD+n\Delta_{7})(1+\delta)+O(\mathsf{u}^{2})$
$\displaystyle=H_{6}+8/3(n+1)\Delta_{8}+O(\mathsf{u}^{2})$

where

$|\Delta_{7}|\leq|B||D|\mathsf{u},\qquad|\Delta_{8}|\leq|B||D|\mathsf{u}.$

Let $H_{7}\coloneqq H_{5}-H_{6}$. By (5.14) and (5.15),

$\displaystyle\widehat{H}_{7}$
$\displaystyle=[H_{5}+(2n+10)\Delta_{6}-H_{6}-8/3(n+1)\Delta_{8}](1+\delta)+O(\mathsf{u}^{2})$
$\displaystyle=H_{7}+(2n+12)\Delta_{9}+8/3(n+2)\Delta_{10}+O(\mathsf{u}^{2})$

where

$|\Delta_{9}|\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u},\qquad|\Delta_{10}|\leq|B||D|\mathsf{u}.$

Then

$\displaystyle\widehat{F}_{\scriptscriptstyle\mathsf{N}}$
$\displaystyle=(1+\delta)[H_{7}+(2n+12)\Delta_{9}+8/3(n+2)\Delta_{10}]/2+O(\mathsf{u}^{2})$
$\displaystyle=F+(n+7)\Delta_{11}+4/3(n+3)\Delta_{12}+O(\mathsf{u}^{2})$

where

$|\Delta_{11}|\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u},\qquad|\Delta_{12}|\leq|B||D|\mathsf{u},$

and from which we obtain (5.8).

Let $H_{8}\coloneqq(A+B/\sqrt{3})(C+D/\sqrt{3})-(A-B/\sqrt{3})(C-D/\sqrt{3})$.
Similar to (5.14), we have

$\widehat{H}_{8}=H_{8}-(2n+10)\Delta_{13}+O(\mathsf{u}^{2})$

where

$|\Delta_{13}|\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u}.$

Then

$\displaystyle\widehat{G}_{\scriptscriptstyle\mathsf{N}}$
$\displaystyle=\sqrt{3}/2[H_{8}-(2n+10)\Delta_{13}](1+\delta)+O(\mathsf{u}^{2})$
$\displaystyle=G+\sqrt{3}(n+6)\Delta_{14}+O(\mathsf{u}^{2})$

where

$|\Delta_{14}|\leq(|A|+|B|/\sqrt{3})(|C|+|D|/\sqrt{3})\mathsf{u},$

from which we obtain (5.9). ∎

If we compute the matrices $F,G$ in Theorem 5.1 using Gauss’s algorithm (5.5)
with floating point arithmetic and let the results be
$\widehat{F}_{\scriptscriptstyle\mathsf{G}}$ and
$\widehat{G}_{\scriptscriptstyle\mathsf{G}}$, then the corresponding error
bounds [15, 17] are

$\displaystyle|F-\widehat{F}_{\scriptscriptstyle\mathsf{G}}|$
$\displaystyle\leq(n+1)(|A||C|+|B||D|)\mathsf{u}+O(\mathsf{u}^{2}),$ (5.16)
$\displaystyle|G-\widehat{G}_{\scriptscriptstyle\mathsf{G}}|$
$\displaystyle\leq(n+4)[(|A|+|B|)(|C|+|D|)+|A||C|+|B||D|]\mathsf{u}+O(\mathsf{u}^{2}).$

When $n\to\infty$, we have $n+c\approx n$ for any constant $c$. Hence the
errors in (5.8) and (5.9) are dominated by

$\displaystyle|F-\widehat{F}_{\scriptscriptstyle\mathsf{N}}|$
$\displaystyle\sim
n\biggl{[}|A||C|+\frac{5}{3}|B||D|+\frac{1}{\sqrt{3}}|A||D|+\frac{1}{\sqrt{3}}|B||C|\biggr{]}\mathsf{u},$
$\displaystyle|G-\widehat{G}_{\scriptscriptstyle\mathsf{N}}|$
$\displaystyle\sim
n\biggl{[}\sqrt{3}|A||C|+|B||C|+|A||D|+\frac{1}{\sqrt{3}}|B||D|\biggr{]}\mathsf{u},$

whereas those in (5.16) are dominated by

$\displaystyle|F-\widehat{F}_{\scriptscriptstyle\mathsf{G}}|$
$\displaystyle\sim n(|A||C|+|B||D|)\mathsf{u},$
$\displaystyle|G-\widehat{G}_{\scriptscriptstyle\mathsf{G}}|$
$\displaystyle\sim n(2|A||C|+2|B||D|+|A||D|+|B||C|)\mathsf{u}.$

For easy comparison suppose the magnitudes of the entries in $A,B,C,D$ are all
approximately $\theta$, then these reduce to

$\displaystyle|F-\widehat{F}_{\scriptscriptstyle\mathsf{N}}|$
$\displaystyle\sim 3.8n^{2}\theta^{2},\qquad$
$\displaystyle|G-\widehat{G}_{\scriptscriptstyle\mathsf{N}}|$
$\displaystyle\sim 4.3n^{2}\theta^{2},$ (5.17)
$\displaystyle|F-\widehat{F}_{\scriptscriptstyle\mathsf{G}}|$
$\displaystyle\sim 2n^{2}\theta^{2},$
$\displaystyle|G-\widehat{G}_{\scriptscriptstyle\mathsf{G}}|$
$\displaystyle\sim 6n^{2}\theta^{2}.$

So Gauss’s algorithm gives an imaginary part that is three times less accurate
than its real part. Note the the imaginary part of Gauss’s algorithm accounts
for all its computational savings; the real part is just the regular
algorithm. On the other hand, our algorithm balances the accuracy of both the
real and imaginary parts by spreading out the computational savings across
both parts.

To quantify this, we use the _max norm_. For a complex matrix
$A+iB\in\mathbb{C}^{n\times n}$, this is

$\lVert A+iB\rVert_{\max}\coloneqq\max\\{\lvert a_{ij}\rvert,\lvert
b_{ij}\rvert\mathrel{\mathop{\mathchar 58\relax}}i,j=1,\dots,n\\}.$ (5.18)

The max norm differs from the usual matrix $\infty$-norm given by maximum row
sum used in [15, 17]. We favor the max norm as it is the strictest measure of
numerical accuracy — a small max norm error implies that each entry is
accurate as opposed to accurate on average.

If we denote the matrices resulting from Gauss’s algorithm and our new
algorithm by

$\widehat{E}_{\scriptscriptstyle\mathsf{G}}\coloneqq\widehat{F}_{\scriptscriptstyle\mathsf{G}}+i\widehat{G}_{\scriptscriptstyle\mathsf{G}},\qquad\widehat{E}_{\scriptscriptstyle\mathsf{N}}\coloneqq\widehat{F}_{\scriptscriptstyle\mathsf{N}}+i\widehat{G}_{\scriptscriptstyle\mathsf{N}}$

respectively, we expect $\lVert
E-\widehat{E}_{\scriptscriptstyle\mathsf{N}}\rVert_{\max}$ to be smaller than
$\lVert E-\widehat{E}_{\scriptscriptstyle\mathsf{G}}\rVert_{\max}$. The
extensive experiments in the next section will attest to this.

##  6\. Experiments for new complex matrix multiplication algorithm

The goal of this section is to provide numerical evidence to show that our new
algorithm (5.6) for complex matrix multiplication is

  * •

nearly as accurate as the regular algorithm (5.4), and

  * •

nearly as fast as Gauss’s algorithm (5.5).

We begin with routine experiments comparing the three algorithms (5.4), (5.5),
(5.6) on random matrices, and move on to three actual applications: matrix
polynomial evaluations, unitary transformations, and the increasingly popular
complex-valued neural networks. The results, we think, show that our new
algorithm can be a realistic replacment for Gauss’s algorithm in engineering
applications.

###  6.1. Speed of the algorithms

We generate random $A+iB,C+iD\in\mathbb{C}^{n\times n}$ with entries of
$A,B,C,D$ drawn uniformly in $[-1,1]$; the results with standard normal are
similar and omitted. We increase $n$ from $2100$ to $7000$ in steps of $100$.
The product $(A+iB)(C+iD)$ is computed numerically with the regular algorithm
(5.4), Gauss’s algorithm (5.5), and our new algorithm (5.6). For each $n$, we
generate ten different matrices and record the average time taken for each
algorithm and plot these in Figure 2, with wall time (in seconds) for vertical
axis and $\log_{10}(n)$ for horizontal axis. The time taken by Matlab’s
internal function for complex matrix multiplication is virtually
indistinguishable from that of the regular algorithm and therefore omitted.

Consistent with the predictions of bilinear complexity, our new algorithm has
roughly the same computation time as Gauss’s algorithm, at roughly $3/4$ the
time taken by the regular algorithm. We will perform more speed experiments in
conjunction with our accuracy experiments in Section 6.2.

Figure 2. Speed of the three algorithms for complex matrix multiplication.

###  6.2. Accuracy of the algorithms

Figure 3. Accuracy and speed of algorithms for complex matrix multiplication.

We generate random $A+iB,C+iD\in\mathbb{C}^{n\times n}$ with $n=64,128,256$
and with condition numbers ranging from $174$ to $3\times 10^{11}$. It is
desirable to limit ourselves to matrices over Gaussian rationals, i.e.,
$\mathbb{Q}+\mathbb{Q}i$, as we will need to compute the exact values of their
products later.

The way we generate such a matrix requires some elaboration. For an
$X\in\mathbb{Z}^{n\times n}$ with a specified
$\kappa_{2}(X)=\kappa\in\mathbb{Z}$. We form a diagonal
$\Lambda\in\mathbb{R}^{n\times n}$ whose diagonal entries are $1$ and $\kappa$
toegether with $n-2$ other random integers between $1$ and $\kappa-1$. We then
form $X=H\Lambda H^{\scriptscriptstyle\mathsf{T}}$ with a random Hadamard
matrix $H\in\mathbb{Z}^{n\times n}$. If $A$ and $B$ are generated in this
manner, then they are dense matrices (important as we do not want sparsity to
unduly influence arithmetic costs) and
$\kappa_{2}(A+iB)=\kappa_{2}(A)=\kappa_{2}(B)=\kappa$ as $(\kappa+\kappa
i)/(1+i)=\kappa$.

We compute the exact value of $(A+iB)(C+iD)$ symbolically with Matlab’s
symbolic toolbox. Given our relatively modest computational resources, this is
the bottleneck for our experiments as this step becomes prohibitively
expensive when $n>256$. In generating the $n=256$ plots in Figure 3, this step
alone took 40 hours on our University’s Research Computing Center servers.

For each pair of complex matrices $A+iB$ and $C+iD$, we compute their product
$\widehat{E}$ using each of the three algorithms (5.4), (5.5), (5.6), and
compare them against the exact result $E$ via the max norm relative error

$\frac{\lVert E-\widehat{E}\rVert_{\max}}{\lVert A+iB\rVert_{\max}\lVert
C+iD\rVert_{\max}}.$

As discussed in [15, 17], it is natural to measure error in matrix
multiplication relative to the norms of the input matrices. We use the max
norm in (5.18) to better capture entrywise accuracy.

The results are plotted in Figure 3: speed plots have wall time in seconds on
the vertical axes; accuracy plots have relative error on the vertical axes;
all plots have $\log_{10}(\kappa)$ on the horizontal axes. We repeat each
experiment ten times: every value on these plots comes from averaging across
the results of ten pairs of random matrices with the same condition number.

Observations from Figure 3: The accuracy of our new algorithm is much higher
than that of Gauss’s algorithm and only slightly worse than that of the
regular algorithm. Gauss’s algorithm also shows a great deal more fluctuation
across varying condition numbers than either our new algorithm or the regular
one. When it comes to speed, our algorithm is closer to that of Gauss’s than
the regular algorithm. These accuracy results attest to Theorem 5.1 and the
discussions around (5.17).

The relative errors and wall times for Matlab’s internal function for complex
matrix multiplication are virtually indistinguishable from those of the
regular algorithm (that we implemented ourselves) and thus omitted. In the
next three sections, we will compare the accuracy and speed of the three
complex matrix multiplication algorithms in more realstic scenarios.

###  6.3. Matrix polynomial evaluations

We evaluate a polynomial $p(x)=\sum_{k=0}^{d}a_{k}x^{k}$ with coefficients
$a_{0},\dots,a_{k}\in\mathbb{R}$ at a $X\in\mathbb{C}^{n\times n}$. This is a
problem that occurs in many tasks involving matrix functions [17, 18]. We
limit ourselves to real coefficients as this is by far most common scenario
[18]; but the complex coefficients case simply reduces to evaluating two real
polynomials $\operatorname{Re}p(x)$ and $\operatorname{Im}p(x)$. The
celebrated Horner’s rule [18, Algorithm 4.3] reduces the problem to one of
repeated matrix multiplications.

Algorithm 1 Compute $p(X)$ via Horner’s rule

1:$a_{0},a_{1},\dots,a_{d}\in\mathbb{R}$, $X\in\mathbb{C}^{n\times n}$

2:$a_{0}I+a_{1}X+\dots+a_{d}X^{d}$

3:$P=X$;

4:$S=a_{0}I+a_{1}X$;

5:for $k=2\mathrel{\mathop{\mathchar 58\relax}}d$ do

6: $P=PX$;

7: $S=S+a_{k}P$;

8:end for

9:return $S$;

We generate random matrices $X\in\mathbb{C}^{256\times 256}$ with condition
numbers from $2^{34}$ to $2^{53}$ as described in Section 6.2. We set $d=5$
and choose random $b_{0},\dots,b_{5}\in(0,1)$ uniformly. We then evaluate
$p(X)$ using Algorithm 1, with Step 4 computed via (5.4), (5.5), and (5.6).

Figure 4. The three algorithms applied to matrix polynomial evaluations.

We measure accuracy in terms of the max norm relative forward error

$\frac{\lVert p(X)-\widehat{p}(X)\rVert_{\max}}{\lVert p(X)\rVert_{\max}},$

using Matlab symbolic toolbox for the exact value of $p(X)$. The results
presented in Figure 4 again show that our new algorithm is nearly as accurate
as the regular algorithm and nearly as fast as Gauss’s algorithm. While our
accuracy tests are again limited by our capacity for symbolic computation
($n=256$ is fine, $n=512$ is beyond reach), our speed tests can go far beyond
(to around $n=4096$), and they show a profile much like Figure 2.

###  6.4. Unitary transforms

Given a unitary matrix $U\in\mathbb{C}^{n\times n}$ and a complex matrix
$X\in\mathbb{C}^{n\times n}$, it may come as a surprise to the reader that
unless $U$ happens to be some special transforms like FFT, DCT, DWT, etc, or
has already been factored into a product of Householder or Givens matrices,
there is no known special algorithm for forming $UX$ that would take advantage
of the unitarity of $U$. Nevertheless, such unitary matrices with no
additional special structure are not uncommon. For instance, the matrix $U$
could come from polar decompositions or matrix sign functions [14, 16, 20],
and computed via iterative methods [14, 16, 20] and thus not in Householder-
or Givens-factored form. Here we will explore the use of algorithms (5.4),
(5.5), (5.6) for unitary transforms $X\mapsto UX$.

We generate the unitary matrix $U\in\mathbb{C}^{256\times 256}$ by QR
factoring complex random matrices with entries in
$\mathcal{U}[0,1]+\mathcal{U}[0,1]i$. Note that a unitary matrix is always
perfectly conditioned. The matrix $X\in\mathbb{C}^{256\times 256}$ is
generated randomly with condition numbers from $2^{34}$ to $2^{53}$ as in
Section 6.3. We compute the exact value $E\coloneqq UX$ symbolically as before
and measure the accuracy of our computed value $\widehat{E}$ by

$\frac{\lVert E-\widehat{E}\rVert_{\max}}{\lVert U\rVert_{\max}\lVert
X\rVert_{\max}}.$

The results, presented in Figure 5, allow us to draw the same conclusion as in
the Section 6.3. Further speed tests up to $n=4096$ again show a profile much
like Figure 2.

Figure 5. The three algorithms applied to unitary transforms.

###  6.5. Complex-valued neural networks

Given the amount of recent attention it has received, a feedforward neural
network probably requires no introduction, and we will just limit ourselves to
the obligatory Figure 6.

[height=4,layerspacing=20mm] [count=4, bias=false, title=input
layer, text=$x_{\hiddenlayer}$[count=4, bias=false, title=hidden
layer 1, text=$h^{(1)}_{\linklayers}$[count=4, bias=false, title=hidden
layer 2, text=$h^{(2)}_{\linklayers}$[count=4, bias=false, title=hidden
layer 3, text=$h^{(3)}_{\linklayers}$[count=4, bias=false, title=hidden
layer 4, text=$h^{(4)}_{\linklayers}$[count=4, bias=false, title=hidden
layer 5, text=$h^{(5)}_{\linklayers}$[count=4, title=output
layer, text=$\widehat{y}_{\linklayers}$

Figure 6. A constant width neural network with input dimension $n=4$ and depth
$d=6$. The edges between adjacent layers are weighted with weight matrices.

A complex-valued neural networks is simply a neural network with complex-
valued weights and is activated by a complex function. It has become
increasingly important and is widely used in signal processing and computer
vision [1, 3, 9, 29, 34, 37]. For simplicity, we consider a $d$-layer constant
width version $f\mathrel{\mathop{\mathchar
58\relax}}\mathbb{C}^{n}\to\mathbb{C}^{n}$ given by

$f(W_{1},\dots,W_{d},\sigma)(x)\coloneqq W_{d}\sigma(W_{d-1}\sigma(\cdots
W_{2}\sigma(W_{1}x)\cdots)),$

with weight matrices $W_{1},\dots,W_{d}\in\mathbb{C}^{n\times n}$ and
activation function $\sigma\mathrel{\mathop{\mathchar
58\relax}}\mathbb{C}\to\mathbb{C}$ applied coordinatewise on $\mathbb{C}^{n}$.
Evaluating a trained complex-valued neural network on multiple inputs
$x_{1},\dots,x_{m}\in\mathbb{C}^{n}$ is an indispensable task when we employ
it to make new predictions. Here we will compare the performance of the three
algorithms (5.4), (5.5), (5.6) when used for this purpose.

For concreteness, we choose a depth of $d=6$ and use the complex ReLU
activation [3, 34]

$\sigma(a+bi)\coloneqq\max(a,0)+\max(b,0)i.$

We generate random weight matrices $W_{1},\dots,W_{6}\in\mathbb{C}^{n\times
n}$ with $n=64$ and $128$, and with condition numbers ranging from $2^{34}$ to
$2^{53}$. We also generate random inputs
$X=[x_{1},\dots,x_{m}]\in\mathbb{C}^{n\times m}$ with entries drawn from
$\mathcal{U}[-\frac{1}{2},\frac{1}{2}]+\mathcal{U}[-\frac{1}{2},\frac{1}{2}]i$,
and with $(m,n)=(25,64)$ or $(50,128)$. The task is then to evaluate

$E\coloneqq f(W_{1},\dots,W_{d},\sigma)(X)\coloneqq
W_{d}\sigma(W_{d-1}\sigma(\cdots W_{2}\sigma(W_{1}X)\cdots)).$

Again we compute its exact value $E$ symbolically, apply the three algorithms
to obtain $\widehat{E}$ numerically, and measure accuracy in terms of the
relative forward error

$\frac{\lVert E-\widehat{E}\rVert_{\max}}{\lVert E\rVert_{\max}}.$

The results, shown in Figure 7, are fully consistent with those in Sections
6.3 and 6.4.

Figure 7. The three algorithms applied to $6$-layer complex neural networks
with complex ReLU activation and widths $64$ and $128$.

##  7\. Conclusion

The notion of bilinear complexity started by Strassen has been a great
motivator for more than five decades of exciting developments in numerical
linear algebra. Its success illustrates the adage that “less is more”.
Bilinear complexity does not capture every operation that underlies the speed
of an algorithm; but by focusing on a single operation (variable
multiplications) and disregarding the rest (e.g., scalar multiplications,
additions), it allows speed to be measured by the number of terms in a
decomposition of a $3$-tensor and the fastest algorithm to be given by a rank
decomposition. This opens a door to other areas of mathematics like algebraic
geometry where such decompositions are studied independent of their
computational relevance.

We hope the notion of bilinear stability proposed in this article would do for
the study of numerical stability what bilinear complexity did for the study of
time complexity. By focusing on a single factor (growth) and disregarding
other factors (e.g., cancellation errors) that play a role in numerical
stability, it allows stability to be measured by the growth factor in a
decomposition of a $3$-tensor and the stablest algorithm to be given by a
nuclear decomposition. Just as tensor rank connects to algebraic geometry,
tensor nuclear norm connects to functional analysis [10, 11, 28]; thus
bilinear stability could potentially open a door to this rich area of
mathematics.

### Acknowledgment

The authors would like to thank Nick Higham and Ke Ye for helpful discussions
and University of Chicago’s Research Computing Center for its services.

## References

  * [1] I. Aizenberg.  Complex-valued neural networks with multi-valued neurons, volume 353 of Studies in Computational Intelligence.  Springer-Verlag, Berlin, 2011.
  * [2] G. Ballard, A. R. Benson, A. Druinsky, B. Lipshitz, and O. Schwartz.  Improving the numerical stability of fast matrix multiplication.  SIAM J. Matrix Anal. Appl., 37(4):1382–1418, 2016.
  * [3] J. Bassey, L. Qian, and X. Li.  A survey of complex-valued neural networks.  arXiv preprint arXiv:2101.12249, 2021.
  * [4] D. Bini and G. Lotti.  Stability of fast algorithms for matrix multiplication.  Numer. Math., 36(1):63–72, 1980/81.
  * [5] D. Bini, G. Lotti, and F. Romani.  Approximate solutions for the bilinear form computational problem.  SIAM J. Comput., 9(4):692–697, 1980.
  * [6] A. Borodin and I. Munro.  The computational complexity of algebraic and numeric problems.  Elsevier Computer Science Library: Theory of Computation Series, No. 1\. American Elsevier Publishing Co., Inc., New York-London-Amsterdam, 1975.
  * [7] R. P. Brent.  Algorithms for matrix multiplication.  March 1970.  Report Stan-CS-70-157, Stanford University.
  * [8] P. Bürgisser, M. Clausen, and M. A. Shokrollahi.  Algebraic complexity theory, volume 315 of Grundlehren der mathematischen Wissenschaften.  Springer-Verlag, Berlin, 1997.
  * [9] B. N. Dash and N. Khare.  Deep complex neural network applications in remote sensing: an introductory review.  In K. I. Ranney and A. M. Raynal, editors, Radar Sensor Technology XXV, volume 11742, pages 34 – 44. International Society for Optics and Photonics, SPIE, 2021.
  * [10] A. Defant and K. Floret.  Tensor norms and operator ideals, volume 176 of North-Holland Mathematics Studies.  North-Holland Publishing Co., Amsterdam, 1993.
  * [11] J. Diestel, J. H. Fourie, and J. Swart.  The metric theory of tensor products.  American Mathematical Society, Providence, RI, 2008.  Grothendieck’s résumé revisited.
  * [12] A. T. Fam.  Efficient complex matrix multiplication.  IEEE Trans. Comput., 37(7):877–879, 1988.
  * [13] S. Friedland and L.-H. Lim.  Nuclear norm of higher-order tensors.  Math. Comp., 87(311):1255–1281, 2018.
  * [14] N. J. Higham.  Computing the polar decomposition—with applications.  SIAM J. Sci. Statist. Comput., 7(4):1160–1174, 1986.
  * [15] N. J. Higham.  Stability of a method for multiplying complex matrices with three real matrix multiplications.  SIAM J. Matrix Anal. Appl., 13(3):681–687, 1992.
  * [16] 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.
  * [17] N. J. Higham.  Accuracy and stability of numerical algorithms.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2002.
  * [18] N. J. Higham.  Functions of matrices.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008.  Theory and computation.
  * [19] A. Karatsuba and Y. Ofman.  Multiplication of many-digital numbers by automatic computers.  Dokl. Akad. Nauk SSSR, 14(145):293–294, 1962.
  * [20] 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.
  * [21] V. V. Kljuev and N. I. Kokovkin-Ščerbak.  On the minimization of the number of arithmetic operations for solving linear algebraic systems of equations.  Ž. Vyčisl. Mat i Mat. Fiz., 5:21–33, 1965.
  * [22] D. E. Knuth.  The art of computer programming. Vol. 2.  Addison-Wesley, Reading, MA, 1998.  Seminumerical algorithms, Third edition [of MR0286318].
  * [23] J. M. Landsberg.  The border rank of the multiplication of $2\times 2$ matrices is seven.  J. Amer. Math. Soc., 19(2):447–459, 2006.
  * [24] J. M. Landsberg.  Geometry and complexity theory, volume 169 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, 2017.
  * [25] L.-H. Lim.  Tensors in computations.  Acta Numer., 30:555–764, 2021.
  * [26] C. Moore and S. Mertens.  The nature of computation.  Oxford University Press, Oxford, 2011.
  * [27] S. Rudich.  Complexity theory: from Gödel to Feynman.  In Computational complexity theory, volume 10 of IAS/Park City Math. Ser., pages 5–87. Amer. Math. Soc., Providence, RI, 2004.
  * [28] R. A. Ryan.  Introduction to tensor products of Banach spaces.  Springer Monographs in Mathematics. Springer-Verlag London, Ltd., London, 2002.
  * [29] S. Scardapane, S. Van Vaerenbergh, A. Hussain, and A. Uncini.  Complex-valued neural networks with nonparametric activation functions.  IEEE Trans. Emerg. Topics Comput., 4(2):140–150, 2018.
  * [30] V. Strassen.  Gaussian elimination is not optimal.  Numer. Math., 13:354–356, 1969.
  * [31] V. Strassen.  Vermeidung von Divisionen.  J. Reine Angew. Math., 264:184–202, 1973.
  * [32] V. Strassen.  Relative bilinear complexity and matrix multiplication.  J. Reine Angew. Math., 375/376:406–443, 1987.
  * [33] V. Strassen.  Algebraic complexity theory.  In Handbook of theoretical computer science, Vol. A, pages 633–672. Elsevier, Amsterdam, 1990.
  * [34] C. Trabelsi, O. Bilaniuk, Y. Zhang, D. Serdyuk, S. Subramanian, J. F. Santos, S. Mehri, N. Rostamzadeh, Y. Bengio, and C. J. Pal.  Deep complex networks.  In International Conference on Learning Representations, 2018.
  * [35] S. Winograd.  On multiplication of $2\times 2$ matrices.  Linear Algebra Appl., 4:381–388, 1971.
  * [36] K. Ye and L.-H. Lim.  Fast structured matrix computations: tensor rank and Cohn-Umans method.  Found. Comput. Math., 18(1):45–95, 2018.
  * [37] H. Zhang, M. Gu, X. Jiang, J. Thompson, H. Cai, S. Paesani, R. Santagati, A. Laing, Y. Zhang, M. Yung, et al.  An optical neural chip for implementing complex-valued neural network.  Nat. Commun., 12(1):1–11, 2021.