Randomization of Approximate Bilinear Computation for Matrix Multiplication

TL;DR

This paper introduces a randomization technique for approximate bilinear matrix multiplication formulas, improving performance in the presence of computational and approximation errors through theoretical analysis and experiments.

Contribution

It proposes a novel randomization method for approximate bilinear formulas and analyzes its effectiveness under different error sources.

Findings

Improved accuracy with randomization in approximate formulas

Enhanced performance when combined with floating point arithmetic errors

Theoretical and experimental validation of the method's benefits

Abstract

We present a method for randomizing formulas for bilinear computation of matrix products. We consider the implications of such randomization when there are two sources of error: One due to the formula itself only being approximately correct, and one due to using floating point arithmetic. Our theoretical results and numerical experiments indicate that our method can improve performance when each of these error sources are present individually, as well as when they are present at the same time.

Equations

$$c_{\_}{ij}=\sum_{\_}{r=1}^{R}w_{\_}{ijr}\bigg{(}\sum_{\_}{k,l=1}^{n}u_{\_}{klr}a_{\_}{kl}\bigg{)}\bigg{(}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}v_{\_}{k^{\prime}l^{\prime}r}b_{\_}{k^{\prime}l^{\prime}}\bigg{)},$$

$$c_{\_}ij=\sum _{\_}{k,l=1}^n\sum _{\_}{k^{\prime },l^{\prime }=1}^n{a}_{\_}kl{b}_{\_}{k}^{\prime }{l}^{\prime }\sum _{\_}{r=1}^R{u}_{\_}klr{v}_{\_}{k}^{\prime }{l}^{\prime }r{w}_{\_}ijr.$$

$$\sum _{\_}{r=1}^R{u}_{\_}klr{v}_{\_}{k}^{\prime }{l}^{\prime }r{w}_{\_}ijr={\delta }_{\_}ki{\delta }_{\_}{l}^{\prime }j{\delta }_{\_}l{k}^{\prime }\text{for all }(k,l,k^{\prime },l^{\prime },i,j)\in [n{]}^6,$$

$$\mathbf{A}=\left[\begin{array}{cccc}{\mathbf{A}}_{\_}11& {\mathbf{A}}_{\_}12& \cdots & {\mathbf{A}}_{\_}1n\\ {\mathbf{A}}_{\_}21& {\mathbf{A}}_{\_}22& \cdots & {\mathbf{A}}_{\_}2n\\ ⋮& ⋮& & ⋮\\ {\mathbf{A}}_{\_}n1& {\mathbf{A}}_{\_}n2& \cdots & {\mathbf{A}}_{\_}nn\end{array}\right].$$

$$\mathcal{X}$$

$$\mathcal{Y}$$

$$f(\mathbf{A},\mathbf{B})_{\_}{ij}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\sum_{\_}{r=1}^{R}w_{\_}{ijr}\bigg{(}\sum_{\_}{k,l=1}^{n}u_{\_}{klr}\mathbf{A}_{\_}{kl}\bigg{)}\bigg{(}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}v_{\_}{k^{\prime}l^{\prime}r}\mathbf{B}_{\_}{k^{\prime}l^{\prime}}\bigg{)}\in\mathbb{R}^{m\times m}\;\;\;\;\text{for }i,j\in[n],$$

$$\mathcal{Y}(\epsilon )+\mathcal{E}(\epsilon )=\mathcal{X},$$

$$\bm{\mathscr{Y}}(\varepsilon)\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\big{(}\sum_{\_}{r=1}^{R}u_{\_}{klr}(\varepsilon)v_{\_}{k^{\prime}l^{\prime}r}(\varepsilon)w_{\_}{ijr}(\varepsilon)\big{)}_{\_}{(k,l,k^{\prime},l^{\prime},i,j)\in[n]^{6}}$$

$$\sum _{\_}{i=1}^{d+1}{\alpha }_{\_}i\mathcal{Y}({\epsilon }_{\_}i)=\mathcal{X},$$

$$\hat{f}(\mathbf{A},\mathbf{B})\stackrel{\text{def}}{=}(1-\kappa {)}^{-1}{\mathbf{M}}_{\_}{1}^{\top }f({\mathbf{M}}_{\_}1\mathbf{A}{\mathbf{M}}_{\_}{2}^{\top },{\mathbf{M}}_{\_}2\mathbf{B}{\mathbf{M}}_{\_}{3}^{\top }){\mathbf{M}}_{\_}3,$$

$$\kappa\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\frac{1}{n^{3}}\sum_{\_}{(i,j,l)\in[n]^{3}}\bigg{(}1-\sum_{\_}{r=1}^{R}u_{\_}{ilr}v_{\_}{ljr}w_{\_}{ijr}\bigg{)}$$

$${\hat{\mathbf{C}}}_{\_}ij=\frac{s_{\_}1(i)s_{\_}3(j)}{1-\kappa }\sum _{\_}{k,l=1}^n\sum _{\_}{k^{\prime },l^{\prime }=1}^n{s}_{\_}1(k)s_{\_}2(l)s_{\_}2(k^{\prime })s_{\_}3(l^{\prime }){\mathbf{A}}_{\_}kl{\mathbf{B}}_{\_}{k}^{\prime }{l}^{\prime }\sum _{\_}{r=1}^R{u}_{\_}{\pi }_{\_}1(k){\pi }_{\_}2(l)r{v}_{\_}{\pi }_{\_}2(k^{\prime }){\pi }_{\_}3(l^{\prime })r{w}_{\_}{\pi }_{\_}1(i){\pi }_{\_}3(j)r.$$

$$\mathbb{E}[\hat{\mathbf{C}}_{\_}{ij}]=(1-\kappa)^{-1}\sum_{\_}{l=1}^{n}\mathbf{A}_{\_}{il}\mathbf{B}_{\_}{lj}\mathbb{E}\Big{[}\sum_{\_}{r=1}^{R}u_{\_}{\pi_{\_}1(i)\pi_{\_}2(l)r}v_{\_}{\pi_{\_}2(l)\pi_{\_}3(j)r}w_{\_}{\pi_{\_}1(i)\pi_{\_}3(j)r}\Big{]}=\sum_{\_}{l=1}^{n}\mathbf{A}_{\_}{il}\mathbf{B}_{\_}{lj}=(\mathbf{A}\mathbf{B})_{\_}{ij},$$

$$\mathbb{E}\Big{[}\sum_{\_}{r=1}^{R}u_{\_}{\pi_{\_}1(i)\pi_{\_}2(l)r}v_{\_}{\pi_{\_}2(l)\pi_{\_}3(j)r}w_{\_}{\pi_{\_}1(i)\pi_{\_}3(j)r}\Big{]}=1-\kappa.$$

$$\displaystyle\|f(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|^{2}=\sum_{\_}{i,j=1}^{n}\bigg{\|}\sum_{\_}{k,l=1}^{n}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}\mathbf{A}_{\_}{kl}\mathbf{B}_{\_}{k^{\prime}l^{\prime}}(y_{\_}{klk^{\prime}l^{\prime}ij}-x_{\_}{klk^{\prime}l^{\prime}ij})\bigg{\|}^{2}$$

$$\displaystyle\leq\sum_{\_}{i,j=1}^{n}\bigg{(}\sum_{\_}{k,l=1}^{n}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}\|\mathbf{A}_{\_}{kl}\|\|\mathbf{B}_{\_}{k^{\prime}l^{\prime}}\||y_{\_}{klk^{\prime}l^{\prime}ij}-x_{\_}{klk^{\prime}l^{\prime}ij}|\bigg{)}^{2}$$

$$\displaystyle\leq\sum_{\_}{i,j=1}^{n}\bigg{(}\sum_{\_}{k,l=1}^{n}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}(\|\mathbf{A}_{\_}{kl}\|\|\mathbf{B}_{\_}{k^{\prime}l^{\prime}}\|)^{2}\bigg{)}\bigg{(}\sum_{\_}{k,l=1}^{n}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}|y_{\_}{klk^{\prime}l^{\prime}ij}-x_{\_}{klk^{\prime}l^{\prime}ij}|^{2}\bigg{)}$$

$$\displaystyle=\bigg{(}\sum_{\_}{k,l=1}^{n}\|\mathbf{A}_{\_}{kl}\|^{2}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}\|\mathbf{B}_{\_}{k^{\prime}l^{\prime}}\|^{2}\bigg{)}\bigg{(}\sum_{\_}{i,j=1}^{n}\sum_{\_}{k,l=1}^{n}\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}|y_{\_}{klk^{\prime}l^{\prime}ij}-x_{\_}{klk^{\prime}l^{\prime}ij}|^{2}\bigg{)}$$

$$=\Vert \mathbf{A}{\Vert }^2\Vert \mathbf{B}{\Vert }^2\Vert \mathcal{Y}-\mathcal{X}{\Vert }^2\le {\mu }^4{\tau }^2,$$

$$\Vert \hat{f}(\mathbf{A},\mathbf{B})-\mathbf{AB}\Vert =\Vert (1-\kappa {)}^{-1}f(\tilde{\mathbf{A}},\tilde{\mathbf{B}})-\tilde{\mathbf{A}}\tilde{\mathbf{B}}\Vert ,$$

$$\Vert \hat{f}(\mathbf{A},\mathbf{B})-\mathbf{AB}\Vert \le \Vert \tilde{\mathbf{A}}\Vert \Vert \tilde{\mathbf{B}}\Vert \Vert (1-\kappa {)}^{-1}\mathcal{Y}-\mathcal{X}\Vert =\Vert \mathbf{A}\Vert \Vert \mathbf{B}\Vert \Vert (1-\kappa {)}^{-1}\mathcal{Y}-\mathcal{X}\Vert \le {\mu }^2(|\eta |\Vert \mathcal{Y}\Vert +\tau )$$

$$\Vert \hat{f}(\mathbf{A},\mathbf{B})-\mathbf{AB}\Vert \le \Vert f(\tilde{\mathbf{A}},\tilde{\mathbf{B}})-\tilde{\mathbf{A}}\tilde{\mathbf{B}}\Vert +|\eta |\Vert f(\tilde{\mathbf{A}},\tilde{\mathbf{B}})\Vert .$$

$$\Vert f(\tilde{\mathbf{A}},\tilde{\mathbf{B}})\Vert \le \Vert \tilde{\mathbf{A}}\Vert \Vert \tilde{\mathbf{B}}\Vert \Vert \mathcal{Y}\Vert \le {\mu }^2\Vert \mathcal{Y}\Vert ,$$

$$|\eta |{\mu }^2\Vert \mathcal{Y}\Vert \le 2{\mu }^2(n^{-5/2}\tau +n^{-1})\tau .$$

$$|\kappa|=\bigg{|}n^{-3}\sum_{\_}{i=1}^{n^{3}}z_{\_}i\bigg{|}\leq n^{-3}\|\mathbf{z}\|_{\_}1\leq n^{-3}\sqrt{n}\|\mathbf{z}\|\leq n^{-5/2}\|\bm{\mathscr{Y}}-\bm{\mathscr{X}}\|,$$

$$\Vert \mathcal{Y}\Vert \le \Vert \mathcal{Y}-\mathcal{X}\Vert +\Vert \mathcal{X}\Vert =\Vert \mathcal{Y}-\mathcal{X}\Vert +n^{3/2}.$$

$$({\hat{F}}^{(q)}(\mathbf{A},\mathbf{B}){)}_{\_}ij\stackrel{\text{def}}{=}(1-\kappa {)}^{-1}{s}_{\_}{1}^{(q)}(i)s_{\_}{3}^{(q)}(j)\sum _{\_}{r=1}^R{w}_{\_}{\pi }_{\_}{1}^{(q)}(i){\pi }_{\_}{3}^{(q)}(j)r{\hat{F}}^{(q-1)}\left({\mathbf{A}}^{(q)},{\mathbf{B}}^{(q)}\right),$$

$${\mathbf{A}}^{(q)}\stackrel{\text{def}}{=}\sum _{\_}{k,l=1}^n{u}_{\_}{\pi }_{\_}{1}^{(q)}(k){\pi }_{\_}{2}^{(q)}(l)r{s}_{\_}{1}^{(q)}(k)s_{\_}{2}^{(q)}(l){\mathbf{A}}_{\_}kl,$$

$${\mathbf{B}}^{(q)}\stackrel{\text{def}}{=}\sum _{\_}{k^{\prime },l^{\prime }=1}^n{v}_{\_}{\pi }_{\_}{2}^{(q)}(k^{\prime }){\pi }_{\_}{3}^{(q)}(l^{\prime })r{s}_{\_}{2}^{(q)}(k^{\prime })s_{\_}{3}^{(q)}(l^{\prime }){\mathbf{B}}_{\_}{k}^{\prime }{l}^{\prime },$$

$$\mathbb{E}[({\hat{F}}^{(Q+1)}(\mathbf{A},\mathbf{B}){)}_{\_}ij]=\mathbb{E}\left[{\mathbb{E}}_{\_}\mathcal{S}[({\hat{F}}^{(Q+1)}(\mathbf{A},\mathbf{B}){)}_{\_}ij]\right]$$

$$\displaystyle=\mathbb{E}\bigg{[}(1-\kappa)^{-1}s_{\_}1^{(Q+1)}(i)s_{\_}3^{(Q+1)}(j)\sum_{\_}{r=1}^{R}w_{\_}{\pi_{\_}1^{(Q+1)}(i)\pi_{\_}3^{(Q+1)}(j)r}$$

$$\displaystyle\;\;\;\;\times\mathbb{E}_{\_}{\mathcal{S}}\Big{[}\hat{F}^{(Q)}\Big{(}\sum_{\_}{k,l=1}^{n}u_{\_}{\pi_{\_}1^{(Q+1)}(k)\pi_{\_}2^{(Q+1)}(l)r}s_{\_}1^{(Q+1)}(k)s_{\_}2^{(Q+1)}(l)\mathbf{A}_{\_}{kl},$$

$$\displaystyle\;\;\;\;\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}v_{\_}{\pi_{\_}2^{(Q+1)}(k^{\prime})\pi_{\_}3^{(Q+1)}(l^{\prime})r}s_{\_}2^{(Q+1)}(k^{\prime})s_{\_}3^{(Q+1)}(l^{\prime})\mathbf{B}_{\_}{k^{\prime}l^{\prime}}\Big{)}\Big{]}\bigg{]}$$

Full text

1 Introduction

Suppose $\mathbf{A},\mathbf{B}\in\mathbb{R}^{n\times n}$. In this paper, we are concerned with formulas for computing $\mathbf{C}=\mathbf{A}\mathbf{B}$ that take the form

[TABLE]

$i,j\in[n]\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\{1,2,\ldots,n\}$, where each $\bm{\mathscr{U}}=(u_{\_}{klr})$, $\bm{\mathscr{V}}=(v_{\_}{k^{\prime}l^{\prime}r})$, and $\bm{\mathscr{W}}=(w_{\_}{ijr})$ is a tensor containing real numbers, and $c_{\_}{ij}$ is the element at position $(i,j)$ in $\mathbf{C}$, with similar notation for elements of $\mathbf{A}$ and $\mathbf{B}$. Such a formula is called a bilinear computation (BC) [7]. We can rewrite the expression above as

[TABLE]

Comparing this to the standard algorithm for matrix multiplication, we can see that for the computation (1) to be exact, $\bm{\mathscr{U}}$, $\bm{\mathscr{V}}$ and $\bm{\mathscr{W}}$ must satisfy

[TABLE]

where $\delta$ is the Kronecker delta (i.e., $\delta_{\_}{ki}=1$ if $k=i$ and [math] otherwise) [8]. If (3) is satisfied, we say that (1) is an exact bilinear computation (EBC). The smallest positive integer $R$ for which there exist $\bm{\mathscr{U}}$, $\bm{\mathscr{V}}$ and $\bm{\mathscr{W}}$ such that (3) holds is referred to as the rank of the computation. We refer to any algorithm of the form (1) for which $R<n^{3}$ as fast, since the asymptotic complexity is smaller than that of standard matrix multiplication, which has complexity $O(n^{3})$. Examples of fast computations of the form (1) include Strassen’s algorithm for $2\times 2$ matrices [30], and Laderman’s algorithm for $3\times 3$ matrices [25]. The formula (1) is also valid if $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ are of size $mn\times mn$ and $a_{\_}{kl}$, $b_{\_}{k^{\prime}l^{\prime}}$ and $c_{\_}{ij}$ are replaced by submatrices $\mathbf{A}_{\_}{kl}$, $\mathbf{B}_{\_}{k^{\prime}l^{\prime}}$ and $\mathbf{C}_{\_}{ij}$ of size $m\times m$, so that e.g.

[TABLE]

This is why fast algorithms of the form (1) can be used recursively to compute the product of larger matrices.

Define 6-way tensors

[TABLE]

The condition in (3) can be written succinctly as $\bm{\mathscr{Y}}=\bm{\mathscr{X}}$. For both matrix and tensor inputs, let $\|\cdot\|$ denote the Frobenius norm, i.e., the square root of the sum of the square of all elements. We are interested in BCs that are only approximately correct, which motivates the following definition.

Definition 1** (Approximate bilinear computation).**

Let $\mathbf{A},\mathbf{B}\in\mathbb{R}^{mn\times mn}$, where $m$ is a positive integer, and let $f$ be defined blockwise via

[TABLE]

where $f(\mathbf{A},\mathbf{B})_{\_}{ij}$, $\mathbf{A}_{\_}{kl}$ and $\mathbf{B}_{\_}{k^{\prime}l^{\prime}}$ are submatrices of $f(\mathbf{A},\mathbf{B})$, $\mathbf{A}$ and $\mathbf{B}$, respectively, of size $m\times m$. We say that $f$ is an approximate bilinear computation (ABC) with parameters $n$ and $\tau$, or $(n,\tau)$-ABC for short, if $\|\bm{\mathscr{Y}}-\bm{\mathscr{X}}\|\leq\tau$.

We present a method for randomizing the computation (7) and consider the implications of this approach when there are two kinds of error present.

• •

Error due to an approximate algorithm. Using an $(n,\tau)$-ABC with $\tau>0$ will introduce error even if exact arithmetic is used. ABCs are interesting since it is sometimes hard or impossible to convert an approximate algorithm found numerically to an exact algorithm. Some examples of such approximate algorithms found numerically appear in [29]. Other papers that search for fast algorithms numerically include [8, 23, 2, 17].

• •

Numerical error due to using floating point arithmetic. Although the standard algorithm for matrix multiplication also incurs numerical error, it is more severe in fast algorithms since the computation of an element $c_{\_}{ij}$ can involve elements from $\mathbf{A}$ and $\mathbf{B}$ other than the vectors $\mathbf{a}_{\_}{i:}$ and $\mathbf{b}_{\_}{:j}$, which are the $i$th row of $\mathbf{A}$ and $j$th column of $\mathbf{B}$, respectively. In exact arithmetic, these additional terms cancel out for exact fast algorithms, but in finite precision those cancellations are typically not exact which can lead to substantial numerical error. These issues are exacerbated when the algorithm itself is only approximately correct. These considerations are especially important when using low precision environments, such as when computing on a GPU [21]. Low precision computation—using 32-bit, 16-bit, and even 8-bit precision numbers—is popular in, e.g., machine learning [19, 31, 20].

We make the following contributions in this paper:

• •

We propose a method for randomizing BCs for matrix multiplication which does not increase the leading order computational complexity of the algorithm.

• •

When exact arithmetic is used, we show that our randomized ABCs compute the correct matrix product in expectation. We also provide some performance guarantees.

• •

We show that these exact arithmetic results largely carry over to a setting when all computations are done in floating point arithmetic.

• •

When floating point arithmetic is used, we provide numerical evidence that randomizing EBCs using our scheme can reduce numerical error and improve robustness to adversarial examples.

1.1 Related work

Bini et al. [4, 6] introduce a concept similar to ABC called Arbitrary Precision Approximating (APA) algorithms for matrix multiplication. For an APA algorithm, the tensor $\bm{\mathscr{X}}$ and $\bm{\mathscr{Y}}$ defined in (5) and (6) satisfy the relationship

[TABLE]

where

[TABLE]

is a function of $\varepsilon$, and $\bm{\mathscr{E}}$ represents the error in the approximation $\bm{\mathscr{Y}}\approx\bm{\mathscr{X}}$. Each entry in $\bm{\mathscr{E}}$ is assumed to be a polynomial with zero constant coefficient, hence for every entry $\bm{\mathscr{E}}_{\_}{klk^{\prime}l^{\prime}ij}(0)=0$. Consequently, an APA algorithm can be made arbitrarily accurate by making $\varepsilon$ small enough. Moreover, as shown in [3], an EBC can be derived from an APA algorithm (see also the related discussion in Section 15.2 of [10]). This EBC takes the form

[TABLE]

where all $\varepsilon_{\_}i$ are distinct, $d:=\max_{\_}{klk^{\prime}l^{\prime}ij\in[n]^{6}}\deg(\bm{\mathscr{E}}_{\_}{klk^{\prime}l^{\prime}ij}(\varepsilon))$ is the maximum degree of the polynomials describing the entries in $\bm{\mathscr{E}}$, and where $[\alpha_{\_}1,\,\ldots,\,\alpha_{\_}{d+1}]$ is chosen as the solution to a certain linear system. Bini [3] uses this approach to convert the APA scheme for $12\times 12$ matrix multiplication in [4] to an exact one, which can be used recursively to get a matrix multiplication algorithm with complexity $O(n^{2.7799})$. Fixing the error in an APA algorithm yields an ABC with some parameter $\tau$. In some of our numerical experiments, we use an ABC which we get by fixing $\varepsilon$ in the $12\times 12$ APA scheme of [4]. We also consider an associated EBC derived via (10) in other experiments. Our paper focuses on ABCs instead of APA algorithms since, in practice, it may be difficult or impossible to express a pair $(\bm{\mathscr{Y}},\bm{\mathscr{E}})$ found numerically in terms of polynomials in $\varepsilon$.

To the best of our knowledge, our paper is the first to consider randomization as a tool for improving BCs for matrix multiplication which are only approximately correct. Various randomized algorithms for the standard matrix multiplication algorithm have been considered in other works; see e.g. [15, 26]. For EBCs in floating point arithmetic, a patent by Castrapel and Gustafson [11] describes a randomized version of Strassen’s algorithm which they claim reduces numerical error. They provide empirical support for this, but no mathematical proof. Our method generalizes their approach by randomly choosing from a wider range of equivalent algorithms. Additionally, our method can be applied to any formula of the form (1).

Early works that analyze the stability of fast algorithms for matrix multiplication include [9, 6, 5]. Other works, such as [16, 24, 14, 1, 12], attempt to improve the stability of fast algorithms using other approaches that do not rely on randomization. In Section 3.2, we compare our proposed method to the rescaling method in [1]. In these experiments, we also consider a restricted version of our method which corresponds to the method in [11].

2 Randomization of bilinear computation for matrix multiplication

In Section 2.1, we first consider a setting in which exact arithmetic is used. In Section 2.2, we then consider a setting in which floating point arithmetic is used.

2.1 In exact arithmetic

We present our randomization scheme in the setting when the input matrices are $mn\times mn$ with $m\geq 1$111We present results for square matrices, but we believe the results can be extended to rectangular matrices.. Accordingly, suppose $f:\mathbb{R}^{mn\times mn}\times\mathbb{R}^{mn\times mn}\rightarrow\mathbb{R}^{mn\times mn}$ is an $(n,\tau)$-ABC. We will now define a randomized version of $f$, denoted by $\hat{f}$, which has the following property: For all $\mathbf{A},\mathbf{B}\in\mathbb{R}^{mn\times mn}$, $\mathbb{E}[\hat{f}(\mathbf{A},\mathbf{B})]=\mathbf{A}\mathbf{B}$. To that end, let $\{s_{\_}i(j)\}_{\_}{(i,j)\in[3]\times[n]}$ be a collection of i.i.d. Rademacher random variables, i.e., each satisfying $\mathbb{P}[s_{\_}i(j)=+1]=\mathbb{P}[s_{\_}i(j)=-1]=1/2$. Moreover, let $\pi_{\_}i:[n]\rightarrow[n]$, $i\in[3]$, be independent random permutation functions, each satisfying $(\forall(j,k)\in[n]^{2})$ $\mathbb{P}[\pi_{\_}i(j)=k]=1/n$. Let $\mathbf{S}_{\_}i\in\mathbb{R}^{mn\times mn}$, $i\in[3]$, be block diagonal matrices with the $j$th nonzero block equal to $s_{\_}i(j)\mathbf{I}_{\_}m$, where $\mathbf{I}_{\_}m$ is the $m\times m$ identity matrix. Also, let $\mathbf{P}_{\_}i\in\mathbb{R}^{mn\times mn}$, $i\in[3]$, be permutation matrices divided into $m\times m$ blocks, with blocks on position $(\pi_{\_}i(j),j)$, $j\in[n]$, equal to $\mathbf{I}_{\_}m$ and all other blocks equal to zero. Define $\mathbf{M}_{\_}i=\mathbf{P}_{\_}i\mathbf{S}_{\_}i$, $i\in[3]$. Note that each $\mathbf{M}_{\_}i$ is orthogonal (i.e., $\mathbf{M}_{\_}i^{-1}=\mathbf{M}_{\_}i^{\top}$). We propose the following definition of $\hat{f}$.

Definition 2** (Randomized approximate bilinear computation).**

Let $f$ be an $(n,\tau)$-ABC. We define a corresponding randomized approximate bilinear computation with parameters $n$, $\tau$ and $\kappa$, or $(n,\tau,\kappa)$-RandABC for short, via

[TABLE]

where

[TABLE]

is assumed to satisfy $\kappa\neq 1$.

Observe that if $f$ was an exact algorithm, then $\kappa=0$ and $\hat{f}(\mathbf{A},\mathbf{B})=\mathbf{A}\mathbf{B}$ since the $\mathbf{M}_{\_}i$’s would cancel out due to orthogonality. Since the cost of applying $\mathbf{M}_{\_}i$ to an $mn\times mn$ matrix is $O(m^{2}n^{2})$, computing $\hat{f}(\mathbf{A},\mathbf{B})$ has the same leading order complexity as computing $f(\mathbf{A},\mathbf{B})$.

Proposition 3**.**

Let $\hat{f}$ be an $(n,\tau,\kappa)$-RandABC with $\tau\geq 0$ and $\kappa\neq 1$ . For all $\mathbf{A},\mathbf{B}\in\mathbb{R}^{mn\times mn}$ we have $\mathbb{E}[\hat{f}(\mathbf{A},\mathbf{B})]=\mathbf{A}\mathbf{B}$.

Proof.

Let $\hat{\mathbf{C}}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\hat{f}(\mathbf{A},\mathbf{B})$. Considering the $(i,j)$th block of $\hat{\mathbf{C}}$, and going through some tedious but straightforward algebra, we get

[TABLE]

If we take the expectation of this equation with respect to the random variables $\{s_{\_}i(j)\}$, most terms will vanish: If $i=k$, $l=k^{\prime}$ and $j=l^{\prime}$ for a given term, then the product of the $s_{\_}i$’s will be 1; otherwise, the expectation of that term will be zero due to independence and the fact that each $\mathbb{E}[s_{\_}i(j)]=0$. Consequently, we have

[TABLE]

where the second equality is true since (12) implies that

[TABLE]

∎

It may seem surprising that Proposition 3 holds for any $\tau\geq 0$, since this means that even an ABC with an arbitrarily large error will be correct in expectation once randomized as in Definition 2. However, as we will see in Proposition 4, in order to be able to guarantee a small error for any realization of $\hat{f}(\mathbf{A},\mathbf{B})$, $\tau$ also needs to be small.

Define $B_{\_}\mu^{(n)}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\left\{\mathbf{M}\in\mathbb{R}^{n\times n}:\|\mathbf{M}\|\leq\mu\right\}$ and $\eta\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}(1-\kappa)^{-1}-1=O(\kappa)$. The following proposition provides performance guarantees for $f$ and $\hat{f}$.

Proposition 4**.**

Consider matrices $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}$.

1. (i)

If $f$ is an $(n,\tau)$-ABC, then $\|f(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|\leq\mu^{2}\tau$. 2. (ii)

If $\hat{f}$ is an $(n,\tau,\kappa)$-RandABC, then $\|\hat{f}(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|\leq\mu^{2}|\eta|\|\bm{\mathscr{Y}}\|+\mu^{2}\tau\lesssim\mu^{2}|\kappa|\|\bm{\mathscr{Y}}\|+\mu^{2}\tau$. 3. (iii)

Moreover, $\sup_{\_}{\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}}\|\hat{f}(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|\leq\sup_{\_}{\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}}\|f(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|+\mu^{2}|\eta|\|\bm{\mathscr{Y}}\|$.

Proof.

We have

[TABLE]

where the first inequality follows from first applying the triangle inequality and then using the sub-multiplicativity of the Frobenius norm, and the second inequality follows from applying the Cauchy–Schwarz inequality. This proves (i). Since the Frobenius norm is invariant under unitary transformations,

[TABLE]

where $\tilde{\mathbf{A}}=\mathbf{M}_{\_}1\mathbf{A}\mathbf{M}_{\_}2^{\top}$, $\tilde{\mathbf{B}}=\mathbf{M}_{\_}2\mathbf{B}\mathbf{M}_{\_}3^{\top}$. Going through the same computations as in the proof of (i), but with the extra $(1-\kappa)^{-1}$ term, we therefore get

[TABLE]

where the equality once again uses the unitary invariance of the Frobenius norm, and the last inequality follows from the definitions of $\eta$ and the triangle inequality. This proves (ii). Fix $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}$. Applying the triangle inequality to (21) and using the definition of $\eta$, we get

[TABLE]

By doing computations almost identical to those in the proof of (i), we get the bound

[TABLE]

since $\tilde{\mathbf{A}},\tilde{\mathbf{B}}\in B_{\_}\mu^{(mn)}$. Combining (23) and (24) and taking supremums appropriately proves (iii). ∎

Points (i) and (ii) in Proposition 4 provide absolute performance guarantees for $f$ and $\hat{f}$, respectively. Note that a tighter version of (i) holds, with $\|\mathbf{A}\|$ replaced by $\sqrt{\sum_{\_}{kl}\|\mathbf{A}_{\_}{kl}\|_{\_}2^{2}}$, where $\|\cdot\|_{\_}2$ is the spectral norm. We keep (i) in its current looser form to make it easier to compare to (ii) and (iii). Point (iii) shows that the worst case performance of $\hat{f}$ is no worse than that of $f$ plus a constant. In fact, that constant, which also appears in (ii), can be bounded as follows.

Proposition 5**.**

If $|\kappa|\leq 1/2$, then

[TABLE]

Proof.

Let $\mathbf{z}\in\mathbb{R}^{n^{3}}$ be the vector with elements $(\sum_{\_}{r=1}^{R}u_{\_}{ilr}v_{\_}{ljr}w_{\_}{ijr}-1)$ for $(i,j,l)\in[n]^{3}$, i.e., containing the elements of $(\bm{\mathscr{Y}}-\bm{\mathscr{X}})$ in positions for which $\bm{\mathscr{X}}$ has an entry 1 (the element order in $\mathbf{z}$ is irrelevant). Note that

[TABLE]

where the first equality follows from (12), the first inequality follows from the triangle inequality and the definition of the 1-norm, and the second inequality is a well known relation (see e.g. Equation (2.2.5) in [18]). Now, note that

[TABLE]

Combining (26), (27), the fact that $|\eta|\leq 2|\kappa|$ when $|\kappa|\leq 1/2$, and the definition of $\tau$ gives us the desired bound. ∎

The upper bound $\mu^{2}\tau$ in Proposition 4 (i) is also an upper bound to the $\sup$ term on the right hand side of the inequality in (iii) of the same proposition. Proposition 5 shows that the size of the additional constant $|\eta|\mu^{2}\|\bm{\mathscr{Y}}\|$ is not much larger than the bound that we already have on this $\sup$ term, and that it will be smaller than that bound if e.g. $n\geq 3$ and $\tau<3^{3/2}/2$.

As is clear from the proof of Proposition 3, the constant $\kappa$ in (12) is used to rescale the computation in (11) so that the output of $\hat{f}$ is correct in expectation. It is important to note that $\kappa$, which can be both positive and negative, is not a measure of error in the algorithm. Indeed, setting $R=1$ and all elements of $\bm{\mathscr{U}}$, $\bm{\mathscr{V}}$ and $\bm{\mathscr{W}}$ to 1 would result in $\kappa=0$, but this would clearly be a very poor algorithm. Although the corresponding randomized algorithm $\hat{f}$ would be correct in expectation, the error guarantees in Proposition 4 (i) and (ii) would be very poor, since $\tau$ would be large.

2.1.1 A recursive algorithm for approximate bilinear computation

We can extend the result in Proposition 3 to a recursive version of $\hat{f}$. We denote the recursive algorithm with $Q$ recursions for multiplication of $mn^{Q}\times mn^{Q}$ matrices by $\hat{F}^{(Q)}$. Let $\{s_{\_}i^{(q)}(j)\}_{\_}{(i,j,q)\in[3]\times[n]\times[Q]}$ be a collection of i.i.d. Rademacher random variables, and let $\pi_{\_}i^{(q)}:[n]\rightarrow[n]$, $(i,q)\in[3]\times[Q]$, be independent random permutation functions, each satisfying $(\forall(j,k)\in[n]^{2})$ $\mathbb{P}[\pi_{\_}i^{(q)}(j)=k]=1/n$. Let $\hat{F}^{(1)}$ be defined exactly as $\hat{f}$ in (11) but based on the random variables $\{s_{\_}i^{(1)}(j)\}_{\_}{(i,j)\in[3]\times[n]}$ and $\{\pi_{\_}i^{(1)}\}_{\_}{i\in[3]}$ and define $\hat{F}^{(q)}:\mathbb{R}^{mn^{q}\times mn^{q}}\times\mathbb{R}^{mn^{q}\times mn^{q}}\rightarrow\mathbb{R}^{mn^{q}\times mn^{q}}$, $q\in\{2,3,\ldots,Q\}$, recursively via

[TABLE]

where

[TABLE]

and $(\hat{F}^{(q)}(\mathbf{A},\mathbf{B}))_{\_}{ij}$ is the subblock of size $mn^{q-1}\times mn^{q-1}$ on position $(i,j)$. If $\mathbf{A}$ and $\mathbf{B}$ are of size $p\times p$ but there is no integer $m$ such that $p=mn^{Q}$, e.g. if $p$ is prime, then we can simply pad the matrices appropriately [22]. If the recursive formula (28) is used $Q$ times to compute the product of two $N\times N$ matrices, where $N\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}mn^{Q}$, and $\hat{F}^{(0)}(\mathbf{A}^{(0)},\mathbf{B}^{(0)})=\mathbf{A}^{(0)}\mathbf{B}^{(0)}$ is computed via the standard matrix multiplication formula, the asymptotic cost of the recursive algorithm is $O(m^{3-\log_{\_}nR}N^{\log_{\_}nR})$. This is the same leading order complexity as the corresponding computation without any randomization. For example, if we insert $m=1$, $n=2$ and $R=7$, we recover the asymptotic cost of Strassen’s algorithm: $O(N^{\log_{\_}27})\approx O(N^{2.81})$.

Proposition 6**.**

For any positive integer $Q$ and for all $\mathbf{A},\mathbf{B}\in\mathbb{R}^{mn^{Q}\times mn^{Q}}$ we have $\mathbb{E}[\hat{F}^{(Q)}(\mathbf{A},\mathbf{B})]=\mathbf{A}\mathbf{B}$.

Proof.

Note that the claim is true for $Q=1$ due to Proposition 3. Now assume it is true for some $Q\geq 1$. We will show that it is also true for $Q+1$, i.e., that $\mathbb{E}[\hat{F}^{(Q+1)}(\mathbf{A},\mathbf{B})]=\mathbf{A}\mathbf{B}$ for $\mathbf{A},\mathbf{B}\in\mathbb{R}^{mn^{Q+1}\times mn^{Q+1}}$. Let $\mathcal{S}$ denote the $\sigma$-algebra generated by the random variables $\{s_{\_}i^{(Q+1)}(j)\}_{\_}{(i,j)\in[3]\times[n]}$ and $\{\pi_{\_}i^{(Q+1)}(j)\}_{\_}{(i,j)\in[3]\times[n]}$, and let $\mathbb{E}_{\_}{\mathcal{S}}[\;\cdot\;]\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\mathbb{E}[\;\cdot\mid\mathcal{S}]$. Then

[TABLE]

due to the smoothing property of expectation (property 10 in [27, p. 348]),

[TABLE]

since each $s_{\_}i^{(Q+1)}(j)$ and $\pi_{\_}i^{(Q+1)}(j)$ is $\mathcal{S}$-measurable,

[TABLE]

due to the induction hypothesis and since all random variables $\{s_{\_}i^{(q)}(j)\}_{\_}{(i,j,q)\in[3]\times[n]\times[Q]}$ and $\{\pi_{\_}i^{(q)}(j)\}_{\_}{(i,j,q)\in[3]\times[n]\times[Q]}$ are independent of $\mathcal{S}$ (and using property 12 in [27, pp. 349–350]),

[TABLE]

by reordering the terms,

[TABLE]

which follows by doing the same analysis as in the proof of Proposition 3. So $\mathbb{E}[\hat{F}^{(Q+1)}(\mathbf{A},\mathbf{B})]=\mathbf{A}\mathbf{B}$. The claim in Proposition 6 now follows by induction. ∎

2.2 In floating point arithmetic

In floating point arithmetic, we cannot achieve guarantees like those in Propositions 3 and 6. This is illustrated in the following example.

Example 7**.**

Let $g:\mathbb{R}^{mn\times mn}\times\mathbb{R}^{mn\times mn}\rightarrow\mathbb{R}^{mn\times mn}$ be a function that computes matrix multiplication according to some EBC in floating point arithmetic. Let $\hat{g}$ be defined analogously to $\hat{f}$ in (11), but in terms of $g$ instead of $f$. Since $g$ is exact, we have $\kappa=0$. We will use $\hat{G}^{(Q)}$ to denote the recursive version of $\hat{g}$, defined analogously to $\hat{F}^{(Q)}$. We will use $G^{(Q)}$ to denote the recursive version of $g$, defined analogously to $\hat{G}^{(Q)}$ but with each $s_{\_}i^{(q)}(j)=1$ and each $\pi_{\_}i^{(q)}(j)=j$, i.e., with no randomness involved so that $G^{(Q)}$ is deterministic. We consider the same setup as in Section 2.7.10 of [18]: Let

[TABLE]

and suppose we are computing on a machine using 2-digit floating point arithmetic. Taking $g$ to be Strassen’s algorithm computed on this machine, we get $\|g(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|\approx 0.0286$. In the definition of $\hat{g}$, there are a total of 64 possible sign functions and 8 possible permutation functions. Each combination of these has the same probability of occurring, so we can readily compute $\mathbb{E}[\hat{g}(\mathbf{A},\mathbf{B})]$. Doing this, we find that $\|\mathbb{E}[\hat{g}(\mathbf{A},\mathbf{B})]-\mathbf{A}\mathbf{B}\|\approx 0.0024$. So we cannot guarantee $\mathbb{E}[\hat{g}(\mathbf{A},\mathbf{B})]=\mathbf{A}\mathbf{B}$ in this finite precision setting. We are not even guaranteed to have $\|\mathbb{E}[\hat{g}(\mathbf{A},\mathbf{B})]-\mathbf{A}\mathbf{B}\|\leq\|g(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|$: Let $\mathbf{A}^{(r)}$ and $\mathbf{B}^{(r)}$ be the same as $\mathbf{A}$ and $\mathbf{B}$, respectively, but with the order of the columns reversed, i.e., $\mathbf{A}^{(r)}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\begin{bmatrix}\mathbf{a}_{\_}{:2}&\mathbf{a}_{\_}{:1}\end{bmatrix}$ and $\mathbf{B}^{(r)}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\begin{bmatrix}\mathbf{b}_{\_}{:2}&\mathbf{b}_{\_}{:1}\end{bmatrix}$. We then get $\|g(\mathbf{A}^{(r)},\mathbf{B}^{(r)})-\mathbf{A}^{(r)}\mathbf{B}^{(r)}\|\approx 0.0001$, but $\|\mathbb{E}[\hat{g}(\mathbf{A}^{(r)},\mathbf{B}^{(r)})]-\mathbf{A}^{(r)}\mathbf{B}^{(r)}\|\approx 0.0024$. Despite this, randomization seems to work remarkably well in practice when an exact algorithm is computed in finite precision arithmetic, as we will see in the numerical experiments.

We now consider ABCs (which include EBCs as a special case) in finite precision arithmetic. As in [18], we use $\operatorname{fl}(x)$ to denote the representation of $x\in\mathbb{R}$ as a floating point number, and $\operatorname{fl}(f(x))$ to denote the result of computing $f(x)$ in floating point arithmetic. When computing the latter, the algorithm used to compute $f$ matters. We make the standard assumption that $\operatorname{fl}(x\operatorname{op}y)=(x\operatorname{op}y)(1+\Delta)$ where $x$ and $y$ are floating point numbers, $\operatorname{op}$ is scalar addition, subtraction or multiplication, and $|\Delta|\leq\varepsilon_{\_}\text{machine}$, where $\varepsilon_{\_}\text{machine}$ is the machine epsilon or unit roundoff. For numerical summations, e.g. $\operatorname{fl}(\sum_{\_}{n=1}^{N}x_{\_}n)$, we assume that the summation is simply done sequentially, e.g.

[TABLE]

We first present a series of results in Propositions 8–12 with the goal of understanding how algorithmic and numerical error combined impact the deterministic and randomized ABCs described by $f$ and $\hat{f}$, respectively. Throughout these results, we make the reasonable assumption that the input matrices $\mathbf{A}$ and $\mathbf{B}$, as well as the tensors $\bm{\mathscr{U}}$, $\bm{\mathscr{V}}$ and $\bm{\mathscr{W}}$ defining the BC, are already stored in floating point format, so that e.g. $\operatorname{fl}(\mathbf{A})=\mathbf{A}$. Proposition 8 provides an upper bound on the numerical error for the (approximate or exact) BC $f$ defined as in (1).

Proposition 8**.**

Suppose $(4n+R-1)\varepsilon_{\_}\textup{machine}\leq 0.01$. For $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(n)}$ and an $(n,\tau)$-ABC $f$ computed according to (7), we have

[TABLE]

where e.g. $\mathbf{U}_{\_}{::r}\in\mathbb{R}^{n\times n}$ is the $r$th frontal slice of $\bm{\mathscr{U}}$, so that $\|\mathbf{U}_{\_}{::r}\|^{2}=\sum_{\_}{k,l}u_{\_}{klr}^{2}$.

Proof.

Throughout this proof, constants $\gamma$, $\theta$, $\lambda$, $\tilde{\gamma}$, $\tilde{\theta}$, $\tilde{\lambda}$, $\phi$ and $\psi$ with subscripts are real numbers of magnitude less than or equal to $\varepsilon_{\_}\text{machine}$. Define

[TABLE]

Then

[TABLE]

etc. More generally,

[TABLE]

Next, define

[TABLE]

Then

[TABLE]

etc. More generally,

[TABLE]

Consequently,

[TABLE]

Similarly,

[TABLE]

Now, define

[TABLE]

We have

[TABLE]

etc. From this, it follows that

[TABLE]

Writing this out in full and rearranging, we have

[TABLE]

Since we have assumed $(4n+R-1)\varepsilon_{\_}\text{machine}\leq 0.01$, it follows from Lemma 2.7.1 in [18] that there exist constants $\varepsilon_{\_}{klk^{\prime}l^{\prime}r}$, for $(k,l,k^{\prime},l^{\prime},r)\in[n]^{4}\times[R]$, such that $|\varepsilon_{\_}{klk^{\prime}l^{\prime}r}|\leq 1.01(4n+R-1)\varepsilon_{\_}\text{machine}$ and

[TABLE]

We have

[TABLE]

where the first and second inequality follows from the Cauchy–Schwarz inequality, and the third inequality follows from the bound on $|\varepsilon_{\_}{klk^{\prime}l^{\prime}r}|$. Summing over all $i,j$, taking a square root, and using that $\|\mathbf{A}\|,\|\mathbf{B}\|\leq\mu$, we arrive at

[TABLE]

∎

Proposition 9 generalizes Proposition 8 to the case when the input matrices are of size $mn\times mn$.

Proposition 9**.**

Suppose $(4n+m+R-2)\varepsilon_{\_}\textup{machine}\leq 1.01$. For $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}$ and an $(n,\tau)$-ABC $f$ computed according to (7), we have

[TABLE]

Proof.

Note that all computations when calculating $\sum_{\_}{k,l}u_{\_}{klr}\mathbf{A}_{\_}{kl}$ and $\sum_{\_}{k^{\prime},l^{\prime}}v_{\_}{k^{\prime}l^{\prime}r}\mathbf{B}_{\_}{k^{\prime}l^{\prime}}$ are done elementwise. The error accumulated for each entry will therefore be the same as in the case when $\mathbf{A}_{\_}{kl}$ and $\mathbf{B}_{\_}{k^{\prime}l^{\prime}}$ are scalars. More precisely,

[TABLE]

where $\circledast$ is the Hadamard (elementwise) product, and $\mathbf{E}_{\_}{klr}\in\mathbb{R}^{m\times m}$ is a matrix with entries which are the product of at most $2n-1$ terms of the form $(1+\varepsilon)$ with $|\varepsilon|\leq\varepsilon_{\_}\text{machine}$. A similar statement is true for $\operatorname{fl}(\sum_{\_}{k^{\prime},l^{\prime}=1}^{n}v_{\_}{k^{\prime}l^{\prime}r}\mathbf{B}_{\_}{k^{\prime}l^{\prime}})$. When computing each

[TABLE]

we now have additional error compared to the scalar case due to the inner product computations. In order to avoid cumbersome notation, we do not write out the following computations in full. Using the result in Section 2.7.6 of [18] for the rounding error in inner products, Equation (54), Lemma 2.7.1 in [18], and letting $(f(\mathbf{A},\mathbf{B})_{\_}{ij})_{\_}{\alpha\beta}$ denote the element on position $(\alpha,\beta)$ in matrix block $(i,j)$, we can compute

[TABLE]

where each $|\Theta_{\_}{klk^{\prime}l^{\prime}r\alpha\beta z}|\leq\hat{\theta}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}1.01(4n+m+R-2)\varepsilon_{\_}\text{machine}$. Rearranging this, we have

[TABLE]

due to the triangle inequality and definition of $\hat{\theta}$,

[TABLE]

where each inequality follows from an application of the Cauchy–Schwarz inequality and the final equality follows from a rearrangement of terms. Finally, since

[TABLE]

and $\|\mathbf{A}\|,\|\mathbf{B}\|\leq\mu$, it follows that

[TABLE]

∎

Proposition 10 shows how the error in Proposition 9 changes when a randomized algorithm $\hat{f}$ is used instead of its deterministic counterpart. In particular, the bound is worse by only a factor $(1-\kappa)^{-1}$.

Proposition 10**.**

Suppose $(4n+m+R-2)\varepsilon_{\_}\textup{machine}\leq 1.01$ and $1-\kappa>0$. For $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}$ and an $(n,\tau,\kappa)$-RandABC $\hat{f}$ computed according to (11), we have

[TABLE]

Proof.

From (13), it follows that

[TABLE]

where

[TABLE]

Note that $\|\hat{\mathbf{U}}_{\_}{::r}\|=\|\mathbf{U}_{\_}{::r}\|$, $\|\hat{\mathbf{V}}_{\_}{::r}\|=\|\mathbf{V}_{\_}{::r}\|$ and $\|\hat{\mathbf{W}}_{\_}{::r}\|=(1-\kappa)^{-1}\|\mathbf{W}_{\_}{::r}\|$. Using this fact, and applying Proposition 9 to the computation in (66), gives us the desired result. ∎

Propositions 11 and 12 give upper bounds for the total error, both due to algorithmic error and numerical rounding, for the deterministic and randomized algorithms, respectively.

Proposition 11**.**

Suppose $(4n+m+R-2)\varepsilon_{\_}\textup{machine}\leq 1.01$. For $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}$ and an $(n,\tau)$-ABC $f$ computed according to (7), we have

[TABLE]

Proof.

Follows by applying the triangle inequality and using Propositions 9 and 4 (i). ∎

Proposition 12**.**

Suppose $(4n+m+R-2)\varepsilon_{\_}\textup{machine}\leq 1.01$. For $\mathbf{A},\mathbf{B}\in B_{\_}\mu^{(mn)}$ and an $(n,\tau,\kappa)$-RandABC $\hat{f}$ computed according to (11), we have

[TABLE]

Proof.

Follows by taking expectations of the inequality in Proposition 10 and applying Jensen’s inequality [28]. ∎

Note that Proposition 12 implies that

[TABLE]

which means that as the numerical precision increases, $\mathbb{E}[\operatorname{fl}(\hat{f}(\mathbf{A},\mathbf{B}))]$ approaches $\mathbf{A}\mathbf{B}$. This is consistent with Proposition 3.

To make the bounds in (68) and (69) easier to compare, using the fact that $(1-\kappa)^{-1}=1+\kappa+O(\kappa^{2})$, we can rewrite (69) as

[TABLE]

The first term on the right hand side of (68) and (71) are identical. Recall from (26) that $|\kappa|\leq n^{-5/2}\|\bm{\mathscr{Y}}-\bm{\mathscr{X}}\|\leq n^{-5/2}\tau$. Consequently, if the quantity

[TABLE]

is not too large, and $m$, $n$ and $R$ are of moderate size, the second term in (71) will be smaller than the second term in (68), showing that the result in Proposition 3 largely carries over to a setting with floating point arithmetic. The following example illustrates this.

Example 13**.**

For the sake of this example, suppose we are using an approximate variant of Strassen’s algorithm. Then $R=7$ and $n=2$. For Strassen’s algorithm, the quantity in (72) is approximately 6, so we will assume that it is less than 10 for our approximate variant of the algorithm. Suppose we are multiplying two $\text{100,000}\times\text{100,000}$ matrices in single precision, so that $m=\text{50,000}$ and $\varepsilon_{\_}\text{machine}\sim 10^{-8}$. Then, the second term in (71) is upper bounded by a quantity which is on the order of

[TABLE]

which is much smaller than the second term in (68). The implementation of Strassen’s algorithm in [21] achieves a speed-up over an efficient implementation of the standard matrix multiplication algorithm for square matrices with as few as 1,536 rows/columns. So multiplication of matrices with 100,000 rows/columns is well beyond the problem size for which fast algorithms can outperform the standard algorithm.

3 Experiments

In this section we present some results from experiments, with additional results provided in Appendix A. We implement all experiments in Matlab with certain parts implemented in C. All our code is available online at https://github.com/OsmanMalik/random-approximate-matrix-multiplication.

In our experiments, we draw the matrices $\mathbf{A},\mathbf{B}$ from different random distributions. By Gaussian matrix, we mean a matrix whose elements are realizations of i.i.d. standard normal random variables. Similarly, a uniform matrix is one whose elements are realizations of i.i.d. $\operatorname{Uniform}(0,1)$ random variables. We also consider three types of random adversarial matrices, which were proposed by [1] and are designed to be challenging for Strassen’s algorithm.

Definition 14** (Adversarial matrices).**

Consider a matrix pair $\mathbf{A},\mathbf{B}\in\mathbb{R}^{n\times n}$. We say that it is type 1 adversarial if

[TABLE]

We say that the matrix pair is type 2 adversarial if

[TABLE]

We say that the matrix pair is type 3 adversarial if

[TABLE]

Here, all the entries are assumed to be independent.

We will also consider the Hilbert matrix $\mathbf{H}\in\mathbb{R}^{n\times n}$, which has entries $h_{\_}{ij}=1/(i+j-1)$, since it appears to be a particularly challenging matrix to the rescaling method which we consider in Section 3.2.

For an ABC, the corresponding randomized computation with $Q$ recursions, $\hat{F}^{(Q)}$, was defined in (28). In this section, we will use $F^{(Q)}$ to denote the deterministic counterpart. It is defined in the same way, but with each $s_{\_}i^{(q)}(j)=1$ and each $\pi_{\_}i^{(q)}(j)=j$, i.e., with no randomness involved. As in Example 7, for EBCs we will use $\hat{G}^{(Q)}$ and $G^{(Q)}$ to denote the corresponding randomized and deterministic computations with $Q$ recursions.

3.1 Approximate algorithm

We create an ABC of the form (7) by taking the tensors $\bm{\mathscr{U}}$, $\bm{\mathscr{V}}$ and $\bm{\mathscr{W}}$ corresponding to Strassen’s algorithm and perturbing them: For each of the three tensors, we add i.i.d. mean zero Gaussian noise with standard deviation $10^{-3}$ to each element in the tensor equal to 1 as well as to five randomly selected elements that are equal to 0. Since we do the computations in double precision, and since we add a considerable amount of noise, these experiments are designed to test Propositions 3–6, i.e., we can ignore any floating point error.

In the first experiment, we draw two Gaussian matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{80\times 80}$ and compute

[TABLE]

for $n\in[10^{4}]$ and $Q\in[3]$. Here, $\hat{F}^{(Q)}_{\_}i$ is the $i$th realization of $\hat{F}^{(Q)}$. Figure 1 shows the results. As expected from Propositions 3 and 6, the quantity in (77) becomes smaller as $n$ increases.

In the second experiment, we draw two Gaussian matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{320\times 320}$ and compute

[TABLE]

for $i\in[100]$ and $Q\in[5]$, i.e., the relative error for each of 100 trials. The box plots in Figure 2 compares the empirical distribution of (78) to the relative error for the deterministic approximate algorithm, i.e., $\|F^{(Q)}(\mathbf{A},\mathbf{B})-\mathbf{A}\mathbf{B}\|/\|\mathbf{A}\mathbf{B}\|$. In this particular case, randomization does not impact the median error and there is very little variation between trials. Figure 3 repeats this experiment, but with $\mathbf{A}=\mathbf{B}=\mathbf{H}$, where $\mathbf{H}$ is the $320\times 320$ Hilbert matrix. In this case, the randomized scheme frequently results in a slightly larger error. In Figure 15–19 in Appendix A we provide additional results for when $\mathbf{A},\mathbf{B}$ are Gaussian, uniform and type 1–3 adversarial. Those results demonstrate that randomization can both increase and decrease the error in this setting. However, as expected from Propositions 4 and 5, the difference in error between the randomized and deterministic variants is typically not substantial.

To see how these results carry over to other ABCs, we also present results in Figure 4 from experiments which use an instance of the $12\times 12$ APA algorithm in [4] with $\varepsilon=\text{1e\textminus 4}$ in the representation (9). Results are shown for single recursion experiments on $12\times 12$ matrices drawn from different distributions. As for the ABC based on the perturbed variant of Strassen’s algorithm, the randomized variant of this ABC sometimes results in a smaller and sometimes in a larger error than the deterministic counterpart.

3.2 Exact algorithm in single precision floating point arithmetic

We first consider Strassen’s algorithm, without any perturbations so that it is exact, when the computations are done in single precision floating point arithmetic. In error computations, we use the double precision product for $\mathbf{A}\mathbf{B}$ computed using the standard algorithm as the true value of the product. Recall that by “standard algorithm” we mean the $O(n^{3})$ algorithm.

In the first experiment, we draw two Gaussian matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{80\times 80}$ and compute the quantity in (77), but with each $\hat{F}_{\_}i^{(Q)}$ replaced by $\hat{G}_{\_}i^{(Q)}$, where $\hat{G}_{\_}i^{(Q)}$ is the $i$th realization of $\hat{G}^{(Q)}$, for $n\in[10^{4}]$ and $Q\in[3]$. Figure 5 shows the results, where we also have included the error for the standard algorithm computed in single precision as a reference. Although it is clear that the randomized algorithms do not converge to the exact correct answer, it seems like their expectations perform better than the standard algorithm. Although the figure only shows the result for a specific random pair $(\mathbf{A},\mathbf{B})$, we get qualitatively similar results every time we draw a new Gaussian matrix pair. The fact that the average of a single, or a few, outcomes of the randomized algorithms performs worse than the standard algorithm is to be expected, since Strassen’s algorithm is more susceptible to numerical error that the standard algorithm. Figure 6 repeats this experiment, but with $\mathbf{A}=\mathbf{B}=\mathbf{H}$, where $\mathbf{H}$ is the $80\times 80$ Hilbert matrix. The results are similar to those in Figure 5. Figures 20–24 in the appendix provide additional results for when $\mathbf{A},\mathbf{B}$ are Gaussian, uniform and type 1–3 adversarial.

In the second experiment, we first compare the deterministic and three randomized versions of Strassen’s algorithm, as well as the rescaling scheme proposed in [1]. The first randomized version uses both random permutations and random signs. The two other randomized versions use only random signs and only random permutations, respectively. We include these two variants to better understand how random signs and random permutations each impact the performance. These three randomized methods will be referred to as “fully randomized,” “random sign” and “random permutation,” respectively. The algorithm that only uses random permutations corresponds to the method suggested in [11].

The purpose of the rescaling scheme in [1] is to improve numerical stability of fast EBCs, and has two steps: Outside scaling and inside scaling. With outside scaling, $\mathbf{C}=\mathbf{A}\mathbf{B}$ is computed via

[TABLE]

where $\mathbf{D}_{\_}\mathbf{A}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\textnormal{diag}(\max_{\_}{j}|a_{\_}{ij}|)$ and $\mathbf{D}_{\_}\mathbf{B}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\textnormal{diag}(\max_{\_}{i}|b_{\_}{ij}|)$. With inside scaling, $\mathbf{C}$ is instead computed via

[TABLE]

where $\mathbf{D}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\textnormal{diag}(\sqrt{\max_{\_}j|b_{\_}{kj}|/\max_{\_}i|a_{\_}{ik}|})$. In exact arithmetic, the scaling matrices in (79) and (80) will cancel out, in which case $\mathbf{C}_{\_}\text{outside}=\mathbf{C}_{\_}\text{inside}=\mathbf{C}$. Both outside and inside scaling can be applied at the same time, as well as multiple times in an alternating fashion. The rescaling scheme that works best in numerical experiments in [1] does outside-inside rescaling twice. We use the same rescaling scheme in our experiments, which we refer to as “Rescaled 2x O-I.” See Section 6 in [1], in particular Algorithm 3, for further details on the rescaling method.

For the experiment, we draw two random matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{320\times 320}$ and compute the quantity in (78), but with each $\hat{F}_{\_}i^{(Q)}$ replaced by $\hat{G}_{\_}i^{(Q)}$, for $i\in[100]$ and $Q\in[5]$. Figures 7–11 show the results for Gaussian, uniform and type 1–3 adversarial matrices. Figure 12 shows the result when $\mathbf{A}$ and $\mathbf{B}$ are both Hilbert matrices. These figures include the error of the standard algorithm computed in single precision as a reference.

When the matrices are Gaussian (Figure 7), all algorithms perform roughly the same with little variation between trials. For uniform matrices (Figure 8), the rescaling method performs about the same as the deterministic method. The fully randomized method has a lower error than the deterministic method, and it seems like this improvement comes from the random signs. For type 1 adversarial matrices (Figure 9), the rescaling method does remarkably well, achieving the same accuracy as the standard algorithm. The fully randomized algorithm outperforms the deterministic algorithm, and it seems like this improvement is coming from the random signs. For type 2 adversarial matrices (Figure 10), the fully randomized method will sometimes perform worse than the deterministic algorithm, but has a lower median error for 2 or more recursions. The rescaling method also improves on the deterministic method, although the median error for the fully randomized method is lower for 4 recursions or more. Type 3 adversarial matrices were specifically proposed in [1] to show a situation when rescaling does not work. This is clear in Figure 11, where the rescaling method has the same error as the deterministic method. Our fully randomized method, however, has a lower error, and it seems like both the random signs and random permutations contribute to this performance improvement. When the matrices are Hilbert matrices (Figure 12), the rescaling method does very poorly, with a much larger error than the deterministic method. Once again, our randomized method reduces the error compared to the deterministic method, with both the random signs and random permutations contributing to the improved performance. Figures 25–29 in the appendix provide additional results for when $\mathbf{A},\mathbf{B}$ are Gaussian, uniform, and type 1–3 adversarial.

In Figure 13, we give the results for some experiments which use an EBC for $12\times 12$ matrix multiplication derived as in [3] from the $12\times 12$ APA algorithm in [4] via the approach discussed in Section 1.1.222There appears to be a few typos in the definition of $w_{\_}r^{(s)}$ in Equation (5.2) in [3], which defines the APA scheme. We encourage the reader to consult our code for a corrected definition. The APA algorithm in question has an error tensor whose entries are polynomials with maximum degree $d=6$. Consequently, 7 distinct values of $\varepsilon$ are required in (10) to derive an exact scheme. For this purpose, we choose $\varepsilon_{\_}1,\varepsilon_{\_}2,\ldots,\varepsilon_{\_}7$ to be $0.1,0.2,\ldots,0.7$. Our results are for single recursion experiments on $12\times 12$ matrices drawn from different distributions. The results are similar to those for the Strassen EBC, except for the Gaussian case where the variability of the randomized algorithms is much higher. The variability of the randomized methods also appears to be somewhat higher on the other matrix types as well. Overall, our randomized methods perform favorably compared to the deterministic method, especially for the adversarial matrices and the Hilbert matrix. The rescaling method once again performs very well on type 1 adversarial matrices, but poorly on the Hilbert matrix.

Although our method has some variability in performance due to being randomized, it seems to reduce the error compared to the deterministic method more reliably on a wider range of matrices than the rescaling method does, especially when more recursions are used. Our method also works well with the Hilbert matrix, which leads to large errors for the rescaling method. One benefit of our randomized method is that it can be done exactly even in low precision arithmetic, since it only involves permutation of rows/columns and the flipping of signs. The rescaling method, on the other hand, involves diagonal matrices with floating point numbers along the diagonals, which adds another potential source for numerical error in the algorithm. These experiments also indicate that both random signs and random permutations may separately help to reduce the error, and that combining the two seems to lower the error further.

3.3 Randomization of ABC derived from APA algorithm

In Section 1.1 we discussed APA algorithms and how to derive EBCs from them by taking a linear combination of a few instances of the APA algorithm following an idea of Bini [3]. Although the EBCs derived in this fashion are exact mathematically, they may suffer from greater numerical error than the standard $O(n^{3})$ algorithm due to cancellation. In this subsection, we compare an instance of the EBC in [3], which is derived from the APA algorithm in [4], to a randomized version of the ABC we get by fixing the error parameter in the same APA algorithm. The goal with these experiments is to see if the expectation of our randomized ABC can perform better than the EBC. For the EBC, which is given by (10), we choose $\varepsilon_{\_}1,\varepsilon_{\_}2,\ldots,\varepsilon_{\_}7$ to be $0.1,0.2,\ldots,0.7$. To get an ABC, we choose $\varepsilon=\text{1.5e\textminus 2}$ in the APA algorithm. This ABC is then randomized as in Definition 2. The experiments are similar to that in Figure 1, but with an added baseline which shows the performance of the EBC. The experiments are done for a single recursion on $12\times 12$ matrices. All computations are done in single precision arithmetic. Figure 14 shows the results, which each subplot showing the outcome for a different kind of matrix. Based on these result, the expectation of the randomized ABC seems to perform better than the EBC.

4 Conclusion

In this paper we have suggested an approach for randomizing formulas for bilinear computation of matrix products which does not increase the asymptotic computational complexity. We have considered the implications of this approach when there are two sources of error: The first due to the algorithm itself being only approximately correct, and the second due to numerical error from using floating point arithmetic. We believe that our results are encouraging, and provide ideas for improving the properties of matrix multiplication when these error sources are present separately or together.

An interesting area for future research is to investigate other methods of randomization (e.g. combining the random sign changes in this paper with the fast Hadamard transform) and see if such a method can further improve the results. It would also be interesting to investigate how fast randomized approximate methods for matrix multiplication can be used in applications. One potential application areas is for computations in neural networks, where some amount of error in the computation is acceptable, and where randomization may help improve robustness of the trained model. It would also be interesting to investigate how our randomization scheme can be used as a component in fast linear algebra algorithms, such as recursive dense matrix inversion; see [13] for details.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant No. ECCS-1810314.

Appendix A Additional experiments

First, we repeat the second experiment we did in Section 3.1 for the Strassen ABC with multiple Gaussian, uniform, and type 1–3 adversarial matrices. We use the same setup as in the main manuscript: For each experiment, we create an approximate algorithm by perturbing Strassen’s algorithm, we draw random matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{320\times 320}$, and then we compute the quantity in (78) for $i\in[100]$ and $Q\in[5]$ and compare it to the relative error for the deterministic approximate algorithm. Figures 15–19 show the results. In the case of Gaussian matrices (Figure 15), the results look very similar to those in the main manuscript, with almost no difference in error between the deterministic and the randomized approximate algorithms. Figures 16–19 (uniform and type 1–3 adversarial) show that randomization can both increase and decrease the error. Note that both the deterministic and randomized approximate algorithms do particularly poorly on type 1 adversarial matrices.

Next, we repeat the first experiment in Section 3.2 with multiple Gaussian, uniform and type 1–3 adversarial matrices. We use the same setup as in the main manuscript: For each experiment, we draw random matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{80\times 80}$ and compute the quantity in (77), but with each $\hat{F}_{\_}i^{(Q)}$ replaced by $\hat{G}_{\_}i^{(Q)}$, for $n\in[10^{4}]$ and $Q\in[3]$. Figures 20–24 show the results. Although using more recursions seems to increase the error for a single realization of the algorithm, it also seems to reduce the error of the expectation of the computed matrix product.

Finally, we repeat the second experiment in Section 3.2 for the Strassen EBC with multiple Gaussian, uniform and type 1–3 adversarial matrices. We use the same setup as in the main manuscript: For each experiment, we draw random matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{320\times 320}$ and compute the quantity in (78), but with each $\hat{F}_{\_}i^{(Q)}$ replaced by $\hat{G}_{\_}i^{(Q)}$, for $i\in[100]$ and $Q\in[5]$. We do the same for variants of the algorithm that only use random permutations or random sign functions. We compare these to the deterministic version $G^{(Q)}$, and also include the error of the standard algorithm computed in single precision as a reference. Figures 25–29 show the results, which look very similar to those presented in the main manuscript.