# skew-sparse matrix multiplication

Qiao-Long Huang1, Ke Ye2, Xiao-Shan Gao2
1 Research Center for Mathematics and Interdisciplinary Sciences, Shandong
University
2 KLMM, UCAS, Academy of Mathematics and Systems Science, Chinese Academy of
Sciences

###### Abstract.

Based on the observation that $\mathbb{Q}^{(p-1)\times(p-1)}$ is isomorphic to
a quotient skew polynomial ring, we propose a new method for
$(p-1)\times(p-1)$ matrix multiplication over $\mathbb{Q}$, where $p$ is a
prime number. The main feature of our method is the acceleration for matrix
multiplication if the product is skew-sparse. Based on the new method, we
design a deterministic algorithm with complexity
$\operatorname{O}(T^{\omega-2}p^{2})$, where $T\leq p-1$ is a parameter
determined by the skew-sparsity of input matrices and $\omega$ is the
asymptotic exponent of matrix multiplication. Moreover, by introducing
randomness, we also propose a probabilistic algorithm with complexity
$\operatorname{O}^{\thicksim}(t^{\omega-2}p^{2}+p^{2}\log\frac{1}{\nu})$,
where $t\leq p-1$ is the skew-sparsity of the product and $\nu$ is the
probability parameter.

##  1\. introduction

Fast algorithms for matrix multiplication are essential ingredients in
numerous fundamental computational problems such as linear system solving
[33], matrix inversion [5], determinant evaluation [24] and Boolean matrix
multiplication [12, 21]. Thus the computational complexity of matrix
multiplication is one of the most important problems in theoretical computer
science and mathematics, which attracts extensive attention from researchers
and practitioners in various areas.

The complexity of matrix multiplication is usually measured by _the asymptotic
exponent_ $\omega$, which is defined by

$\omega\coloneqq\liminf_{n\to\infty}\log_{n}M_{n},$

where $M_{n}$ is the total cost of operations for $n\times n$ matrix
multiplication. It is clear from the definition that $2\leq\omega\leq 3$. For
a relatively long time, it was believed that $\omega=3$. The seminal work [29]
of Strassen for the first time proved that $\omega\leq 2.81$, after which the
upper bound of $\omega$ was further improved by a series of works [20, 27, 3,
30]. Using a powerful technique called the “laser method” developed in [9, 34,
14], Alman and Williams were able to obtain the state-of-the-art upper bound
$\omega\leq 2.37286$ [1]. Unfortunately, it was proved that the upper bound of
$\omega$ one can obtain by the “laser method” is at least $2.3078$ [2],
indicating the necessity of developing a new method.

Besides the complexity of matrix multiplication for general matrices, the same
problem for structured matrices is also of great importance in both
theoretical studies [19, 28, 36, 8, 35, 25] and practical applications[6, 18,
26, 31, 17, 11]. The most commonly considered structured matrices are
circulant, Toeplitz, Hankel, Cauchy, symmetric, sparse matrices and their
variants.

In this paper, we discuss fast matrix multiplication for a new type of
structured matrices called _skew-sparse_ matrices, which are defined by sparse
skew polynomials. It is noticeable that the complexity of fast sparse matrix
multiplication algorithms [19, 28, 36] can be even higher than
$\operatorname{O}(n^{2.37286})$ when applied to dense matrices. Our algorithm
for skew-sparse matrix multiplication, as a comparison, has complexity at most
$\operatorname{O}(n^{\omega})$.

##  2\. preliminaries

In this section, we provide some preliminaries on cyclotomic field, normal
basis and skew polynomials.

###  2.1. cyclotomic field

Let $p=n+1$ be a prime and let $\beta$ be a $p$-th primitive root of unity.
Hence we have

$\beta^{p-1}+\cdots+\beta+1=0.$

The extension field $\mathbb{Q}(\beta)$ of $\mathbb{Q}$ is called the
_cyclotomic field_ generated by $\beta$. It is clear that
$[\mathbb{Q}(\beta):\mathbb{Q}]=p-1$.

A _primitive root of $\mathbb{Z}_{p}$_ is an integer $r\in\\{0,1,\dots,p-1\\}$
such that $\langle r\rangle=\mathbb{Z}^{\times}_{p}$. Given a primitive root
of $\mathbb{Z}_{p}$, we consider a $\mathbb{Q}$-homomorphism

$\sigma:\mathbb{Q}(\beta)\to\mathbb{Q}(\beta),\quad\beta\mapsto\beta^{r}.$ (1)

Since $r$ is a primitive root of $\mathbb{Z}_{p}$, we have

$\\{1,2,\dots,p-1\\}\equiv\\{1,r,\dots,r^{p-2}\\}\pmod{p}.$

This implies that $\sigma$ permutes $\beta,\beta^{2},\dots,\beta^{p-1}$, which
are roots of the irreducible polynomial $x^{p-1}+x^{p-2}+\cdots+x+1$. Thus
$\sigma$ is an automorphism of $\mathbb{Q}(\beta)$ over $\mathbb{Q}$, i.e.,
$\sigma\in\operatorname{Aut}(\mathbb{Q}(\beta)/\mathbb{Q})$.

###  2.2. normal basis

We recall that given a Galois extension $F\subseteq K$ with the Galois group
$G$. A _normal basis_ is a basis of $K/F$, which is of the form
$\\{g(\beta):g\in G\\}$ for some $\beta\in K$. We notice that the Galois group
of $\mathbb{Q}(\beta)/\mathbb{Q}$ is
$\operatorname{Gal}(\mathbb{Q}(\beta)/\mathbb{Q})=\langle\sigma\rangle$, where
$\sigma$ is the automorphism defined in (1). For each $i=1,\dots,p-1$, we
denote

$v_{i}\coloneqq\sigma^{i}(\beta)=\beta^{r^{i-1}}.$ (2)

It is straightforward to verify that $v_{i+1}=\sigma(v_{i})$. According to the
discussion in Subsection 2.1 and the construction of $v_{i}$’s, we have the
following:

######  Lemma 2.1.

[15, page 4] $\\{v_{1},v_{2},\dots,v_{p-1}\\}$ is a normal basis of
$\mathbb{Q}(\beta)/\mathbb{Q}$.

Since each $v_{i}$ is a power of $\beta$, with respect to the basis
$\\{v_{1},v_{2},\dots,v_{p-1}\\}$, there exist fast algorithms [4] to compute
multiplication, division, addition and subtraction in $\mathbb{Q}(\beta)$.
Moreover, by identifying $\mathbb{Q}(\beta)$ with $\mathbb{Q}[z]/(H(z))$,
where $H(z)=z^{p-1}+\cdots+1$, we have the following observation regarding the
complexity of arithmetic operations in $\mathbb{Q}(\beta)$:

######  Lemma 2.2.

[22] One arithmetic operation in $\mathbb{Q}(\beta)$ equals to
$\operatorname{O}^{\thicksim}(p)$ arithmetic operations in $\mathbb{Q}$.

###  2.3. skew polynomials

Let $x$ be an indeterminate and let $\sigma$ be the automorphism defined in
(1). We denote by $\mathbb{Q}(\beta)[x;\sigma]$ the ring
$(\mathbb{Q}(\beta)[x],+,\ast_{\sigma})$, where $+$ is the usual polynomial
addition and $\ast_{\sigma}$ is induced by
$x\ast_{\sigma}c\coloneqq\sigma(c)x$ and $c\ast_{\sigma}x\coloneqq cx$. More
precisely, we have

$\displaystyle\sum_{i=0}^{d}a_{i}x^{i}+\sum_{i=0}^{d}b_{i}x^{i}$
$\displaystyle=\sum_{i=0}^{d}(a_{i}+b_{i})x^{i},$
$\displaystyle\left(\sum_{j=0}^{d}a_{j}x^{j}\right)\ast_{\sigma}\left(\sum_{k=0}^{e}b_{k}x^{k}\right)$
$\displaystyle=\sum_{\ell=0}^{d+e}\sum_{s=0}^{\ell}a_{s}\sigma^{s}(b_{\ell-s})x^{\ell}.$

Here in the second formula, we adopt the convention that $a_{i}=0$ (resp.
$b_{i}=0$) if $i>d$ (resp. $i>e$). $\mathbb{Q}(\beta)[x;\sigma]$ is called
_the ring of skew polynomials_. According to [7, Lemma 1.4], there exists an
isomorphism of algebras:

$\displaystyle\varepsilon:\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$
$\displaystyle\to\mathrm{End}_{\mathbb{Q}}(\mathbb{Q}(\beta)),$ (3)
$\displaystyle\sum_{i=0}^{p-2}a_{i}x^{i}$
$\displaystyle\mapsto\sum_{i=0}^{d}a_{i}\sigma^{i}.$

Here an element in $\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$ is written as a
polynomial of degree at most $p-2$ and we adopt this convention in the rest of
the paper.

####  2.3.1. sumset of skew polynomials

In the sequel, we need the notion of the sparsity of a skew polynomial. Given
a polynomial
$f=\sum_{i=1}^{t}a_{i}x^{e_{i}}\in\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$
where $0\leq e_{1}<\cdots<e_{t}\leq p-2$ and all $a_{i}\neq 0,i=1,\dots,t$,
the set $\mathrm{supp}(f)\coloneqq\\{e_{1},\dots,e_{t}\\}$ is called the
_support_ of $f$. Consequently, we define the _sparsity_ of $f$ to be
$\\#f\coloneqq\\#\mathrm{supp}(f)=t$.

Given two polynomials $f$ and $g$ in
$\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$, we define the _sumset of $f$ and
$g$_ by

$\mathbb{S}(f,g)\coloneqq\\{e_{f}+e_{g}:e_{f}\in\mathrm{supp}(f),e_{g}\in\mathrm{supp}(g)\\}\pmod{p-1}.$
(4)

######  Lemma 2.3.

We have the following:

  1. (i)

$\mathbb{S}(f,g)$ contains $\mathrm{supp}(f\ast_{\sigma}g)$;

  2. (ii)

$\\#(f\ast_{\sigma}g)\leq\\#\mathbb{S}(f,g)\leq\\#f\cdot\\#g$ and
$\\#\mathbb{S}(f,g)\leq p-1$;

  3. (iii)

the strict inequality $\\#(f\ast_{\sigma}g)<\\#\mathbb{S}(f,g)$ holds only in
the occurrence of coefficient cancellations.

####  2.3.2. evaluation of a skew polynomial

According to (3), a skew polynomial
$f\in\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$ defines a $\mathbb{Q}$-linear
mapping $\epsilon(f):\mathbb{Q}(\beta)\rightarrow\mathbb{Q}(\beta)$ obtained
by evaluating $f$ at $\sigma$. For simplicity, we denote $\epsilon(f)(b)$ by
$f(b)$ for each $b\in\mathbb{Q}(\beta)$.

Next we briefly discuss a fast approach to simultaneously evaluate $f(b)$ for
$b\in\mathbb{Q}(\beta)$. Since $\mathbb{Q}(\beta)$ is a $p-1$ dimensional
vector space over $\mathbb{Q}$, we may write

$b=\sum_{j=1}^{p-1}b_{j}v_{j},$

where $b_{j}\in\mathbb{Q}$ and $\\{v_{1},\dots,v_{p-1}\\}$ is the normal basis
defined in (2). In particular, we have $f(v_{i})=\sum_{j=1}^{p-1}a_{ij}v_{j}$
for some $a_{ij}\in\mathbb{Q}$, $1\leq i,j\leq p-1$. Thus with respect to the
normal basis $\\{v_{j}\\}_{j=1}^{p-1}$, the $\mathbb{Q}$-linear map
$\epsilon(f)$ can be represented by the $(p-1)\times(p-1)$ matrix over
$\mathbb{Q}$:

$\epsilon(f)=\begin{bmatrix}a_{11}&a_{12}&\cdots&a_{1,p-1}\\\
a_{21}&a_{22}&\cdots&a_{2,p-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{p-1,1}&a_{p-1,2}&\cdots&a_{p-1,p-1}\\\ \end{bmatrix}.$

The above observation enables us to write (3) in a more explicit way. Namely,
we have

$\displaystyle\varphi:\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$
$\displaystyle\to\mathbb{Q}^{(p-1)\times(p-1)}$ (5) $\displaystyle f$
$\displaystyle\mapsto(a_{ij})_{i,j=1}^{p-1}$

where $f(v_{i})=\sum_{j=1}^{p-1}a_{ij}v_{j},1\leq i\leq p-1$.

Accordingly, the evaluation of $f$ at
$b=\sum_{j=1}^{p-1}b_{j}v_{j}\in\mathbb{Q}(\beta)$ can be written as

$f(b)=\begin{bmatrix}b_{1}&\cdots&b_{p-1}\end{bmatrix}\begin{bmatrix}f(v_{1})\\\
\vdots\\\ f(v_{p-1})\\\
\end{bmatrix}=\begin{bmatrix}b_{1}&\cdots&b_{p-1}\end{bmatrix}\begin{bmatrix}a_{1,1}&\cdots&a_{1,p-1}\\\
\vdots&\ddots&\vdots\\\ a_{p-1,1}&\cdots&a_{p-1,p-1}\\\
\end{bmatrix}\begin{bmatrix}v_{1}\\\ \vdots\\\ v_{p-1}\\\ \end{bmatrix}.$

It is straightforward to verify that for
$b_{1},\dots,b_{t}\in\mathbb{Q}(\beta)$, we have

$\left[\begin{array}[]{c}f(b_{1})\\\ f(b_{2})\\\ \vdots\\\ f(b_{t})\\\
\end{array}\right]=\begin{bmatrix}b_{1,1}&\cdots&b_{1,p-1}\\\
\vdots&\ddots&\vdots\\\ b_{t,1}&\cdots&b_{t,p-1}\\\
\end{bmatrix}\begin{bmatrix}a_{1,1}&\cdots&a_{1,p-1}\\\
\vdots&\ddots&\vdots\\\ a_{p-1,1}&\cdots&a_{p-1,p-1}\\\
\end{bmatrix}\begin{bmatrix}v_{1}\\\ \vdots\\\ v_{p-1}\\\ \end{bmatrix},$ (6)

where $b_{i}=\sum_{j=1}^{p-1}b_{ij}v_{j}$ and $1\leq i\leq t$. We notice that
in (6), the two matrices have size of $t\times(p-1)$ and $(p-1)\times(p-1)$
respectively. Thus we have the following:

######  Lemma 2.4.

[16] The evaluation of $f(b_{1}),\dots,f(b_{t})$ can be done by
$\operatorname{O}(t^{\omega-2}p^{2})$ operations in $\mathbb{Q}$, where
$\omega$ is the exponent of the matrix multiplication.

####  2.3.3. interpolation of a skew polynomial

Assume $f\in\mathbb{Q}(\beta)[x;\sigma]$ and
$\mathrm{supp}(f)=\\{e_{1},\dots,e_{t}\\}$. We discuss how to find
coefficients of $f$ from $f(1),f(v_{1}),\dots,f(v_{1}^{t-1})$.

By assumption, $f$ has the form $f=\sum_{i=1}^{t}c_{i}x^{e_{i}}$. Since

$\sigma^{k}(v^{i}_{1})=(\sigma^{k}(v_{1}))^{i}=(v_{k+1})^{i},\quad
k\in\mathbb{N},$

we have

$f(v^{i}_{1})=c_{1}(v_{e_{1}+1})^{i}+\cdots+c_{t}(v_{e_{t}+1})^{i},\quad
i=0,1,\dots,t-1.$

Here indexes $e_{1}+1,\dots,e_{t}+1$ are are taken modulo $p$. This implies

$\left[\begin{array}[]{c}f(v^{0}_{1})\\\ f(v^{1}_{1})\\\ \vdots\\\
f(v^{t-1}_{1})\\\
\end{array}\right]=\left[\begin{array}[]{cccc}1&1&\cdots&1\\\
v_{e_{1}+1}&v_{e_{2}+1}&\cdots&v_{e_{t}+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
v^{t-1}_{e_{1}+1}&v^{t-1}_{e_{2}+1}&\cdots&v^{t-1}_{e_{t}+1}\\\
\end{array}\right]\left[\begin{array}[]{c}c_{1}\\\ c_{2}\\\ \vdots\\\ c_{t}\\\
\end{array}\right].$ (7)

To recover $c_{1},\dots,c_{t}$, we need to solve the linear system (7), whose
coefficient matrix is Vandermonde. Therefore we obtain the result that
follows.

######  Lemma 2.5.

[16] The interpolation of $f$ can be done by $\operatorname{O}^{\sim}(tn)$
operations in $\mathbb{Q}$.

##  3\. deterministic matrix multiplication via skew polynomials

###  3.1. relation between matrices and skew polynomials

Let
$\varphi:\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)\to\mathbb{Q}^{(p-1)\times(p-1)}$
be the $\mathbb{Q}$-algebra isomorphism defined in (5). In this subsection, we
discuss algorithms computing $\varphi$ and $\varphi^{-1}$. According to the
definition of $\varphi$ and properties of normal basis, it is easy to
establish the following relation between skew polynomials and matrices.

######  Lemma 3.1.

Let
$f=\sum_{j=0}^{p-2}\mu_{j+1}x^{j}\in\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$
be a skew polynomial and let $C=(c_{ij})_{i,j=1}^{p-1}$ be a
$(p-1)\times(p-1)$ matrix over $\mathbb{Q}$. Then $\varphi(f)=C$ if and only
if

$\begin{bmatrix}v_{1}&v_{2}&\cdots&v_{p-1}\\\ v_{2}&v_{3}&\cdots&v_{1}\\\
\vdots&\vdots&\ddots&\vdots\\\ v_{p-1}&v_{1}&\cdots&v_{p-2}\\\
\end{bmatrix}\begin{bmatrix}\mu_{1}\\\ \mu_{2}\\\ \vdots\\\ \mu_{p-1}\\\
\end{bmatrix}=\begin{bmatrix}c_{11}&c_{12}&\cdots&c_{1,p-1}\\\
c_{21}&c_{22}&\cdots&c_{2,p-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
c_{p-1,1}&c_{p-1,2}&\cdots&c_{p-1,p-1}\\\ \end{bmatrix}\begin{bmatrix}v_{1}\\\
v_{2}\\\ \vdots\\\ v_{p-1}\\\ \end{bmatrix},$ (8)

where $\\{v_{1},\dots,v_{p-1}\\}$ is the normal basis defined in (2).

####  3.1.1. from matrices to skew polynomials

Given a $(p-1)\times(p-1)$ matrix $C=(c_{ij})_{i,j=1}^{p-1}$ over
$\mathbb{Q}$, we may compute the skew polynomial $f\coloneqq\varphi^{-1}(C)$.

######  Lemma 3.2.

Let $\\{v_{i}\\}_{i=1}^{p-1}$ be as above. We denote

$V\coloneqq\begin{bmatrix}v_{1}&v_{2}&\cdots&v_{p-1}\\\
v_{2}&v_{3}&\cdots&v_{1}\\\ \vdots&\vdots&\ddots&\vdots\\\
v_{p-1}&v_{1}&\cdots&v_{p-2}\end{bmatrix},\quad
W\coloneqq\begin{bmatrix}\frac{1}{v_{1}}-1&\frac{1}{v_{2}}-1&\cdots&\frac{1}{v_{p-1}}-1\\\
\frac{1}{v_{2}}-1&\frac{1}{v_{3}}-1&\cdots&\frac{1}{v_{1}}-1\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{1}{v_{p-1}}-1&\frac{1}{v_{1}}-1&\cdots&\frac{1}{v_{p-2}}-1.\end{bmatrix}$

Then we have $VW=pI_{p-1}$.

###### Proof.

We proceed by directly computing the $(i,j)$-th entry $a_{ij}$ of $VW$. We
want to prove that

$a_{ij}=p\delta_{ij},$

where $\delta_{ij}$ is the Kronecker delta. We observe that the $i$-th row of
$V$ is $(v_{i},v_{i+1},\dots,v_{i+n-1})$ and $j$-th column of $W$ is
$(\frac{1}{v_{j}}-1,\frac{1}{v_{j+1}}-1,\dots,\frac{1}{v_{j+p-2}}-1)^{T}$.
Thus we have

$a_{ij}=\sum_{k=0}^{p-2}v_{i+k}\left(\frac{1}{v_{j+k}}-1\right)=\sum_{k=0}^{p-2}\frac{v_{i+k}}{v_{j+k}}-\sum_{k=0}^{p-2}v_{i+k}=\sum_{k=0}^{p-2}\sigma^{k}\left(\frac{v_{i}}{v_{j}}\right)-\sum_{k=0}^{p-2}v_{k}.$

Here indexes are taken modulo $p$. Since $\sum_{j=0}^{p-1}\beta^{j}=0$, we
conclude that for each $i\in\mathbb{Z}$,

$\sum_{k=0}^{p-2}v_{i+k}=\sum_{k=0}^{p-2}v_{k}=-1.$

If $i=j$, then $a_{ij}=\left(\sum_{k=0}^{p-2}\sigma^{k}(1)\right)+1=p$. If
$i\neq j$, then $\frac{v_{i}}{v_{j}}=\beta^{r^{i-1}-r^{j-1}}\neq 0$. As $r$ is
a primitive root of $\mathbb{Z}_{p}$, there exists an
$m\in\\{1,2,\dots,p-1\\}$ such that $r^{m-1}=r^{i-1}-r^{j-1}\pmod{p}$. Hence
$\frac{v_{i}}{v_{j}}=v_{m}$ and

$a_{ij}=\left(\sum_{k=0}^{p-2}\sigma^{k}(v_{m})\right)+1=\left(\sum_{k=0}^{p-2}v_{m+k}\right)+1=0.$

∎

By a combination of Lemmas 3.1 and 3.2, we derive the following formula for
$\varphi^{-1}(C)$.

######  Lemma 3.3.

Let $C=(c_{ij})_{i,j=1}^{p-1}$ be a $(p-1)\times(p-1)$ matrix over
$\mathbb{Q}$. Then $\varphi^{-1}(C)=\sum_{j=0}^{p-2}\mu_{j+1}x^{j}$ where

$\begin{bmatrix}\mu_{1}\\\ \mu_{2}\\\ \vdots\\\
\mu_{p-1}\end{bmatrix}=\frac{1}{p}\begin{bmatrix}\frac{1}{v_{1}}-1&\frac{1}{v_{2}}-1&\cdots&\frac{1}{v_{p-1}}-1\\\
\frac{1}{v_{2}}-1&\frac{1}{v_{3}}-1&\cdots&\frac{1}{v_{1}}-1\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{1}{v_{p-1}}-1&\frac{1}{v_{1}}-1&\cdots&\frac{1}{v_{p-2}}-1\end{bmatrix}\begin{bmatrix}c_{11}&c_{12}&\cdots&c_{1,p-1}\\\
c_{21}&c_{22}&\cdots&c_{2,p-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
c_{p-1,1}&c_{p-1,2}&\cdots&c_{p-1,p-1}\end{bmatrix}\begin{bmatrix}v_{1}\\\
v_{2}\\\ \vdots\\\ v_{p-1}\end{bmatrix}.$ (9)

######  Proposition 3.4.

The skew polynomial $\varphi^{-1}(C)$ can be computed by
$\operatorname{O}^{\sim}(p^{2})$ operations in $\mathbb{Q}$.

###### Proof.

Since (9) only involves the $(p-1)\times(p-1)$ matrix vector product, the
complexity is at most $\operatorname{O}^{\sim}(p^{2})$. ∎

####  3.1.2. from skew polynomials to matrices

Assume $f\in\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$ has the form

$f:=\mu_{1}+\mu_{2}x+\mu_{3}x^{2}+\cdots+\mu_{p-1}x^{p-2}.$ (10)

Now we present the procedure $\varphi$ that transforms the skew polynomial
into the matrix in pseudocode.

Algorithm 1 matrices to skew polynomials

1:skew polynomial $f=\sum_{j=0}^{p-2}\mu_{j+1}x^{j}$ and normal basis
$\\{v_{1},\dots,v_{p-1}\\}$ of $\mathbb{Q}(\beta)/\mathbb{Q}$.

2:the matrix $\varphi(f)=(c_{ij})_{i,j=1}^{p-1}$.

3: compute  $\begin{bmatrix}\beta_{1}\\\ \beta_{2}\\\ \vdots\\\ \beta_{p-1}\\\
\end{bmatrix}=\begin{bmatrix}v_{1}&v_{2}&\cdots&v_{p-1}\\\
v_{2}&v_{3}&\cdots&v_{1}\\\ \vdots&\vdots&\ddots&\vdots\\\
v_{p-1}&v_{1}&\cdots&v_{p-2}\\\ \end{bmatrix}\begin{bmatrix}\mu_{1}\\\
\mu_{2}\\\ \vdots\\\ \mu_{p-1}\\\ \end{bmatrix}.$

4: compute $c_{ij}\in\mathbb{Q},1\leq i,j\leq p-1$ such that
$\begin{bmatrix}\beta_{1}\\\ \beta_{2}\\\ \vdots\\\ \beta_{p-1}\\\
\end{bmatrix}=\begin{bmatrix}c_{11}&c_{12}&\cdots&c_{1,p-1}\\\
c_{21}&c_{22}&\cdots&c_{2,p-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
c_{p-1,1}&c_{p-1,2}&\cdots&c_{p-1,p-1}\\\ \end{bmatrix}\begin{bmatrix}v_{1}\\\
v_{2}\\\ \vdots\\\ v_{p-1}\\\ \end{bmatrix},$

######  Proposition 3.5.

Algorithm 1 returns $\varphi(f)$ using $\operatorname{O}^{\sim}(p^{2})$
operations in $\mathbb{Q}$.

###### Proof.

Since the matrix in Step 3 is a circulant matrix, we compute the product in
time $\operatorname{O}^{\sim}(p)$ over $\mathbb{Q}(\beta)$, which is
$\operatorname{O}^{\sim}(p^{2})$ operations in $\mathbb{Q}$. ∎

###  3.2. matrix multiplication via skew polynomials

Now we are ready to present our algorithm for matrix multiplication. The main
idea is to reduce the matrix multiplication to the skew polynomial
multiplication.

Algorithm 2 matrix multiplication via skew polynomials

1:matrices $A,B\in\mathbb{Q}^{(p-1)\times(p-1)}$, where $p$ is a prime.

2:the product $AB$.

3:find a $p$-th primitive root of unity $\beta$.

4:find a primitive root $r$ of $\mathbb{Z}_{p}$ and compute
$v_{i}=\beta^{r^{i-1}},i=1,\dots,p-1$.

5:compute $\varphi^{-1}(A),\varphi^{-1}(B)$ by (9) w.r.t the normal basis
$\\{v_{1},\dots,v_{p-1}\\}$.

6:compute the sumset
$\mathbb{S}(\varphi^{-1}(A),\varphi^{-1}(B))=\\{e_{1},e_{2},\dots,e_{t}\\}$.

7:compute
$\mathrm{Eval}_{i}:=\left(\varphi^{-1}(A)\ast_{\sigma}\varphi^{-1}(B)\right)(v^{i}_{1}),i=0,1,\dots,t-1$.

8:solve  $\begin{bmatrix}\mathrm{Eval}_{0}\\\ \mathrm{Eval}_{1}\\\ \vdots\\\
\mathrm{Eval}_{t-1}\\\ \end{bmatrix}=\begin{bmatrix}1&1&\cdots&1\\\
v_{e_{1}+1}&v_{e_{2}+1}&\cdots&v_{e_{t}+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
v^{t-1}_{e_{1}+1}&v^{t-1}_{e_{2}+1}&\cdots&v^{t-1}_{e_{t}+1}\end{bmatrix}\begin{bmatrix}c_{1}\\\
c_{2}\\\ \vdots\\\ c_{t}\end{bmatrix}$ and set $f\coloneqq
c_{1}x^{e_{1}}+\cdots+c_{t}x^{e_{t}}$.

9:compute $\varphi(f)$ by Algorithm 1.

######  Proposition 3.6.

Algorithm 2 computes $AB$ using $\operatorname{O}^{\sim}(t^{\omega-2}p^{2})$
operations in $\mathbb{Q}$, where $t$ is the cardinality of the sumset of
$\varphi^{-1}(A)$ and $\varphi^{-1}(B)$.

###### Proof.

If we denote by $m$ (resp. $\ast_{\sigma}$) the matrix (resp. skew-polynomial)
multiplication, then we have the following commutative diagram, which
validates Algorithm 2.

${\mathbb{Q}^{(p-1)\times(p-1)}\times\mathbb{Q}^{(p-1)\times(p-1)}}$${\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)\times\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)}$${\mathbb{Q}^{(p-1)\times(p-1)}}$${\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)}$$\scriptstyle{m}$$\scriptstyle{\varphi^{-1}\times\varphi^{-1}}$$\scriptstyle{\ast_{\sigma}}$$\scriptstyle{\varphi}$

Next we analyse the complexity of Algorithm 2. If $r$ is a primitive root of
$p$, then $r^{1},r^{2},\dots,r^{p-1}\pmod{p}$ must generate distinct integers
from $1$ to $(p-1)$. Checking the membership of a number in
$\\{2,\dots,p-2\\}$ has complexity at most $\operatorname{O}(p^{2})$ over
$\mathbb{Q}$. By Proposition 3.4, the complexity of Step 5 in Algorithm 2 is
$\operatorname{O}^{\sim}(p^{2})$ over $\mathbb{Q}$. It is straightforward to
verify that the complexity of Step 6 is $\operatorname{O}(p^{2})$ over
$\mathbb{Q}$. The cost of Step 7 is
$\operatorname{O}^{\sim}(t^{\omega-2}p^{2})$ $\mathbb{Q}$-operations by Lemma
2.4. In Step 8, since the coefficient matrix is a Vandermonde matrix, solving
this linear system costs $\operatorname{O}^{\sim}(tp)$ $\mathbb{Q}$-operations
[23]. By Proposition 3.5, Step 9 costs $\operatorname{O}^{\thicksim}(p^{2})$
operations in $\mathbb{Q}$. ∎

###  3.3. skew-sparse matrix multiplication

Suppose $p$ is a prime number and $r$ is a primitive root of $\mathbb{Z}_{p}$.
Let $1\leq k\leq p-1$ be the integer such that $r^{k-1}\equiv p-1\pmod{p}$.
For each $1\leq j\neq k\leq p-1$, we denote by $s_{j}$ the integer between $1$
and $(p-1)$ such that $r^{s_{j}-1}\equiv r^{j-1}+1\pmod{p}$. By Lemma 2.1,
$\\{v_{j}\\}_{j=1}^{p-1}$ is a normal basis of $\mathbb{Q}(\beta)/\mathbb{Q}$
where $v_{i}=\beta^{r^{j-1}},1\leq j\leq p-1$.

We construct two $(p-1)\times(p-1)$ matrices

$X\coloneqq\begin{bmatrix}0&1&0&\cdots&0\\\ 0&0&1&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&\cdots&1\\\
1&0&0&\cdots&0\end{bmatrix},\quad Y\coloneqq\begin{bmatrix}E_{s_{1}}\\\
\vdots\\\ E_{s_{k-1}}\\\ J\\\ E_{s_{k+1}}\\\ \vdots\\\ E_{s_{p-1}}\\\
\end{bmatrix},$

where $E_{i}$ is the row vector whose elements are all $0$ except the $i$-th,
which is equal to $1$, and $J$ is the row vector whose elements are all equal
to $-1$.

######  Lemma 3.7.

The set

$\\{X^{i}Y^{j}:0\leq i\leq p-2,1\leq j\leq p-1\\}$

is a $\mathbb{Q}$-basis of $\mathbb{Q}^{(p-1)\times(p-1)}$. In particular, $X$
and $Y$ are generators of $\mathbb{Q}^{(p-1)\times(p-1)}$ as a
$\mathbb{Q}$-algebra.

###### Proof.

Let $\beta$ be a $p$-th primitive root of unity. According to (5), there
exists an isomorphism of algebras:

$\varphi:\mathbb{Q}(\beta)[x,\sigma]/(x^{p-1}-1)\xrightarrow{\sim}\mathbb{Q}^{(p-1)\times(p-1)}.$

It is straightforward to verify that $\varphi(x)=X$ and $\varphi(\beta)=Y$ and
our claim follows easily. ∎

######  Corollary 3.8.

$\sum_{j=1}^{p-1}Y^{j}=-I$.

###### Proof.

Since $\beta^{p-1}+\cdots+\beta+1=0$, we have
$Y^{p-1}+\cdots+Y+I=\varphi(\beta^{p-1}+\cdots+\beta+1)=0$. ∎

######  Proposition 3.9.

The circulant matrices are in a one-to-one correspondence with elements in
$\mathbb{Q}(\beta)[x,\sigma]/(x^{p-1}-1)$ with coefficients in $\mathbb{Q}$.

###### Proof.

Assume $f=\sum_{j=0}^{p-2}\mu_{j+1}x^{j}$ is a skew polynomial with
coefficients in $\mathbb{Q}$. Then

$\varphi(f)=\sum_{j=0}^{p-2}\mu_{j+1}\varphi(x^{j})=\sum_{j=0}^{p-2}\mu_{j+1}X^{j}.$

Since $X$ is circulant, $X^{j}$ is also circulant for each $j=0,\dots,p-2$.
This implies that $\varphi(f)$ is a circulant matrix.

Conversely, we observe that each circulant matrix $C$ can be uniquely written
as $C=\sum_{j=0}^{p-2}c_{j+1}X^{j}$ with $c_{i}\in\mathbb{Q},j=1,\dots,p-1$,
since $\\{X^{j}\\}_{j=0}^{p-2}$ is a basis of the space of all
$(p-1)\times(p-1)$ circulant matrices over $\mathbb{Q}$. This implies that
$\varphi^{-1}(C)=\sum_{j=0}^{p-2}c_{j+1}x^{j}$ is a skew polynomial whose
coefficients are in $\mathbb{Q}$. ∎

For each $0\leq i\leq p-2$, we define a $(p-1)$-dimensional subspace of
matrices:

$\mathcal{L}_{i}\coloneqq\operatorname{span}\left\\{X^{i}Y^{j}:0\leq j\leq
p-2\right\\}.$

######  Definition 3.10.

Given a matrix $A\in\mathbb{Q}^{(p-1)\times(p-1)}$, the _skew-sparsity_ of $A$
is the smallest positive integer $s$ such that $A\in\bigoplus_{i\in
I}\mathcal{L}_{i}$ for some $I\subseteq\\{0,\dots,p-2\\}$ with $|I|=s$.
Equivalently, the skew-sparsity of $A$ can also be defined as the sparsity of
$\varphi^{-1}(A)$, where $\varphi$ is the map defined in (5).

As a direct consequence of Proposition 3.6, we have the following:

######  Theorem 3.11.

Let $p$ be a prime number and let $A,B$ be $(p-1)\times(p-1)$ matrices over
$\mathbb{Q}$. If $A\in\bigoplus_{i\in I}\mathcal{L}_{i}$ and
$B\in\bigoplus_{k\in K}\mathcal{L}_{k}$ for some subsets
$I,K\subseteq\\{0,\dots,p-2\\}$, then the product $AB$ can be computed by an
algorithm with complexity $\operatorname{O}^{\sim}(T^{\omega-2}p^{2})$ where
$T$ is the cardinality of $I+K\pmod{p-1}$. In particular, if there exists some
$\epsilon>0$ such that

$|I||K|\leq O(p^{1-\epsilon/(\omega-2)}),$

then the product $AB$ can be computed by an algorithm with complexity
$\operatorname{O}(p^{\omega-\epsilon})$.

###  3.4. Structure of $\mathcal{L}_{i}$

This subsection is devoted to the discussion of the structure of
$\mathcal{L}_{i}$. By definition, we have

$\mathcal{L}_{i}=X^{i}\mathcal{L}_{0},\quad i=0,\dots,p-2.$

Observing that for any
$A=\begin{bmatrix}a_{1},a_{2},\dots,a_{p-1}\end{bmatrix}^{\scriptscriptstyle\mathsf{T}}\in\mathbb{Q}^{(p-1)\times(p-1)}$,
we have

$XA=\begin{bmatrix}a_{2},\dots,a_{p-1},a_{1}\end{bmatrix}^{\scriptscriptstyle\mathsf{T}}.$

This implies that the effect of left multiplication by $X$ is the cyclic shift
up of rows. Thus it suffices to consider the structure of $\mathcal{L}_{0}$,
which consists of $(p-1)\times(p-1)$ matrices over $\mathbb{Q}$ corresponding
to elements in $\mathbb{Q}(\beta)\subseteq\mathbb{Q}(\beta)[x,\sigma]$, via
the isomorphism $\varphi$ defined in (5).

According to (8), the matrix $Y$ satisfies the relation

$\begin{bmatrix}v_{1}\\\ \vdots\\\ v_{p-1}\\\
\end{bmatrix}\beta=Y\begin{bmatrix}v_{1}\\\ \vdots\\\ v_{p-1}\\\
\end{bmatrix}.$

By induction, we may further derive

$\begin{bmatrix}v_{1}\\\ \vdots\\\ v_{p-1}\\\
\end{bmatrix}\beta^{i}=Y^{i}\begin{bmatrix}v_{1}\\\ \vdots\\\ v_{p-1}\\\
\end{bmatrix},\quad i=0,\dots,p-2.$ (11)

With (11), we are able to characterize $\mathcal{L}_{0}$. To that end, we
notice that an element $a\in\mathbb{Q}(\beta)$ can be written as

$a=c_{1}\beta+c_{2}\beta^{2}+\cdots+c_{p-1}\beta^{p-1}$

for some $c_{1},\dots,c_{p-1}\in\mathbb{Q}$. We introduce two permutations
$s,q\in\mathfrak{S}_{p-1}$ of $\\{1,\dots,p-1\\}$ defined by

$\displaystyle r^{s(i)-1}+i$ $\displaystyle\equiv 0\pmod{p},\quad
i=1,\dots,p-1,$ $\displaystyle r^{q(i)-1}-i$ $\displaystyle\equiv
0\pmod{p},\quad i=1,\dots,p-1.$

######  Lemma 3.12.

For any $1\leq i\leq p-1$, the $s(i)$-th row of $Y^{j}$ is

$(Y^{j})_{s(i)}=\left\\{\begin{aligned} &E_{q(j-i)},\quad if\quad j-i>0\\\
&E_{q(p+j-i)},\quad if\quad j-i<0\\\ &J,\quad\quad\quad if\quad
j-i=0\end{aligned}\right.$

where $J=(-1,\dots,-1)$ and $E_{k}$ is the row vector whose entries are all
zero except the $k$-th, which is $1$.

###### Proof.

By definition of $\varphi$, elements in the $s(i)$-th row of $Y^{j}$ are
$\mathbb{Q}$-coefficients in the expansion of
$v_{s(i)}\beta^{j}\in\mathbb{Q}(\beta)$ with respect to the basis
$\\{v_{1},\dots,v_{p-1}\\}$. Since $v_{s(i)}\beta^{j}=\beta^{r^{s(i)-1}+j}$
and $r^{s(i)-1}\equiv-i\pmod{p}$, we have $v_{s(i)}\beta^{j}=\beta^{j-i}$. We
conclude the argument by the following three cases:

  * •

if $j-i>0$, then $\beta^{j-i}=v_{q(j-i)}$, thus $(Y^{j})_{s(i)}=E_{q(j-i)}$;

  * •

if $j-i<0$, then $\beta^{j-i}=\beta^{p+j-i}=v_{q(p+j-i)}$ and hence
$(Y^{j})_{s(i)}=E_{q(p+j-i)}$;

  * •

if $j-i=0$, then $\beta^{0}=-v_{1}-\cdots-v_{p-1}$ and $(Y^{j})_{s(i)}=J$.

∎

######  Corollary 3.13.

For any $1\leq i\leq p-1$ and $c_{1},\dots,c_{p-1}\in\mathbb{Q}$, the
$s(i)$-th row of $\varphi(\sum_{j=1}^{p-1}c_{j}\beta^{j})$ is

$\left(\sum_{j=1}^{i-1}c_{j}E_{q(p+j-i)}\right)+c_{i}J+\left(\sum_{j=i+1}^{p-1}c_{j}E_{q(j-i)}\right).$

###### Proof.

By the linearity of $\varphi$, we have
$\varphi(\sum_{j=1}^{p-1}c_{j}\beta^{j})=\sum_{j=1}^{p-1}c_{j}Y^{j}$ and the
conclusion follows directly from Lemma 3.12. ∎

Next we consider

$A=\begin{bmatrix}E_{s(1)}\\\ E_{s(2)}\\\ \vdots\\\
E_{s(n)}\end{bmatrix},\quad B=\begin{bmatrix}E_{q(1)}\\\ E_{q(2)}\\\ \vdots\\\
E_{q(n)}\end{bmatrix}^{\scriptscriptstyle\mathsf{T}}.$ (12)

Since $r^{s(i)-1}+i\equiv 0\pmod{p}$ and $r^{q(i)-1}\equiv i\pmod{p}$, we have

$r^{s(i)}+r^{q(i)}=0\pmod{p},$

We notice that $r$ is a primitive root of $\mathbb{Z}_{p}$, thus
$r^{\frac{p-1}{2}}\equiv-1\pmod{p}$ and

$s(i)\equiv q(i)+\frac{p-1}{2}\pmod{p-1}.$ (13)

According to (13), we obtain the next lemma describing a relation between
matrices $A$ and $B$.

######  Lemma 3.14.

Let $A,B$ be as defined in (12). We have $AB=\begin{bmatrix}0&0&\cdots&1\\\
0&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 1&0&\cdots&0\\\ \end{bmatrix}$.

###### Proof.

By definition, we have $AB=(c_{ij})_{(p-1)\times(p-1)}$ where
$c_{ij}=\delta_{s(i)q(j)},1\leq i,j\leq p-1$ and $\delta_{kl}$ is the
Kronecker delta. Equation (13) implies

$c_{ij}=\delta_{q(i)+\frac{p-1}{2},q(j)},$

where $q(i)+\frac{p-1}{2}$ is regarded as an element in $\mathbb{Z}_{p-1}$. We
observe that $q(i)+\frac{p-1}{2}\equiv q(j)\pmod{p-1}$ if and only if
$r^{q(i)+\frac{p-1}{2}}\equiv r^{q(j)}\pmod{p}$. Since
$r^{\frac{p-1}{2}}\equiv-1\pmod{p}$ and $r^{q(k)}\equiv k\pmod{p}$ for any
$1\leq k\leq p-1$, we may further conclude that $q(i)+\frac{p-1}{2}\equiv
q(j)\pmod{p-1}$ if and only if $-i\equiv j\pmod{p}$. Thus we obtain

$c_{ij}=\begin{cases}1,\quad\text{if $p$ divides $(i+j)$},\\\
0,\quad\text{otherwise}.\end{cases}$

∎

Given $c_{1},\dots,c_{p-1}\in\mathbb{Q}$, we define

$P(c_{1},\dots,c_{p-1})\coloneqq\begin{bmatrix}0&c_{p-1}&\cdots&c_{3}&c_{2}\\\
c_{1}&0&\cdots&c_{4}&c_{3}\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
c_{p-3}&c_{p-4}&\cdots&0&c_{p-1}\\\ c_{p-2}&c_{p-3}&\cdots&c_{1}&0\\\
\end{bmatrix},\quad
Q(c_{1},\dots,c_{p-1})\coloneqq\begin{bmatrix}c_{1}&c_{1}&\cdots&c_{1}\\\
c_{2}&c_{2}&\cdots&c_{2}\\\ \vdots&\vdots&\ddots&\vdots\\\
c_{p-2}&c_{p-2}&\cdots&c_{p-2}\\\ c_{p-1}&c_{p-1}&\cdots&c_{p-1}\\\
\end{bmatrix}.$

For instance, if $p=7$ then we have

$P(c_{1},\dots,c_{6})=\begin{bmatrix}0&c_{6}&c_{5}&c_{4}&c_{3}&c_{2}\\\
c_{1}&0&c_{6}&c_{5}&c_{4}&c_{3}\\\ c_{2}&c_{1}&0&c_{6}&c_{5}&c_{4}\\\
c_{3}&c_{2}&c_{1}&0&c_{6}&c_{5}\\\ c_{4}&c_{3}&c_{2}&c_{1}&0&c_{6}\\\
c_{5}&c_{4}&c_{3}&c_{2}&c_{1}&0\\\ \end{bmatrix},\quad
Q(c_{1},\dots,c_{6})=\begin{bmatrix}c_{1}&c_{1}&c_{1}&c_{1}&c_{1}&c_{1}\\\
c_{2}&c_{2}&c_{2}&c_{2}&c_{2}&c_{2}\\\ c_{3}&c_{3}&c_{3}&c_{3}&c_{3}&c_{3}\\\
c_{4}&c_{4}&c_{4}&c_{4}&c_{4}&c_{4}\\\ c_{5}&c_{5}&c_{5}&c_{5}&c_{5}&c_{5}\\\
c_{6}&c_{6}&c_{6}&c_{6}&c_{6}&c_{6}\end{bmatrix}.$

######  Theorem 3.15.

The subspace $\mathcal{L}_{0}\subseteq\mathbb{Q}^{(p-1)\times(p-1)}$ consists
of matrices of the form

$A^{-1}\left(P(c_{1},\dots,c_{p-1})-Q(c_{1},\dots,c_{p-1})\right)A,$

where $c_{1},\dots,c_{p-1}\in\mathbb{Q}$. Moreover, for each $0\leq i\leq
p-2$, the subspace $\mathcal{L}_{i}$ is obtained by shifting up rows of
elements in $\mathcal{L}_{0}$ by $i$. Here for a given
$C=(c_{kl})_{k,l=1}^{p-1}\in\mathbb{Q}^{(p-1)\times(p-1)}$, the matrix
obtained by shifting up rows of $C$ by $i$ is simply
$(c_{k-i,j})_{k,l=1}^{p-1}$ and $(k-i)$ should be understood as
$(k-i)\pmod{p-1}$.

###### Proof.

Let $A,B$ be matrices defined in (12). Then for any
$c_{1},\dots,c_{p-1}\in\mathbb{Q}$, it is straightforward to verify the
equation

$A\varphi\left(\sum_{j=1}^{p-1}c_{j}\beta^{j}\right)B=\left(P(c_{1},\dots,c_{p-1})-Q(c_{1},\dots,c_{p-1})\right)\begin{bmatrix}0&0&\cdots&1\\\
0&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 1&0&\cdots&0\\\ \end{bmatrix}.$

According to Lemma 3.14, we obtain the desired characterization of
$\mathcal{L}_{0}$. ∎

We conclude this section by a remark on the multiplication of matrices in
$\mathcal{L}_{0}$. By Theorem 3.15, an element in $\mathcal{L}_{0}$ can be
written as $A(P(c_{1},\dots,c_{p-1})-Q(c_{1},\dots,c_{p-1}))A^{-1}$ for some
$c_{1},\dots,c_{p-1}\in\mathbb{Q}$. On the one hand, we notice that in
particular, $P(c_{1},\dots,c_{p-1})$ is a Toeplitz matrix and
$Q(c_{1},\dots,c_{p-1})$ is a rank one matrix. We recall that the product of
any matrix with a rank one matrix can be evaluated by
$\operatorname{O}(p^{2})$ operations, while the exsiting optimal algorithm for
the product of a matrix and a Toeplitz matrix has complexity
$\operatorname{O}^{\sim}(p^{2})$ [17]. Therefore, the product of two matrices
in $\mathcal{L}_{0}$ can be evaluated by $\operatorname{O}^{\sim}(p^{2})$
operations. Unfortunately, this method can not efficiently multiply elements
in $\bigoplus_{i\in I}\mathcal{L}_{i}$. On the other hand, according to
Theorem 3.11, not only does Algorithm 2 has complexity
$\operatorname{O}^{\sim}(p^{2})$ as well when applied to elements in
$\mathcal{L}_{0}$, but it is also efficient on $\bigoplus_{i\in
I}\mathcal{L}_{i}$ for any $I\subseteq\\{0,\dots,p-2\\}$.

##  4\. randomized matrix multiplication via skew polynomials

According to the proof of Proposition 3.6, we notice that the complexity of
Algorithm 2 is dominated by the cost of multiplying two skew polynomials. Here
the support of the product is estimated by the sumset, which is guaranteed by
Lemma 2.3 (i). However, this only supplies a rough estimate in the occurrence
of cancellations of coefficients. For instance, if we take $f(x)=1-x$ and
$g(x)=\sum_{j=0}^{k}x^{j}$ where $0\leq k\leq p-3$, then obviously we have
$\\#(f\ast_{\sigma}g)=2$ while $\\#\mathbb{S}(f,g)=k+2$. This example also
indicates that the sparsity of the product can be much smaller than the
cardinality of the sumset.

The goal of this section is to present a randomized algorithm for matrix
multiplication. In comparison to Algorithm 2, this randomized algorithm gets
an acceleration by multiplying two skew-polynomials with the sparsity of their
product known a priori. Essentially the randomized algorithm consists of two
parts. One is the estimation of the sparsity of the product and the other is
the interpolation of the product. It is worth to remark that the randomness in
this algorithm is caused by the probabilistic method we use to estimate the
sparsity of the product.

###  4.1. testing $AB\stackrel{{\scriptstyle?}}{{=}}C$

In our randomized algorithm, we need to verify whether our result is correct.
To do that, it suffices to test whether $AB=C$ for arbitrary matrices $A,B$
and $C$. We recall the celebrated Freivalds’ algorithm [13] in Algorithm 3.

Algorithm 3 Freivalds’ algorithm

1:matrices $M,A,B\in\mathbb{Q}^{n\times n}$ and $\mu\in(0,1)$.

2:“equal” or “not equal”.

3: compute $k:=\lceil\log_{2}\frac{1}{\mu}\rceil$.

4:for $i=1,\dots,k$ do

5: randomly and uniformly choose $y=(y_{1},\dots,y_{n})\in\\{0,1\\}^{n}$.

6: compute $\xi\coloneqq My$ and $\eta\coloneqq ABy$.

7: if $\xi\neq\eta$, then return “not equal”.

8:end for

9:return “equal”.

By a straightforward calculation one can obtain the following error analysis.

######  Lemma 4.1.

[10, Theorem 1.4] Let $M,A,B\in\mathbb{Q}^{n\times n}$ be matrices. Assume
$M\neq AB$. If $y\in\\{0,1\\}^{n}$ is uniformly and randomly chosen, then

$\mathbf{Pr}(My\neq ABy)\geq\frac{1}{2}.$

######  Theorem 4.2.

If $M=AB$, Algorithm 3 returns “equal”. If $M\neq AB$, it returns “not equal”
with probability at least $(1-\mu)$, and it returns “equal” with probability
at most $\mu$. Assuming that we may obtain a random bit with bit-cost
$\operatorname{O}(1)$, the cost of Algorithm 3 is
$\operatorname{O}(n^{2}\log\frac{1}{\mu})$ over $\mathbb{Q}$.

###### Proof.

If $M=AB$, then $\xi=\eta$ and thus Algorithm 3 returns “equal”. If $M\neq
AB$, we observe that Algorithm 3 returns “not equal” if and only if in at
least one iteration of Step 4, $\xi\neq\eta$. In each iteration of Step 4,
Lemma 4.1 indicates that the probability of $\xi=\eta$ is at most $1/2$. Since
all iterations are independent, the probability for Algorithm 3 to return “not
equal” is thus at least

$1-(\frac{1}{2})^{k}=1-(\frac{1}{2})^{\lceil\log_{2}\frac{1}{\mu}\rceil}\geq
1-(\frac{1}{2})^{\log_{2}\frac{1}{\mu}}=1-\mu.$

Next we analyse the complexity. By assumption, the cost of obtaining $kn$
random bits is $\operatorname{O}(kn)$, while computing $\xi$ and $\eta$ in
each iteration of Step 4 costs $\operatorname{O}(n^{2})$ operations over
$\mathbb{Q}$. Since $k=\lceil\log_{2}\frac{1}{\mu}\rceil$, we may conclude
that the total complexity is $\operatorname{O}(n^{2}\log\frac{1}{\mu})$ over
$\mathbb{Q}$. ∎

###  4.2. interpolation

In this subsection, we discuss the interpolation of a skew polynomial
$f\in\mathbb{Q}(\beta)[x;\sigma]/(x^{p-1}-1)$ when an upper bound of its
sparsity is given. Let $T\geq t\coloneqq\\#(f)$ be a positive integer. We want
to recover $f$ from its evaluations $f(v^{\ell}_{1}),\ell=0,1\dots,2T-1$.

Since the sparsity of $f$ is at most $T$, we can write
$f(x)=\sum_{i=1}^{T}c_{i}x^{e_{i}}$. It is clear that there are exactly $t$
nonzero coefficients in $f$. We consider an auxiliary polynomial
$\Lambda_{T}(z)$ defined as follows:

$\Lambda_{T}(z)\coloneqq\prod_{i=1}^{T}(z-v_{e_{j}+1})=z^{T}+\lambda_{T-1}z^{T-1}+\dots+\lambda_{1}z+\lambda_{0}.$
(14)

For any $\ell=0,\dots,2T-1$, we have

$a_{\ell}\coloneqq
f(v^{\ell}_{1})=\sum_{i=1}^{T}c_{i}\sigma^{e_{i}}(v^{\ell}_{1})=\sum_{i=1}^{T}c_{i}v^{\ell}_{e_{i}+1}.$

We observe that for each $j=0,\dots,T-1$, we have

$\displaystyle\sum_{i=1}^{T}c_{i}v_{e_{i}+1}^{j}\Lambda_{T}(v_{e_{i}+1})$
$\displaystyle=\left(\sum_{k=0}^{T-1}\lambda_{k}(c_{1}v_{e_{1}+1}^{k+j}+c_{2}v_{e_{2}+1}^{k+j}+\dots+c_{T}v_{e_{T}+1}^{k+j})\right)+(c_{1}v_{e_{1}+1}^{T+j}+c_{2}v_{e_{2}+1}^{T+j}+\dots+c_{T}v_{e_{T}+1}^{T+j})$
$\displaystyle=\left(\sum_{k=0}^{T-1}a_{j+k}\lambda_{k}\right)+a_{j+T}.$

According to (14), we have $\Lambda_{T}(v_{e_{i}+1})=0$ for each $1\leq i\leq
T$. This implies that

$\left(\sum_{k=0}^{T-1}a_{j+k}\lambda_{k}\right)+a_{j+T}=0,\quad 0\leq j\leq
T-1.$

We now have the Toeplitz system $A_{T}\vec{\lambda}=\vec{b}$ where

$A_{T}=\begin{bmatrix}a_{T-1}&a_{T}&\cdots&a_{2T-2}\\\
a_{T-2}&a_{T-1}&\cdots&a_{2T-3}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{0}&a_{1}&\cdots&a_{T-1}\end{bmatrix},\quad\vec{\lambda}=\begin{bmatrix}\lambda_{0}\\\
\lambda_{1}\\\ \vdots\\\ \lambda_{T-1}\\\
\end{bmatrix},\quad\vec{b}=-\begin{bmatrix}a_{2T-1}\\\ a_{2T-2}\\\ \vdots\\\
a_{T}\\\ \end{bmatrix}.$ (15)

The matrix $A_{T}$ has rank $t$ which can be seen from the factorization:

$A_{T}=\begin{bmatrix}v_{e_{1}+1}^{T-1}&v_{e_{2}+1}^{T-1}&\cdots&v_{e_{T}+1}^{T-1}\\\
v_{e_{1}+1}^{T-2}&v_{e_{2}+1}^{T-2}&\cdots&v_{e_{T}+1}^{T-2}\\\
\vdots&\vdots&\ddots&\vdots\\\
1&1&\cdots&1\end{bmatrix}\begin{bmatrix}c_{1}&0&\cdots&0\\\
0&c_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&c_{T}\end{bmatrix}\begin{bmatrix}1&v_{e_{1}+1}&\cdots&v_{e_{1}+1}^{T-1}\\\
1&v_{e_{2}+1}&\cdots&v_{e_{2}+1}^{T-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
1&v_{e_{T}+1}&\cdots&v_{e_{T}+1}^{T-1}\end{bmatrix}.$ (16)

Since the $v_{i}$’s are pairwise distinct, the two Vandermonde matrices in
(16) are nonsingular. By assumption, there are exactly $t$ nonzero $c_{i}$’s,
thus the diagonal matrix in (16) is has rank $t$. In particular, $A_{t}$ is
nonsingular.

By choosing the first $t$ evaluations of $f$, we obtain the following
transposed Vandermonde system $V\vec{c}=\vec{a}$ for the coefficients of $f$,
where

$V=\begin{bmatrix}1&1&\cdots&1\\\
v_{e_{1}+1}&v_{e_{2}+1}&\cdots&v_{e_{t}+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
v_{e_{1}+1}^{t-1}&v_{e_{2}+1}^{t-1}&\cdots&v_{e_{t}+1}^{t-1}\end{bmatrix},\quad\vec{c}=\begin{bmatrix}c_{1}\\\
c_{2}\\\ \vdots\\\ c_{t}\end{bmatrix},\quad\vec{a}=\begin{bmatrix}a_{0}\\\
a_{1}\\\ \vdots\\\ a_{t-1}\end{bmatrix}.$ (17)

Now we are ready to present Algorithm 4 for the interpolation of a skew
polynomial with a given upper bound on its sparsity. Main ingredients of
Algorithm 4 are:

  * •

compute the roots $v_{e_{1}+1},\dots,v_{e_{t}+1}$ of $\Lambda_{t}(z)$ from the
$\\{a_{i}\\}_{i=0}^{2T-1}$;

  * •

determine the exponents $e_{1},\dots,e_{t}$ from
$\\{v_{e_{j}+1}\\}_{j=1}^{t}$;

  * •

compute the coefficients $c_{1},\dots,c_{t}$ by $\\{a_{k}\\}_{k=0}^{t-1}$ and
$\\{v_{e_{j}+1}\\}_{j=1}^{t}$.

Algorithm 4 Interpolation with an upper bound of sparsity

1:evaluations $\\{a_{\ell}=f(v^{\ell}_{1})\\}_{\ell=0}^{2T-1}$ and upper bound
$T\geq\\#(f)$.

2:coefficients of $f$.

3: compute $t=\operatorname{rank}(A_{T})$.

4: form $A_{t},\vec{b}_{t}$ by (15) and solve the Toeplitz system
$A_{t}\vec{\lambda}_{t}=\vec{b}_{t}$ for $\vec{\lambda}_{t}$.

5: construct $\Lambda_{t}(z)$ by (14) and compute the roots
$v_{e_{1}+1},\dots,v_{e_{t}+1}$ of $\Lambda_{t}(z)$ by fast multiple
evaluations at $v_{1},v_{2},\dots,v_{p-1}$.

6: form $V$ and $\vec{a}$ by (17) and solve the transposed Vandermonde system
$V\vec{c}=\vec{a}$ for $\vec{c}$.

######  Lemma 4.3.

Algorithm 4 computes coefficients of $f$ by
$\operatorname{O}^{\thicksim}(p^{2})$ arithmetic operations over $\mathbb{Q}$.

###### Proof.

The correctness follows from the discussion before Algorithm 4, thus it is
sufficient to analyse the complexity. According to [23], the cost of Steps 3,
4 and 6 are respectively $\operatorname{O}^{\thicksim}(T)$ over
$\mathbb{Q}(\beta)$. In Step 5, it requires $\operatorname{O}^{\thicksim}(p)$
operations over $\mathbb{Q}(\beta)$ by [32]. Since $T\leq p-1$, the total
complexity of Algorithm 4 is $\operatorname{O}^{\thicksim}(p^{2})$ over
$\mathbb{Q}$. ∎

###  4.3. Monte Carlo matrix multiplication via skew polynomials

Our randomized algorithm for matrix multiplication is obtained by assembling
ingredients discussed in Subsections 3.1, 4.1 and 4.2. We record our algorithm
in Algorithm 5.

Algorithm 5 Monte Carlo matrix multiplication

1:$A,B\in\mathbb{Q}^{(p-1)\times(p-1)}$ and $\nu\in(0,1)$

2:$AB$.

3: find a $p$-th primitive root of unity $\beta$ and a generator $r$ of
$\mathbb{Z}_{p}$.

4: compute $v_{i}\coloneqq\beta^{r^{i-1}},i=1,\dots,p-1$.

5: compute $\varphi^{-1}(A),\varphi^{-1}(B)$ by (9) w.r.t the normal basis
$\\{v_{1},\dots,v_{p-1}\\}$.

6: set $T=1$ and $\tau=0$.

7:while $\tau=0$ do

8: compute
$a_{\ell}\coloneqq\varphi^{-1}(A)\ast_{\sigma}\varphi^{-1}(B)(v^{i}_{1}),\ell=0,1,\dots,2T-1$.

9: interpolate $f$ with input $(\\{a_{\ell}\\}_{\ell=0}^{2T-1},T)$ by
Algorithm 4.

10: compute $M\coloneqq\varphi(f)$ by Algorithm 1 w.r.t. the normal basis
$\\{v_{1},\dots,v_{p-1}\\}$.

11: test whether $M=AB$ by calling Algorithm 3 with
$\mu=\frac{\nu}{\lceil\log_{2}(p-1)\rceil}$.

12: if the test returns “equal” then set $\tau=1$; if the test returns “not
equal” then set $T\coloneqq 2T$.

13:end while

######  Theorem 4.4.

Given $A,B\in\mathbb{Q}^{(p-1)\times(p-1)}$ and $\nu\in(0,1)$, Algorithm 5
computes $AB$ correctly with probability at least $(1-\nu)$. The cost of
Algorithm 5 is
$\operatorname{O}^{\thicksim}(t^{\omega-2}p^{2}+p^{2}\log\frac{1}{\nu})$ over
$\mathbb{Q}$, where $t$ is the skew-sparsity of $AB$.

###### Proof.

We denote by $t=\\#\varphi^{-1}(AB)$ the skew-sparsity of $AB$. Assume $k$ is
a positive integer such that $2^{k-1}<t\leq 2^{k}$. If Algorithm 5 computes
$AB$ correctly, then according to Lemma 4.3, the while loop in Step 7 of
Algorithm 5 must terminate at the $s$-th iteration, where $s\leq k$. We denote
by $M_{i}$ (resp. $\tau_{i},T_{i}$) the matrix (resp. numbers) obtained in
Step 10 (resp. Step 12) at the $i$-th iteration, $i=1,\dots,s$. Then clearly
we have

  * •

$T_{i}=2^{i},i=1,\dots,s$;

  * •

$M_{i}\neq AB,i=1,\dots,s-1$ and $M_{s}=AB$;

  * •

$\tau_{1}=\cdots=\tau_{s-1}=0$ and $\tau_{s}=1$.

By Theorem 4.2, we obtain

$\mathbf{Pr}\left(M_{i}\neq
AB,\tau_{i}=0,i=1,\dots,s-1\right)\geq\left(1-\frac{\nu}{\lceil\log_{2}(p-1)\rceil}\right)^{s-1}\geq
1-\frac{(s-1)\nu}{\lceil\log_{2}(p-1)\rceil}\geq 1-\nu.$

Therefore, the probability of Algorithm 5 computing $AB$ correctly is at least
$(1-\nu)$.

It remains to analyse the complexity. First of all, if $r$ is a primitive root
of $\mathbb{Z}_{p}$, then $r^{1},r^{2},\dots,r^{p-1}\pmod{p}$ must be distinct
integers from $1$ to $(p-1)$. Since the complexity of checking the membership
of a number in $\\{2,\dots,p-1\\}$ is $\operatorname{O}(p^{2})$ over
$\mathbb{Q}$, the total cost of Steps 3 and 4 is $\operatorname{O}(p^{2})$. By
Proposition 3.4, the complexity of Step 5 is
$\operatorname{O}^{\thicksim}(p^{2})$ over $\mathbb{Q}$.

Next we count the complexity of each iteration in Step 7. According to the
analysis in Section 2.3.2, the complexity of Step 8 is
$\operatorname{O}^{\sim}(t^{\omega-2}p^{2})$ over $\mathbb{Q}$. By Lemma 4.3
and Proposition 3.5, the total cost of Steps 9 and 10 is
$\operatorname{O}^{\thicksim}(p^{2})$ operations in $\mathbb{Q}$. Theorem 4.2
impies that Step 11 requires
$\operatorname{O}^{\thicksim}(n^{2}\log\frac{1}{\nu})$ operations in
$\mathbb{Q}$.

Lastly, since there are only at most $\log_{2}(p-1)$ iterations, the
complexity of Algorithm 5 is
$\operatorname{O}^{\thicksim}(t^{\omega-2}p^{2}+p^{2}\log\frac{1}{\nu})$ over
$\mathbb{Q}$. ∎

##  5\. Conclusion

In this paper, we propose a new method for $(p-1)\times(p-1)$ matrix
multiplication over $\mathbb{Q}$ via skew polynomials, where $p$ is a prime
number. Via the sparsity of skew polynomials, we are able to define skew-
sparse matrices and we obtain an explicit characterization of such matrices.
Based on the new method, we design a deterministic algorithm for matrix
multiplication, which attains an acceleration if the product is skew-sparse.
We also propose a randomized algorithm which can further accelerate the matrix
multiplication.

## References

  * [1] J. Alman and V. V. Williams.  A refined laser method and faster matrix multiplication.  In D. Marx, editor, Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 - 13, 2021, pages 522–539. SIAM, 2021.
  * [2] A. Ambainis, Y. Filmus, and F. Le Gall.  Fast matrix multiplication: Limitations of the coppersmith-winograd method.  In Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing, STOC ’15, page 585?593, New York, NY, USA, 2015. Association for Computing Machinery.
  * [3] D. Bini.  Relations between exact and approximate bilinear algorithms. Applications.  Calcolo, 17(1):87–97, 1980.
  * [4] P. Bürgisser, M. Clausen, and M. A. Shokrollahi.  Algebraic complexity theory.  Ibm Journal of Research and Development, 25(5):825–832, 2010.
  * [5] J. R. Bunch and J. E. Hopcroft.  Triangular factorization and inversion by fast matrix multiplication.  Mathematics of Computation, 28(125):231–236, 1974.
  * [6] J. Carrier, L. Greengard, and V. Rokhlin.  A fast adaptive multipole algorithm for particle simulations.  SIAM J. Sci. Statist. Comput., 9(4):669–686, 1988.
  * [7] X. Caruso and J. Le Borgne.  Fast multiplication for skew polynomials.  In Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation, ISSAC ’17, pages 77–84, New York, NY, USA, 2017. Association for Computing Machinery.
  * [8] L. Chiantini, J. D. Hauenstein, C. Ikenmeyer, J. M. Landsberg, and G. Ottaviani.  Polynomials and the exponent of matrix multiplication.  Bull. Lond. Math. Soc., 50(3):369–389, 2018.
  * [9] D. Coppersmith and S. Winograd.  Matrix multiplication via arithmetic progressions.  J. Symb. Comput., 9(3):251–280, 1990.
  * [10] M. Cryan.  Probability and computing: randomized algorithms and probabilistic analysis [book review of MR2144605].  Bull. Symbolic Logic, 12(2):304–308, 2006.
  * [11] C. De Sa, A. Gu, R. Puttagunta, C. Ré, and A. Rudra.  A two-pronged progress in structured dense matrix vector multiplication.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1060–1079. SIAM, Philadelphia, PA, 2018.
  * [12] M. J. Fischer and A. R. Meyer.  Boolean matrix multiplication and transitive closure.  In 12th Annual Symposium on Switching and Automata Theory (swat 1971), pages 129–131. IEEE, 1971.
  * [13] R. Freivalds.  Probabilistic machines can use less running time.  In Information processing 77 (Proc. IFIP Congr., Toronto, Ont., 1977), pages 839–842. IFIP Congr. Ser., Vol. 7, 1977.
  * [14] F. L. Gall.  Powers of tensors and fast matrix multiplication.  In K. Nabeshima, K. Nagasaka, F. Winkler, and Á. Szántó, editors, International Symposium on Symbolic and Algebraic Computation, ISSAC ’14, Kobe, Japan, July 23-25, 2014, pages 296–303. ACM, 2014.
  * [15] S. Gao.  Normal bases over finite fields. University of Waterloo Waterloo, Canada, 1993.
  * [16] M. Giesbrecht, Q.-L. Huang, and É. Schost.  Sparse multiplication for skew polynomials.  In Proceedings of the 45th International Symposium on Symbolic and Algebraic Computation, pages 194–201, 2020.
  * [17] G. H. Golub and C. F. Van Loan.  Matrix computations.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, fourth edition, 2013.
  * [18] L. Greengard and V. Rokhlin.  A fast algorithm for particle simulations [ MR0918448 (88k:82007)].  volume 135, pages 279–292. 1997.  With an introduction by John A. Board, Jr., Commemoration of the 30th anniversary {of J. Comput. Phys.}.
  * [19] F. G. Gustavson.  Two fast algorithms for sparse matrices: Multiplication and permuted transposition.  ACM Trans. Math. Softw., 4(3):250–269, 1978.
  * [20] J. E. Hopcroft and L. R. Kerr.  On minimizing the number of multiplications necessary for matrix multiplication.  SIAM J. Appl. Math., 20:30–36, 1971.
  * [21] Ian and Munro.  Efficient determination of the transitive closure of a directed graph.  Information Processing Letters, 1971.
  * [22] C. E. Kaltofen.  On fast multiplication of polynomials over arbitrary algebras.  Acta Informatica, 1991.
  * [23] E. Kaltofen and L. Yagati.  Improved sparse multivariate polynomial interpolation algorithms.  In Symbolic and Algebraic Computation, International Symposium ISSAC’88, Rome, Italy, July 4-8, 1988, Proceedings, pages 467–474. Springer, 1988.
  * [24] W. Keller-Gehrig.  Fast algorithms for the characteristics polynomial.  Theoretical computer science, 36:309–317, 1985.
  * [25] L.-H. Lim and K. Ye.  Ubiquity of the exponent of matrix multiplication.  In Proceedings of the 45th International Symposium on Symbolic and Algebraic Computation, ISSAC ’20, page 8?11, New York, NY, USA, 2020. Association for Computing Machinery.
  * [26] V. Olshevsky and A. Shokrollahi.  Matrix-vector product for confluent Cauchy-like matrices with application to confluent rational interpolation.  In Proceedings of the Thirty-Second Annual ACM Symposium on Theory of Computing, pages 573–581. ACM, New York, 2000.
  * [27] V. Y. Pan.  Strassen’s algorithm is not optimal: Trililnear technique of aggregating, uniting and canceling for constructing fast algorithms for matrix operations.  In 19th Annual Symposium on Foundations of Computer Science, Ann Arbor, Michigan, USA, 16-18 October 1978, pages 166–176. IEEE Computer Society, 1978.
  * [28] A. Schorr.  Fast algorithm for sparse matrix multiplication.  Inf. Process. Lett., 15(2):87–89, 1982.
  * [29] V. Strassen.  Gaussian elimination is not optimal.  Numerische Mathematik, 13(4):354–356, 1969.
  * [30] V. Strassen.  The asymptotic spectrum of tensors and the exponent of matrix multiplication.  In 27th Annual Symposium on Foundations of Computer Science, Toronto, Canada, 27-29 October 1986, pages 49–54. IEEE Computer Society, 1986\.
  * [31] C. F. Van Loan.  The ubiquitous Kronecker product.  volume 123, pages 85–100. 2000.  Numerical analysis 2000, Vol. III. Linear algebra.
  * [32] J. Von Zur Gathen and J. Gerhard.  Modern computer algebra.  Cambridge university press, 1999.
  * [33] D. Wiedemann.  Solving sparse linear equations over finite fields.  IEEE transactions on information theory, 32(1):54–62, 1986.
  * [34] V. V. Williams.  Multiplying matrices faster than coppersmith-winograd.  In H. J. Karloff and T. Pitassi, editors, Proceedings of the 44th Symposium on Theory of Computing Conference, STOC 2012, New York, NY, USA, May 19 - 22, 2012, pages 887–898. ACM, 2012.
  * [35] K. Ye and L.-H. Lim.  Fast structured matrix computations: tensor rank and Cohn-Umans method.  Found. Comput. Math., 18(1):45–95, 2018.
  * [36] R. Yuster and U. Zwick.  Fast sparse matrix multiplication.  ACM Trans. Algorithms, 1(1):2–13, 2005.