Nearly Optimal Sparse Polynomial Multiplication

TL;DR

This paper presents a new nearly optimal algorithm for multiplying sparse polynomials efficiently, improving upon previous methods by handling cancellations and structural sparsity more effectively.

Contribution

The paper introduces a clean, nearly optimal algorithm for sparse polynomial multiplication that advances previous work by addressing cancellations and structural sparsity.

Findings

Algorithm runs in near-optimal time proportional to input and output size

Handles cancellations in polynomial coefficients effectively

Improves efficiency over previous algorithms in structural sparsity cases

Abstract

In the sparse polynomial multiplication problem, one is asked to multiply two sparse polynomials f and g in time that is proportional to the size of the input plus the size of the output. The polynomials are given via lists of their coefficients F and G, respectively. Cole and Hariharan (STOC 02) have given a nearly optimal algorithm when the coefficients are positive, and Arnold and Roche (ISSAC 15) devised an algorithm running in time proportional to the "structural sparsity" of the product, i.e. the set supp(F)+supp(G). The latter algorithm is particularly efficient when there not "too many cancellations" of coefficients in the product. In this work we give a clean, nearly optimal algorithm for the sparse polynomial multiplication problem.

Equations

$$(x\star y{)}_i=\sum _{j,j^{\prime }\in [n]:(j+j^{\prime })\text{leavevmode}\mathrm{mod}\text{leavevmode}n=i}{x}_j{y}_{j^{\prime }}.$$

$$k{\log }^{2}N\cdot \log (k{\log }^{2}N)\cdot \log \log k+$$

$$s\cdot \log k\log N\cdot (\log \log k{)}^2+$$

$$\tilde{O}({\log }^{3}N)\cdot \log \Vert x\star y{\Vert }_0.$$

$$\sum _{j\in [N]:j\text{leavevmode}\mathrm{mod}\text{leavevmode}p=i}(x\star y{)}_j\cdot g^j.$$

$$h_m(i)=i\text{leavevmode}\mathrm{mod}\text{leavevmode}m.$$

$${\left({\mathcal{P}}_m(x)\right)}_i=\sum _{j\in [N]:h_m(j)=i}{x}_j\cdot {\omega }^j,\forall i\in [m].$$

$$O(\left(\Vert x{\Vert }_0+\Vert y{\Vert }_0+\Vert w{\Vert }_0\right)\log N+m\log m).$$

$${\mathcal{P}}_m((x\star y)-w)={\mathcal{P}}_m(x\star y)-{\mathcal{P}}_m(w),$$

$$|{\left({\mathcal{P}}_m(x\star y)\right)}_i-{\left({\mathcal{P}}_m(x)\star {\mathcal{P}}_m(y)\right)}_i|\le 1/\mathrm{poly}(N).$$

$${\left({\mathcal{P}}_m(x\star y)\right)}_i={\left({\mathcal{P}}_m(x)\star {\mathcal{P}}_m(y)\right)}_i.$$

$${\left({\mathcal{P}}_m(x)\star {\mathcal{P}}_m(y)\right)}_i=$$

$$\sum _{\ell ,{\ell }^{\prime }\in [m]:(\ell +{\ell }^{\prime })\text{leavevmode}\mathrm{mod}\text{leavevmode}m=i}({\mathcal{P}}_m(x){)}_{\ell }\cdot ({\mathcal{P}}_m(y){)}_{{\ell }^{\prime }}=$$

$$\sum _{\ell ,{\ell }^{\prime }\in [m]:(\ell +{\ell }^{\prime })\text{leavevmode}\mathrm{mod}\text{leavevmode}m=i}\left(\sum _{j,j^{\prime }\in [n]:h_m(j)=\ell ,h_m(j^{\prime })={\ell }^{\prime }}{x}_j{y}_{j^{\prime }}{\omega }^{j+j^{\prime }}\right)=$$

$$\sum _{\ell ,{\ell }^{\prime }\in [m],j,j^{\prime }\in [N]:h_m(j)=\ell ,h_m(j^{\prime })={\ell }^{\prime },(\ell +{\ell }^{\prime })\text{leavevmode}\mathrm{mod}\text{leavevmode}m=i}{x}_j{y}_{j^{\prime }}{\omega }^{j+j^{\prime }}=$$

$$\sum _{j,j^{\prime }\in [N]:(h_m(j)+h_m(j^{\prime }))\text{leavevmode}\mathrm{mod}\text{leavevmode}m=i}{x}_j{y}_{j^{\prime }}{\omega }^{j+j^{\prime }}=$$

$$\sum _{j,j^{\prime }\in [N]:h_m((j+j^{\prime })\text{leavevmode}\mathrm{mod}\text{leavevmode}m))=i}{x}_j{y}_{j^{\prime }}{\omega }^{j+j^{\prime }}=$$

$$\sum _{j^{\prime \prime }\in [N],h_m(j^{\prime \prime }\text{leavevmode}\mathrm{mod}\text{leavevmode}m)=i}\sum _{j,j^{\prime }\in [N]:j+j^{\prime }=j^{\prime \prime }}{x}_j{y}_{j^{\prime }}{\omega }^{j^{\prime \prime }}=$$

$$\sum _{j^{\prime \prime }\in [N],h_m(j^{\prime \prime }\text{leavevmode}\mathrm{mod}\text{leavevmode}m)=i}{\omega }^{j^{\prime \prime }}\left(\sum _{j,j^{\prime }\in [N]:j+j^{\prime }=j^{\prime \prime }}{x}_j{y}_{j^{\prime }}\right)=$$

$$\sum _{j^{\prime \prime }\in [N],h_m(j^{\prime \prime })=i}{\omega }^{j^{\prime \prime }}{\left(x\star y\right)}_{j^{\prime \prime }},$$

$$\left\{0,\dots ,\lceil N\right\},$$

$$\left\{\lceil \frac{N}{2}\rceil +1,\dots ,N\right\}$$

$$\left\{N+1,\dots ,\lceil \frac{3}{2}N\rceil \right\}$$

$$\left\{\lceil \frac{3}{2}N\rceil +1,\dots ,2N-1\right\}$$

$$\forall j^{\prime }\in \mathrm{supp}((x\star y)-w)\setminus \{j\}:h_p(j^{\prime })\ne h_p(j),$$

$$\mathbb{P}\left[h_p(j)=h_p(j^{\prime })\right]\le 1/B.$$

$${\mathcal{E}}_j:\{\exists j^{\prime }\in \mathrm{supp}((x\star y)-w)\setminus \{j\}:h_p(j)=h_p(j^{\prime })\}.$$

$$\mathbb{E}\left[\sum _{j\in \mathrm{supp}((x\star y)-w)}{X}_j\right]\le (2C^{-2})\Vert (x\star y)-w{\Vert }_0.$$

$$\Vert z-((x\star y)-w){\Vert }_0\le (5\gamma )\Vert (x\star y)-w){\Vert }_0.$$

$$({\mathcal{P}}_p(x\star y){)}_i-({\mathcal{P}}_p(w){)}_i=({\mathcal{P}}_p(x\star y-w){)}_i=$$

$$\sum _{j\in [N]:h_p(j)=i}((x\star y)-w{)}_j{\omega }^j=$$

$$\sum _{j\in [N]:h_p(j)=i\text{leavevmode}\mathrm{and}\text{leavevmode}((x\star y)-w{)}_j\ne 0}((x\star y)-w{)}_j{\omega }^j$$

$$({\mathcal{P}}_p(x\star y){)}_{i^{\ast }}-({\mathcal{P}}_p(w){)}_{i^{\ast }}=((x\star y{)}_{j^{\ast }}-w_{j^{\ast }}){\omega }^{j^{\ast }}.$$

$$\left(\genfrac{}{}{0ex}{}{5\log (1/\delta )}{(5/2)\log (1/\delta )}\right)p^{(5/2)\log (1/\delta )}\le 2^{5\log (1/\delta )}{q}^{(5/2)\log (1/\delta )}\le \delta ,$$

$$\gamma \cdot 4\log (1/\delta )\cdot \Vert (x\star y)-w{\Vert }_0+\log (1/\delta )\Vert (x\star y)-w{\Vert }_0.$$

$$1+\beta \cdot \frac{1}{4}\cdot 5\log (1/\delta )\Vert (x\star y)-w{\Vert }_0>$$

$$\gamma \cdot 4\log (1/\delta )\cdot \Vert (x\star y)-w{\Vert }_0+\log (1/\delta )\Vert (x\star y)-w{\Vert }_0$$

$$(B{\log }^{2}N\cdot \log (B{\log }^{2}N)+$$

Full text

Nearly Optimal Sparse Polynomial Multiplication

Vasileios Nakos Saarland University and Max Planck Institute for Informatics, Saarland Informatics Campus, Saarbrücken, Germany. vnakos@mpi-inf.mpg.de. This work is part of the project TIPEA that has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 850979). Part of this work was done while the author was a ph.D. student at Harvard University and supported by NSF CAREER award CCF-1350670. vnakos@mpi-inf.mpg.de,

Saarland University and Max Planck Institute for Informatics.

In the sparse polynomial multiplication problem, one is asked to multiply two sparse polynomials $f$ and $g$ in time that is proportional to the size of the input plus the size of the output. The polynomials are given via lists, $F$ and $G$, of their coefficients. Cole and Hariharan (STOC 02) have given a nearly optimal algorithm when the coefficients are positive, and Arnold and Roche (ISSAC 15) devised an algorithm running in time proportional to the “structural sparsity” of the product, i.e. the set $\mathrm{supp}(F)+\mathrm{supp}(G)$. The latter algorithm is particularly efficient when there are not ”too many cancellations” of coefficients in the product. In this work we give a clean, nearly optimal algorithm for the sparse polynomial multiplication problem.

1 Introduction

Multiplying two polynomials is a fundamental computational primitive, with multiple applications in computer science. Using the Fast Fourier Transform, one can perform multiplication of polynomials stored as vectors of floating point numbers in time $O(n\log n)$, where $n$ is a bound on the largest degree.

An important and natural question is whether, and under which circumstances, a faster algorithm can be invented. Researchers have tried to obtain algorithms that beat the $O(n\log n)$-time bound, when the two polynomials are sparse, i.e. the number of non-zero terms in each polynomial is at most $k$. Interestingly, some ideas from the relevant literature have found applications in computer algebra packages such as Maple, Mathematica and Singular, including ways to represent and store polynomials [Maz01, MP14, MP15, GR16].

When two polynomials have at most $s$ coefficients, the trivial algorithm gives $O(s^{2}\log n\log s)$ time, which is already and improvement for $s\leq\sqrt{n/\log n}$. With the rise of big data and the exploitation of inherent sparsity in the data sets by researchers from compressed sensing [CRT06b, CRT06a, CT06, HIKP12a] and machine learning, the question of whether two sparse polynomials can be rapidly multiplied becomes a fundamental and natural question. The desirable goal is to obtain an algorithm that is output-sensitive, i.e. runs in nearly linear time with respect to $k+s$, where $k$ is the number of non-zero coefficients in the product. A result of Cole and Hariharan [CH02] obtains an algorithm that runs in $O(k\log^{2}n)$ time, when the coefficients of the two polynomials are non-negative. A data structure for carefully allocating and de-allocating memory has been designed in [Yan98], trying to tackle the problem of memory handling can be the main bottleneck in complexity of sparse multiplication in practical scenarios. The aforementioned algorithm is based on a heap, an idea which also lead to implementations developed in [MP07, MP09, MP14, MP15]. The authors in [MP09] develop a parallel algorithm for multiplying sparse distributed polynomials, where each core uses a heap of pointers to multiply parts of polynomials, exploiting its L3 cache. The same authors in [MP14] have created a data structure suitable for the Maple kernel, that allows for obtains significant performance in many Maple library routines.

When the support of the product is known or structured, work in [Roc08, Roc11, VDHL12, VDHL13] indicates how to perform the multiplication fast. Using techniques from spare interpolation, Arnold and Roche [AR15] have given an algorithm that runs in time that is nearly linear in the “structural sparsity” of the product, i.e. the sumset of the supports of the two polynomials. When there are not “too many” cancellations, this is roughly the same as the size of the support of the product, and the above algorithm is quite efficient. However, in the presence of a considerable amount of cancellations in the product, the aforementioned algorithm becomes sub-optimal. Removing this obstacle seems to be the final step, and has been posed as an open problem in the excellent survey of [Roc18].

In this paper, we resolve the aforementioned open question, giving an algorithm that is near-linear in the size of the input plus the size of the output. This can be considered near-optimal. Moreover, assuming optimality of the Fast Fourier Transform, one cannot expect an algorithm faster than $O((k+s)\log(s+k))$, and our algorithm reaches quite close to this barrier. We note that one can use the big hammers of filter-based (sparse) Fourier sampling [GMS05, HIKP12a, HIKP12b, Kap16, Kap17] along with semi-equispaced Fourier transforms [DR93] to obtain nearly optimal algorithms for sparse polynomial multiplication, but these algorithms are much more complicated and technical, and thus likely to be worse in a real-life scenario. Furthermore, puting all the pieces together is a non-trivial task for someone who is not an expert in Sparse Fourier transform, as how to solve sparse polynomial multiplication using those techniques has not been explicitly written down in the literature. In contrast, our algorithm is much less technically demanding and we feel accessible even to a graduate level student.

2 Preliminaries

We will be concerned with polynomials with integer coefficients. This suffices for most applications, since numbers in a machine are represented using floating point arithmetic. We denote by $\mathbb{Z}_{n}$ the ring of residue modulo $n$ and by $[n]$ the set $\{0,1,\ldots n-1\}$. For a vector $x\in\mathbb{R}^{n}$ we define $\mathrm{supp}(x)=\{i\in[n]:x_{i}\neq 0\}$ and $\|x\|_{0}=|\mathrm{supp}(x)|$. All vectors in this paper are zero-indexed.

We define the cyclical convolution of two vectors $x,y\in\mathbb{R}^{n}$ as the $n$-dimensional vector $x\star y$ such that

[TABLE]

It is easy to see that cyclical convolution is immediately related with polynomial multiplication: if $f(z)=\sum_{j=0}^{n}\alpha_{j}z^{j}$ and $g(z)=\sum_{j=0}^{n}\beta_{j}z^{j}$, we have that the polynomial $(f\cdot g)(z)=\sum_{j=0}^{2n}\gamma_{j}z^{j}$ satisfies $\gamma=\alpha\star\beta$, where $\gamma=(\gamma_{0},\gamma_{1},...\ldots,\gamma_{N})\in\mathbb{Z}^{N},\alpha=(\alpha_{0},\alpha_{1},\ldots,\alpha_{n},0,\ldots,0)\in\mathbb{Z}^{N},\beta=(\beta_{0},\beta_{1},\ldots,\beta_{n},0,\ldots,0)\in\mathbb{Z}^{N}$, for $N=2n$.

Throughout the paper we assume that we work on a machine where the word size is $w=C\log n=\Theta(\log n)$ for some sufficiently large constant $C$, all coefficients of the polynomials fit in a word, and elementary operations between two integers given as part of the input can be done in $O(1)$ time. For a complex number $z$ we denote by $|z|$ its magnitude, and by $\mathrm{arg}(\phi)$ its phase. We also use $\widetilde{O}(f)$ throughtout the paper to denote $f\cdot\mathrm{poly}(\log f)$.

3 Result

The contribution of this work is the following.

Theorem 3.1**.**

Let $u,v\in\mathbb{Z}^{n}$, given as lists of their non-zero coordinates along with their values. Let $N=2n$. Define $x=(u_{0},u_{1},\ldots,u_{n-1},0,\ldots,0)\in\mathbb{Z}^{N}$, as well as $y=(v_{0},v_{1},\ldots,v_{n-1},0,\ldots,0)\in\mathbb{Z}^{N}$. Let $s=\|x\|_{0}+\|y\|_{0}$ (size of input), and $k=\|x\star y\|_{0}+4$ (size of output). Then, with probability $99/100$, we can compute a list which contains the non-zero coefficients and values of $x\star y$, in time $O((s+k)\cdot\mathrm{poly}(\log N))$. In particular, the running time is (ignoring constant factors)

[TABLE]

4 Overview of the Algorithm

Our algorithm resembles the iterative framework employed in the renowned sublinear-time sparse recovery algorithm of [GLPS10]. Similarly to that paper, our algorithm implements a routine which carefully hashes every coordinate $i\in[N]$ to an appropriate number of buckets, and from each bucket identifies the location and values of at least a constant fraction of the non-zero coordinates in $x\star y$. Let $r$ be the $k$-sparse vector obtained. Considering vector $(x\star y)-r$, we can perform a similar “hash and identify” procedure to peel of a constant fraction of the coordinates left. Repeating in $O(\log k)$ phases, we can guarantee that we’ve found all coordinates and cleaned up all mistakes introduced.

The main contribution of our work lies in effectively hashing the support of $x\star y$ and identifying the locations and the values of the non-zero coordinates, given access only to $x$ and $y$. Let us assume for now that $k$ is given to us as a promise. It is a folklore fact (see for example Section 8 in [CL15]) that if one takes a random prime $p\in[10k\log^{2}N]$, folds $x$ by forming a $p$-length vector by summing all $x_{j}$ for $j$ which are equal modulo $p$, folds $y$ in the same way, and then computes the cyclical convolution of the folded vectors, then most coordinates in the support of $x\star y$ are isolated in the resulting vector. However, this trick only gives the values of the non-zero coordinates in $x\star y$, and it is unclear how to extract the corresponding indices. What we observe is that for any (complex) number $g$, the following thing holds. If we first multiply each $x_{j}$ by $g^{j}$ to obtain vector $x^{\prime}$, then fold $x^{\prime}$ as before, do the same thing with $y$, and then compute the cyclical convolution of the folded vectors, we can prove that the vector we obtain has in its $i$th entry the number

[TABLE]

The nice thing now is that if some $j\in\mathrm{supp}(x\star y)$ is isolated under the hash function $h(x)=x\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ p$, we will have $(x\star y)_{j}\cdot g^{j}$ in hand, rather than only $(x\star y)_{j}$ that we had before. From this we would like to infer both $(x\star y)_{j}$ and $j$. It is easy to see that even for $(x\star y)_{j}=1$ this is essentially a discrete logarithm problem. By choosing $g\neq 1$ to be a real number, this would result in $N$-digit numbers, which take too long to manipulate. The right solution is to choose $g=\omega$, where $\omega$ is an $(2\cdot N)$-th root of unity rounded to fit in a word. In this way, every $g^{j}$ can still be written down using $\Theta(\log n)$ digits, and we can infer integer $j\in[N]$ from $g^{j}$ by performing a ternary search over the corresponding quarter of the complex circle. Note that we chose $(2N)$-th root of unity in order to infer both $(x\star y)_{j}$ and $j$ from $(x\ast y)_{j}\cdot g^{j}$: if the phase of the aforementioned complex number is less than $\frac{\pi}{2}$, then $(x\star y)_{j}$ is positive (and vice versa), otherwise it is negative. In the latter case, to read the correct location one needs to decrease the found $j$ by $N$.

Of course, the above discussion assumes that we know $k$. We can guess $k$ by doubling and invoking standard fingerprinting techniques. Putting carefully everything together results in the desired guarantee claimed in 3.1.

5 Algorithm and Proof

We proceed by building the tools needed for the proof of theorem 3.1. In what folows $\omega$ is an $(2N)$-th root of unity; our algorithm should treat it as rounded in order to fit in a word. We remind that our word size is $\Omega(\log n)$, which is the word size needed just to write down the largest exponent of the polynomial.

The following operator will be particularly important for our algorithm.

Definition 5.1**.**

Let $x\in\mathbb{Z}^{N}$. For each $m$ define the function $h_{m}:[N]\rightarrow[m]$ by

[TABLE]

Moreover, define $\mathcal{P}_{m}(x)\in\mathbb{Z}^{m}$ to be such that

[TABLE]

We remind that $\omega$ is a $(2N)$-th root of unity rounded to fit in a word.

The following claim lies at the crux of our argument.

Lemma 5.2**.**

Given vectors $x,y,w\in\mathbb{Z}^{N}$ the vector $\mathcal{P}_{m}((x\star y)-w)$ up to $1/\mathrm{poly}(N)$ error can be computed in time

[TABLE]

Proof.

First, note that

[TABLE]

since $\mathcal{P}_{m}$ is a linear operator. We compute $\mathcal{P}_{m}(x),\mathcal{P}_{m}(y),\mathcal{P}_{m}(w)$ in time $O(m+\left(\|x\|_{0}+\|y\|_{0}+\|w\|_{0}\right)\log N)$ by computing $\omega^{j}$ for all $j\in\mathrm{supp}(x)\cup\mathrm{supp}(y)\cup\mathrm{supp}(w)$ using Taylor expansion of sine and cosine functions and keeping the first $\Theta(\log N)$ digits. We then compute, via FFT in time $O(m\log m)$, the vector $\mathcal{P}_{m}(x)\star\mathcal{P}_{m}(y)$. We claim that

[TABLE]

If we could write down the roots of unity using using an infinite number of bits we would have

[TABLE]

The only catch is that we have to round every root of unity to $\Theta(\log N)$ bits. We will prove the latter inequality and notice that the precision issue arises only in the usage of the identiy $\omega^{j}\cdot\omega^{j^{\prime}}=\omega^{j+j^{\prime}}$; when rounding to $\Theta(\log N)$ bits this introduces a negligible error of at most $1/\mathrm{poly}(N)$ which allows our argument to go through. Thus, we focus on proving the latter claim, via the following series of inequalities:

[TABLE]

where (1) to (2) follows by defition of convolution, (2) to (3) by definition of the $\mathcal{P}_{B}$ operator, (3) to (4) by expanding the product in (2), (5) to (6) by the trivial fact each element in $[n]$ is mapped to some element in $[m]$ via $h_{m}$, (6) to (7) by the fact that $(h_{m}(a)+h_{m}(b))\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m=(a\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m+b\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m)\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m=(a+b)\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m=h(a+b)$, (7) to (8) by introducing the auxilliary variabe $j^{\prime\prime}=j+j^{\prime}$, (8) to (9) by the fact that $\omega^{j^{\prime\prime}}$ can be pulled outside of the inner sum since in that scope $j^{\prime\prime}$ is fixed, and (9) to (10) since $h_{m}(j^{\prime\prime}\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m)=(j^{\prime\prime}\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m)\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m=(j^{\prime\prime}\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ m)=h_{m}(j^{\prime\prime})$ and the fact that the inner sum is the definition of convolution evaluated at point $j^{\prime\prime}$. ∎

In what follows $C$ is some sufficiently large absolute constant.

The following Lemma is important, since we are dealing with numbers with finite precision.

Lemma 5.3**.**

Let $a,b\in[2N]$ with $a\neq b$. If $\omega$ is rounded such that it fits in a word, then $|\omega^{a}-\omega^{b}|=\Omega(1/N)$.

Proof.

The quantity is minimized when $a=b+1$. For $N$ sufficiently large, it can then be approximated by an arc of length $2\pi/(2N)$, and since the word size $w$ is $\Omega(\log N)$ we get the desired result. ∎

The follows Lemma is a crucial building block of our algorithm.

Lemma 5.4**.**

Given $\omega^{j}$ for $j\in[2N]$ with $\omega$ rounded to fit a word, one can find $j$ in time $O(\log^{2}N)$.

Proof.

From the pair (real part of $\omega^{j}$,imaginary part of $\omega^{j}$) we can find in which of the four following sets $j$ lies in

[TABLE]

since each one corresponds to an arc of length approximately $\pi/4$ of the complex circle. After detecting the set (equivalently the corresponding arc of the complex circle) one can perform a standard ternary search to find $j$. Due to Lemma 5.3 $O(\log N)$ iterations suffice to find $j$. In order for to implement ternary search, at every time we need to compute $\omega^{\ell}$, for some $\ell\in[2N]$, which can be done in $\Theta(\log N)$ time up to error $1/\mathrm{poly}(N)$ by performing a Taylor expansion on $\omega^{\ell}$. This gives the desired result. We note that if we are allowed to precompute all $(2N)$-th roots of unity we may obtain running time $O(\log N)$. ∎

Before proceeding, we first note the following catch in Line 10. In fact, $v$ will be a complex number with integer magnitude plus (due to rounding errors) a small $1/\mathrm{poly}(n)$ error. Thus we can round $v$ to the closest integer to read off exactly the desired value.

The following Lemma resembles standard isolation-type arguments which appear in the sublinear-time sparse recovery literature. $C$ is a large enough absolute constant.

Lemma 5.5**.**

Let an integer $B$ such that $B>C\cdot\|(x\star y)-w\|_{0}$, and let $p$ be chosen at random from $[CB\log^{2}n]$. Then, with probability $1-q$, there exist at least $(1-\gamma)\|(x\star y)-w)\|_{0}$ indices $j\in\mathrm{supp}((x\star y)-w)$ such that

[TABLE]

where $2C^{-2}/q=\gamma$.

Proof.

Let $j,j^{\prime}\in\mathrm{supp}((x\star y)-w)$, with $j\neq j^{\prime}$. The hash function $h_{p}$ is not pairwise independent, but the following property, which suffices for our purpose, holds

[TABLE]

To see that, observe first that in order for $h_{p}(j)=h_{p}(j^{\prime})$ to hold, it must be the case that $p$ is a divisor of $j-j^{\prime}$. Since $j-j^{\prime}\leq N$ there there can be at most $\lceil\log N\rceil$ prime divisors of $j-j^{\prime}$, otherwise $j-j^{\prime}$ would be at least $2^{\left\lceil\log N\right\rceil+1}>N$. By the Prime Number Theorem, there exist at least $(C/2)B\log N$ primes in $[CB\log^{2}N]$, and hence a random prime will be one of the divisors of $j-j^{\prime}$ with probability $2/(CB)$. All the primes in $[CB\log^{2}N]$ can be found in $O(B\log^{2}N\cdot\log\log(B\log^{2}N))$ time using Eratosthene’s sieve, and then we can sample uniformly at random from them.

Le the random variable $X_{j}$ be the indicator variable of the event

[TABLE]

Its expected value is $\mathbb{E}\left[X_{j}\right]=\mathbb{P}\left[\mathcal{E}_{j}\leavevmode\nobreak\ \mathrm{holds}\right]\leq(\|(x\star y)-w)\|_{0}-1)\cdot(2/CB)\leq 2C^{-2}$, by a union-bound. We have that

[TABLE]

By Markov’s inequality, with probability $1-q$ there exist at most $\gamma\|(x\star y)-w\|_{0}$ indices $j\in\mathrm{supp}((x\star y)-w)$ such that $X_{j}=1$, if $2C^{-2}/q=\gamma$. This finishes the proof of the claim. ∎

We proceed by analyzing the iterative loop in Algorithm 2.

Lemma 5.6**.**

Let the constants $C,\gamma,q$ be as in Lemma 5.5 with $q\leq 2^{-12/5}$, and assume that $B>C\|(x\star y)-w\|_{0}$. If $(x\star y)-w$ is not the zero vector, then with probability $1-\delta$ the subroutine $\textsc{Locate}(x,y,w,B,\delta)$ returns a vector $z$ such that

[TABLE]

Proof.

Fix $t\in[5\log(1/\delta)]$, and assume that $q,C,\gamma$ satisfy $2C^{-2}/p=\gamma$ We have that

[TABLE]

The condition of Lemma 5.5 holds, so with probability $1-q$ its conclusion also holds. Condition on that event and consider the at least $(1-\gamma)\|(x\star y)-w\|_{0}$ indices in $\|(x\star y)-w\|_{0}$, for which the conclusion of Lemma 5.5 holds. Fix such an index $j^{*}$ and let $i^{*}=h_{p}(j^{*})$. Due to the isolation property, we have that

[TABLE]

Now, due to Lemma 5.3 subroutine $\textsc{Locate}(x,y,w,B,\delta)$ will infer $j^{*}$ correctly from $(\mathcal{P}_{\sigma,B}(x\star y))_{i^{*}}-(\mathcal{P}_{\sigma,B}w)_{i^{*}}$, as well as $(x\star y)_{j^{*}}-w_{j^{*}}$. We will say $j^{*}$ is recognised in repetition $t$.

For the rest of the proof, unfix $t$. Since the conclusion of Lemma 5.5 holds with probability $1-q$, the number of $t\in[5\log(1/\delta)]$ for which the conclusion of the Lemma holds is at least $4\log(1/\delta)$ with probability $1-\delta$ since

[TABLE]

as long as $q\leq 2^{-12/5}\cdot$.

Let us call for convenience the above pairs good. Thus, with probability $1-\delta$ the number of pairs $(j,t)$ for which $j$ is not recognised in repetition $t$ is at most

[TABLE]

Hence there exist at most $\beta=4\gamma\|(x\star y)-w\|_{0}$ indices which are recognized in less than $(3/4)\cdot 5\log(1/\delta)$ repetitions, otherwise the number of not good pairs $(j,t)$ is at least

[TABLE]

which does not hold for $\|(x\star y)-w\|_{0}>0$. Moreover, there can be at most $\gamma\|(x\star y)-w\|_{0}$ indices that do not belong in $\mathrm{supp}((x\star y)-w)$, and which were mistakenly inserted into $z$. This gives in total the factor of $5\gamma$. ∎

Lemma 5.7**.**

Let $\gamma<1/10$, and let also $B$ be an integer such that $B>C\|(x\star y)\|_{0}$. Then the routine $\textsc{HashAndIterate}(x,y,B,\delta)$ returns an $\|x\star y\|_{0}$-sparse vector $r$ such that $r=x\star y$, with probability $1-\delta$. Moreover, the running time is

[TABLE]

Proof.

It is an easy induction to show that at each step $\|(x\star y)-w^{(r)}\|_{0}\leq(5\gamma)^{r}\|x\star y\|_{0}$, with probability $1-\delta r/\log B$, so the total failure probability is $\delta$. Conditioned on the previous events happending, we have $x\star y-w^{(\lceil\log B\rceil)}$ is the all-zeros vector since $\|x\star y\|_{0}\leq(4\gamma)^{\lceil\log B\rceil}\|(x\star y)-w\|_{0}<1$. This gives that $w^{(\lceil\log B\rceil)}=x\star y$.

The running time for $\textsc{Locate}(x,y,w,B_{r},\delta_{r})$, since $\|w\|_{0}\leq 2B$ at all times is (ignoring constant factors for ease of exposition)

[TABLE]

due to Lemma 5.2 and Lemma 5.4.

So the total running time of $\textsc{HashAndIterate}(x,y,B,\delta)$ becomes, by summing over all $\lceil\log B\rceil$ rounds (ignoring constant factors for ease of exposition)

[TABLE]

∎

The following Lemma is a standard claim which follows by the fact that a degree $n$ polynomial over the prime field $\mathbb{Z}_{p}$ has at most $n$ roots. We give a sketch of the proof.

Lemma 5.8**.**

There exists a procedure $\textsc{EqualityTesting}(x,y,w)$, which runs in time $O(\|x\|_{0}+\|y\|_{0}+\|w\|_{0})\log N\log(1/\delta)+\widetilde{O}(\log^{2}N\cdot\log(N/\delta)\cdot\log(1/\delta))$, and answers whether $x\star y=w$ with probability $1-\delta$.

Proof.

Let $c^{\prime}$ large enough. We pick a random prime $p\in[c^{\prime}N,2c^{\prime}N]$, by picking a random number in that interval and running the Miler-Rabin primality test with target failure probability $\delta/3$. We form polynomials $f_{x},f_{y},f_{w}$ that have $x,y,w$ as their coefficients respectively. We then pick $\Theta(\log(3/\delta))$ random elements in $\mathbb{Z}_{p}$ and check whether $(f_{x}(r)\cdot f_{y}(r))\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ p=f_{w}(r)\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ p$ or not. We return Yes if this is the case for all chosen, and No otherwise. To evalute each each of the polynomials we need time (number of coefficients) $\cdot\log N$, in order to perform repetated squaring. The Miller-Rabin test takes time $\widetilde{O}(\log^{2}N\cdot\log(1/\delta))$, and sampling a prime with probability $1-\delta/3$ needs time $\Theta(\log(N/\delta))$.

∎

We are now ready to prove our main theorem.

Proof.

Let $c$ be a sufficiently small constant and $C$ a sufficiently large constant. For $r=1,2,\ldots$, one by one we set $B_{r}\leftarrow C\cdot 2^{r}$ and $\delta_{r}=c\cdot r^{-2}$, run $\textsc{HashAndIterate}(x,y,B_{r},\delta_{r})$ to obtain $z$, and feed it to $\textsc{EqualityTesting}(x,y,z,cr^{-2})$. We stop when the latter procedure returns Yes. The total failure probability thus is at most

[TABLE]

Conditioned on the aforementioned event happening, the total running time is (ignoring constants)

[TABLE]

The running time in the first two lines follows by invoking Lemma 5.7 and a straightforward summation over all $O(\log(\|x\star y\|_{0}))$ rounds, and the third line is the cost of invoking Lemma 5.8 in every round.

The final expression now follows by recalling that $N=2n$. ∎