# Error correction in fast matrix multiplication and inverse

Daniel S. Roche
Computer Science Department
United States Naval Academy
Annapolis, Maryland, U.S.A

###### Abstract

We present new algorithms to detect and correct errors in the product of two
matrices, or the inverse of a matrix, over an arbitrary field. Our algorithms
do not require any additional information or encoding other than the original
inputs and the erroneous output. Their running time is softly linear in the
number of nonzero entries in these matrices when the number of errors is
sufficiently small, and they also incorporate fast matrix multiplication so
that the cost scales well when the number of errors is large. These algorithms
build on the recent result of Gasieniec et al. [2017] on correcting matrix
products, as well as existing work on verification algorithms, sparse low-rank
linear algebra, and sparse polynomial interpolation.

##  1 Introduction

Efficiently and gracefully handling computational errors is critical in modern
computing. Such errors can result from short-lived hardware failures,
communication noise, buggy software, or even malicious third-party servers.

The first goal for fault-tolerant computing is _verification_. Freivalds
[1979] presented a linear-time algorithm to verify the correctness of a single
matrix product. Recently, efficient verification algorithms for a wide range
of computational linear algebra problems have been developed [Kaltofen et al.,
2011, Dumas and Kaltofen, 2014, Dumas et al., 2016, 2017].

A higher level of fault tolerance is achieved through _error correction_ ,
which is the goal of the present study. This is a strictly stronger model than
verification, where we seek not only to identify when a result is incorrect,
but also to compute the correct result if it is “close” to the given,
incorrect one. Some recent error-correction problems considered in the
computer algebra literature include Chinese remaindering [Goldreich et al.,
1999, Khonji et al., 2010, Böhm et al., 2015], system solving [Boyer and
Kaltofen, 2014, Kaltofen et al., 2017], and function recovery [Comer et al.,
2012, Kaltofen and Yang, 2013].

A generic approach to correcting errors _in transmission_ would be to use an
error-correcting code, writing the result with some redundancy to support
transmission over a noisy channel. But this approach requires extra bandwidth,
is limited to a fixed number of errors, and does not account for malicious
alterations or mistakes in the computation itself.

Instead, we will use the information in the problem itself to correct errors.
This is always possible by re-computing the result, and so the goal is always
to correct a small number of errors in (much) less time than it would take to
do the entire computation again. Our notion of error correction makes it
possible to correct an unbounded number of mistakes that occur either in
computation or transmission, either at random or by a malicious party.

Some other applications of error correction have nothing to do with
alterations or mistakes. One example is sparse matrix multiplication, which
can be solved with an error-correction algorithm by simply assuming the
“erroneous” matrix is zero. Another example is computing a result over
$\mathbb{Z}$ using Chinese remaindering modulo small primes, where the entries
have varying bit lengths. The small entries are determined modulo the first
few primes, and the remaining large entries are found more quickly by applying
an error correction algorithm using (known) small entries.

###  1.1 Our results

#### Matrix multiplication with errors

We consider two computational error correction problems. The first problem is
the same as in Gasieniec et al. [2017], correcting a matrix product. Given
$A,B,C\in\mathsf{F}^{n\times n}$, the task is to compute the unique error
matrix $E\in\mathsf{F}^{n\times n}$ such that $AB=C-E$.

Algorithm 5 solves this problem using

$\widetilde{O}\mathopen{}\left(t+k\cdot\min\left(\left\lceil\tfrac{t}{r}\right\rceil,n/\min(r,\tfrac{k}{r})^{3-\omega}\right)\right)\mathclose{}$

field operations, where

  * •

$t=\\#A+\\#B+\\#C$ is the number of nonzero entries in the input matrices, at
most $O(n^{2})$;

  * •

$k=\\#E$ is the number of errors in the given product $C$; and

  * •

$r\leq n$ is the number of distinct rows (or columns) in which errors occur,
whichever is larger.

| Spread-out errors | Compact errors
---|---|---
| $r=\min(n,k)$ | $r=\sqrt{k}$
| $k\leq n$ | $k>n$ |
Sparse | $t+k+n$ | $k\left\lceil\tfrac{t}{n}\right\rceil+n$ | $\sqrt{k}t+k+n$
Dense | $t+kn$ | $k^{\omega-2}n^{4-\omega}$ | $t+k^{(\omega-1)/2}n$
$k^{0.38}n^{1.63}$ | $t+k^{0.69}n$
Table 1: Soft-oh cost of Algorithm 5 in different cases. The two rows indicate
whether sparse or dense rectangular multiplication is used as a subroutine.

An even more detailed complexity statement, incorporating rectangular matrices
and a failure probability, can be found in Theorem 7.1. To understand the
implications of this complexity, consider six cases, as summarized in Table 1.

The cost of our algorithm depends on the locations of the errors. If there are
a constant number of errors per row or column, we can use sparse
multiplication to achieve softly-linear complexity
$\widetilde{O}\mathopen{}\left(t+n+k\right)\mathclose{}$. This is similar to
[Gasieniec et al., 2017] except that we do not require the input matrices to
be dense.

If $k>n$ errors are spread out among all rows and columns, the worst-case cost
using dense multiplication is
$\widetilde{O}\mathopen{}\left(k^{0.38}n^{1.63}\right)\mathclose{}$. As the
number of errors grows to $n^{2}$, this cost approaches that of normal matrix
multiplication without error correction.

In the other extreme, if the errors are isolated to a square (but not
necessarily contiguous) submatrix, then the worst-case cost using dense
multiplication is
$\widetilde{O}\mathopen{}\left(t=k^{0.69}n\right)\mathclose{}$. Again this
scales to $n^{\omega}$ as $k\to n^{2}$, but the cost is better for $k<n^{2}$.
This situation may also make sense in the context of distributed computing,
where each node in a computing cluster might be assigned to compute one
submatrix of the product.

Our algortihm is never (asymptotically) slower than that of [Gasieniec et al.,
2017] or the $O(n^{\omega})$ cost of naïve recomputation, but it is faster in
many cases: for example when the inputs are sparse, when there are many
errors, and when the errors are compactly located.

#### Matrix inverse with errors

The second problem we consider is that of correcting computational errors in a
matrix inverse. Formally, given $A,B\in\mathsf{F}^{n\times n}$, the task is to
compute the unique error matrix $E\in\mathsf{F}^{n\times n}$ such that
$A^{-1}=B+E$.

Algorithm 6 shows how to solve this problem using

$\widetilde{O}\mathopen{}\left(t+nk/\min(r,\tfrac{k}{r})^{3-\omega}+r^{\omega}\right)\mathclose{}$

field operations, where the parameters $t,k,n$ are the same as before, and $r$
is the number of rows or columns which contain errors, whichever is _smaller_.
This is the same complexity as our algorithm for matrix multiplication with
errors, plus an additional $r^{\omega}$ term.

The additional $r^{\omega}$ term makes the situation here less fractured than
with the multiplication algorithm. In the case of spread-out errors, the cost
is $\widetilde{O}\mathopen{}\left(t+nk+k^{\omega}\right)\mathclose{}$, which
becomes simply
$\widetilde{O}\mathopen{}\left(n^{2}+k^{\omega}\right)\mathclose{}$ when the
input matrices are dense. The case of compact errors is better, and the
complexity is just
$\widetilde{O}\mathopen{}\left(t+nk^{(\omega-1)/2}\right)\mathclose{}$, the
same as that of our multiplication algorithm

As before, the cost of this algorithm approaches the $O(n^{\omega})$ cost of
the naïve inverse computation as the number of errors $k$ grows to $n^{2}$.

#### Domains

The algorithms we present work over any field and require only that a field
element of order larger than $n$ can be found. Because our algorithms
explicitly take $O(n)$ time to compute the powers of this element anyway, the
usual difficulties with efficiently generating high-order elements in finite
fields do not arise.

Over the integers $\mathbb{Z}$ or rationals $\mathbb{Q}$, the most efficient
approach would be to work modulo a series of primes and use Chinese
remaindering to recover the result. A similar approach could be used over
polynomial rings.

Over the reals $\mathbb{R}$ or complex numbers $\mathbb{C}$, our algorithms
work only under the unrealistic assumption of infinite precision. We have not
considered the interesting question of how to recover from the two types of
errors (noise and outliers) that may occur in finite precision.

###  1.2 Related work

Let $\omega$ be a constant between 2 and 3 such that two $n\times n$ matrices
can be multiplied using $O(n^{\omega})$ field operations. In practice,
$\omega=3$ for dimensions $n$ up to a few hundred, after which Strassen-
Winograd multiplication is used giving $\omega=<2.81$. The best asymptotic
algorithm currently gives $\omega<2.3728639$ [Le Gall, 2014].

Even though the practical value of $\omega$ is larger than the asymptotic one,
matrix multiplication routines have been the subject of intense implementation
development for decades, and highly-tuned software is readily available for a
variety of architectures [Dumas et al., 2008, Wang et al., 2013]. Running
times based on $n^{\omega}$ have practical as well as theoretical
significance; the $\omega$ indicates where an algorithm is able to take
advantage of fast low-level matrix multiplication routines.

Special routines for multiplying sparse matrices have also been developed.
Writing $t$ for the number of nonzero entries in the input, the standard row-
or column-wise sparse matrix multiplication costs $O(tn)$ field operations.
Yuster and Zwick [2005] improved this to
$\widetilde{O}\mathopen{}\left(n^{2}+t^{0.697}n^{1.197}\right)\mathclose{}$.
Later work also incorporates the number $k$ of nonzeros in the product matrix.
Lingas [2009] gave an output-sensitive algorithm with running time
$O(k^{0.186}n^{2})$. Amossen and Pagh [2009] have another with complexity
$\widetilde{O}\mathopen{}\left(t^{0.667}k^{0.667}+t^{0.862}k^{0.408}\right)\mathclose{}$,
an improvement when the outputs are not too dense. Pagh [2013] developed a
different approach with complexity
$\widetilde{O}\mathopen{}\left(t+kn\right)\mathclose{}$; that paper also
includes a nice summary of the state of the art. Our multiplication with
errors algorithm can also be used for sparse multiplication, and it provides a
small improvement when the input and output are both not _too_ sparse.

The main predecessor to our work is the recent algorithm of Gasieniec et al.
[2017], which can correct up to $k$ errors in the product of two $n\times n$
matrices using $\widetilde{O}\mathopen{}\left(n^{2}+kn\right)\mathclose{}$
field operations. Their approach is based in part on the earlier work of Pagh
[2013] and makes clever use of hashing and fast Fourier transforms in order to
achieve the stated complexity.

Computing the inverse of a matrix has the same asymptotic complexity
$O(n^{\omega})$ as fast matrix multiplication [Bunch and Hopcroft, 1974]. In
practice, using recent versions of the FFPACK library [Dumas et al., 2008] on
an Intel desktop computer, we find that a single $n\times n$ matrix inversion
is only about 33% more expensive than a matrix-matrix product of the same
size.

##  2 Overview

#### Multiplication with errors

A high-level summary of our multiplication algorithm is as follows:

  1. 1.

Determine which rows in the product contain errors.

  2. 2.

Write the unknown sparse matrix of errors as $E=C-AB$. Remove the rows with no
errors from $E,C,$ and $A$, so we have $E^{\prime}=C^{\prime}-A^{\prime}B$.

  3. 3.

Treating the rows of $E^{\prime}$ as sparse polynomials, use structured linear
algebra and fast matrix multiplication to evaluate each row polynomial at a
small number of points.

  4. 4.

Use sparse polynomial interpolation to recover at least half of the rows from
their evaluations.

  5. 5.

Update $C$ and iterate $O(\log k)$ times to recover all $k$ errors.

To determine the rows of $E$ which are nonzero in Step 1, we apply a simple
variant of Frievalds’ verification algorithm: for a random vector
$\mathbf{v}$, the nonzero rows of $E\mathbf{v}$ are probably the same as the
nonzero rows of $E$ itself.

The evaluation/interpolation approach that follows is reminiscent of that in
Gasieniec et al. [2017], with two important differences. First, using sparse
polynomial techniques rather than hashing and FFTs means the complexity scales
linearly with the number of nonzeros $t$ in the input rather than the dense
matrix size $n^{2}$. Second, by treating the rows separately rather than
recovering all errors at once, we are able to incorporate fast matrix
multiplication so that our worst-case cost when all entries are erroneous is
never more than $O(n^{\omega})$, the same as that of naïve recomputation.

Section 4 provides details on the Monte Carlo algorithm to determine the rows
and columns where errors occur (Step 1 above), Section 5 presents a variant of
sparse polynomial interpolation that suits our needs for Step 4, and Section 6
explains how to perform the row polynomial evaluations in Step 3. A certain
rectangular matrix multiplication in this step turns out to dominate the
overall complexity. We connect all these pieces of the multiplication
algorithm in Section 7.

#### Inverse with errors

Our algorithm for correcting errors in a matrix inverse is based on the
multiplication algorithm, but with two complications. Writing $A$ for the
input matrix, $B$ for the erroneous inverse, and $E$ for the (sparse) matrix
of errors in $B$, we have:

$AE=I-AB,\qquad EA=I-BA.$

We want to compute $E\mathbf{v}$ for a random vector $\mathbf{v}$ in order to
determine which rows of $E$ are nonzero. This is no longer possible directly,
but using the second formulation as above, we can compute $EA\mathbf{v}$.
Because $A$ must be nonsingular, $A\mathbf{v}$ has the same distribution as a
random vector, so this approach still works.

The second complication is that removing zero rows of $E$ does not change the
dimension of the right-hand side, but only removes corresponding _columns_
from $A$ on the left-hand side. We make use of the recent result of Cheung et
al. [2013], which shows how to quickly find a maximal set of linearly
independent rows in a sparse matrix. This allows us to find a small submatrix
$X$ of $A$ such that $XE^{\prime}=I^{\prime}-A^{\prime}B$, with
$E^{\prime},I^{\prime},A^{\prime}$ being compressed matrices as before.
Because the size of $X$ depends on the number of errors, it can be inverted
much more quickly than re-computing the entire inverse of $A$.

The dominating cost in most cases is the same rectangular matrix product as
needed as in matrix multiplication with error correction. However, here we
also have the extra cost of finding $X$ and computing its inverse, which
dominates the complexity when there are a moderate number of errors and their
locations are spread out.

This algorithm — which depends on the subroutines of Sections 4, 5 and 6 — is
presented in detail in Section 8. .

##  3 Notation and preliminaries

The soft-oh notation $\widetilde{O}\mathopen{}\left(\cdots\right)\mathclose{}$
is the same as the usual big-oh notation but ignoring sub-logarithmic factors:
$f\in\widetilde{O}\mathopen{}\left(g\right)\mathclose{}$ if and only if $f\in
O\mathopen{}\left(g\log^{O(1)}g\right)\mathclose{}$, for some runtime
functions $f$ and $g$.

We write $\mathbb{N}$ for the set of nonnegative integers and $\mathbb{F}_{q}$
for the finite field with $q$ elements.

For a finite set $\mathcal{S}$, the number of elements in $\mathcal{S}$ is
denoted $\\#\mathcal{S}$, and $\mathbb{P}(\mathcal{S})$ is the powerset of
$\mathcal{S}$, i.e., the set of subsets of $\mathcal{S}$.

The number of nonzero entries in a matrix $M\in\mathsf{F}^{m\times n}$ is
written as $\\#M$. Note that $\\#M\leq mn$, and if $M$ has rank $k$ then
$\\#M\geq k$.

For a univariate polynomial $f\in\mathsf{F}[x]$,
$\mathsf{supp}(f)\subseteq\mathbb{N}$ denotes the exponents of nonzero terms
of $f$. Note that $\\#\mathsf{supp}(f)\leq\deg f+1$.

We assume that the number of field operations to multiply two (dense)
polynomials in $\mathsf{F}[x]$ with degrees less than $n$ is
$\widetilde{O}\mathopen{}\left(n\right)\mathclose{}$. This is justified by the
generic algorithm from Cantor and Kaltofen [1991] which gives
$O\mathopen{}\left(n\log n\operatorname{loglog}n\right)\mathclose{}$, or more
result results of Harvey et al. [2017] which improve this to
$O\mathopen{}\left(n\log n\,8^{\log^{*}n}\right)\mathclose{}$ for finite
fields.

As discussed earlier, we take $\omega$ with $2\leq\omega\leq 3$ to be any
feasible exponent of matrix multiplication. By applying fast matrix
multiplication in blocks, it is also possible to improve the speed of
rectangular matrix multiplication in a straightforward way. In particular:

######  Fact 3.1.

The product of any $m\times\ell$ matrix times a $\ell\times n$ matrix can be
found using $O\mathopen{}\left(m\ell
n/\min(m,\ell,n)^{3-\omega}\right)\mathclose{}$ ring operations.

Note that it might be possible to improve this for certain values of
$m,\ell,n$ using Le Gall [2012], but that complexity is much harder to state
and we have not investigated such an approach.

We restate two elementary facts on sparse matrix multiplication.

######  Fact 3.2.

For any $M\in\mathsf{F}^{m\times n}$ and $\mathbf{v}\in\mathsf{F}^{n}$, the
matrix-vector product $M\mathbf{v}$ can be computed using $O(\\#M)$ field
operations.

######  Corollary 3.3.

For any $A\in\mathsf{F}^{m\times\ell}$ and $B\in\mathsf{F}^{\ell\times n}$,
their product $AB$ can be computed using $O(\\#A\cdot n)$ field operations.

##  4 Identifying nonzero rows and columns

Our algorithms for error correction in matrix product and inverse computation
both begin by determining which rows and columns contain errors. This is
accomplished in turn by applying a random row or column vector to the unknown
error matrix $E$.

For now, we treat the error matrix $E$ as a _black box_ , i.e., an oracle
matrix which we can multiply by a vector on the right or left-hand side. The
details of the black box construction differ for the matrix multiple and
inverse problem, so we delay those until later.

The idea is to apply a row or column vector of random values to the unknown
matrix $E$, and then examine which entires in the resulting vector are zero.
This is exactly the same idea as the classic Monte Carlo verification
algorithm of Freivalds [1979], which we extend to which rows or columns in a
matrix are nonzero. This black-box procedure is detailed in Algorithm 1
FindNonzeroRows.

Input: Black box for right-multiplication by an unknown matrix
$M\in\mathsf{F}^{m\times n}$ and error bound
$\epsilon\in\mathbb{R},0<\epsilon<1$

Output: Set $\mathcal{J}\subseteq\\{0,\ldots,m-1\\}$ that with probability at
least $1-\epsilon$ contains the index of all nonzero rows in $M$

1
$\ell\leftarrow\left\lceil\log_{\\#\mathsf{F}}\tfrac{m}{\epsilon}\right\rceil$

2 $V\leftarrow$ matrix in $\mathsf{F}^{n\times\ell}$ with uniform random
entries from $\mathsf{F}$

3 $B\leftarrow MV$ using the black box

return _Indices of rows of $B$ which are not all zeros_

Algorithm 1 FindNonzeroRows$(V\mapsto MV,\epsilon)$

The running time of FindNonzeroRows depends entirely on the size of
$\mathsf{F}$ and the cost of performing one black-box apply over the field
$\mathsf{F}$. The correctness depends on the parameter $\epsilon$.

######  Lemma 4.1.

FindNonzeroRows always returns a list of nonzero row indices in $M$ With
probability at least $1-\epsilon$, it returns the index of _every_ nonzero row
in $M$.

###### Proof.

Consider a single row $\mathbf{u}^{T}$ of $M$. If $\mathbf{u}$ is zero, then
the corresponding row of $MV$ will always be zero; hence every index returned
must be a nonzero row.

Next assume $\mathbf{u}^{T}$ has at least one nonzero entry, say at index $i$,
and consider one row $\mathbf{v}\in\mathsf{F}^{n}$ of the random matrix $V$.
Whatever the other entries in $\mathbf{u}$ and $\mathbf{v}$ are, there is
exactly one value for the $i$’th entry of $v$ such that
$\mathbf{u}^{T}\mathbf{v}=0$. So the probability that
$\mathbf{u}^{T}\mathbf{v}=0$ is $\tfrac{1}{\\#\mathsf{F}}$, and by
independence the probability that $\mathbf{u}^{T}V=\mathbf{0}^{1\times\ell}$
is $1/(\\#\mathsf{F})^{\ell}$.

Because $M$ has at most $m$ nonzero rows, the union bound tells us that the
probability _any one_ of the corresponding rows of $MV$ is zero is at most
$m/(\\#F)^{\ell}$. From the definition of $\ell$, this ensures a failure
probability of at most $\epsilon$. ∎

Because we mostly have in mind the case that $\mathsf{F}$ is a finite field,
we assume in FindNonzeroRows and in the preceding proof that it is possible to
sample uniformly from the entire field $\mathsf{F}$. However, the same results
would hold if sampling from any fixed subset of $\mathsf{F}$ (perhaps
increasing the size of $\ell$ if the subset is small).

##  5 Batched low-degree sparse interpolation

Our algorithms use techniques from sparse polynomial interpolation to find the
locations and values of erroneous entries in a matrix product or inverse.

Consider an unknown univariate polynomial $f\in\mathsf{F}[x]$ with degree less
than $n$ and at most $s$ nonzero terms. The nonzero terms of $f$ can be
uniquely recovered from evaluations $f(1),f(\theta),\ldots,f(\theta^{2s-1})$
at $2s$ consecutive powers [Ben-Or and Tiwari, 1988, Kaltofen and Yagati,
1989].

Doing this efficiently in our context requires a few variations to the typical
algorithmic approach. While the state of the art seen in recent papers such as
[Javadi and Monagan, 2010, Kaltofen, 2010, Arnold et al., 2015, Huang and Gao,
2017] achieves softly-linear complexity in terms of the sparsity $t$, there is
either a higher dependency on the degree, a restriction to certain fields, or
both.

Here, we show how to take advantage of two unique aspects of the problem in
our context. First, we will interpolate a _batch_ of sparse polynomials at
once, on the same evaluation points, which allows for some useful
precomputation. Secondly, the maximum degree of any polynomial we
interpolation is $n$, the column dimension of the unknown matrix, and we are
allowed linear-time in $n$.

The closest situation in the literature is that of van der Hoeven and Lecerf
[2013], who show how to recover a sparse polynomial if the exponents are
known. In our case, we can tolerate a much larger set containing the possible
exponents by effectively amortizing the size of this exponent set over the
group of batched sparse interpolations.

Recall the general outline of Prony’s method for sparse interpolation of
$f\in\mathsf{F}[x]$ from evaluations $(f(\theta^{i}))_{0\leq i<2s}$:

  1. 1.

Find the minimum polynomial $\Gamma\in\mathsf{F}[z]$ of the sequence of
evaluations.

  2. 2.

Find the roots $v_{1},\ldots,v_{s}$ of $\Gamma$.

  3. 3.

Compute the discrete logarithm base $\theta$ of each $v_{i}$ to determine the
exponents of $f$.

  4. 4.

Solve a transposed Vandermonde system built from $v_{i}$ and the evaluations
to recover the coefficients of $f$.

Steps 1 and 4 can be performed over any field using fast multipoint evaluation
and interpolation in $\widetilde{O}\mathopen{}\left(s\right)\mathclose{}$
field operations [Kaltofen and Yagati, 1989]. However, the other steps depend
on the field and in particular Step 3 can be problematic over large finite
fields.

To avoid this issue, we first pre-compute all possible roots $\theta^{e_{i}}$
for any minpoly $\Gamma_{i}$ in any of the batch of $r$ sparse interpolations
that will be performed. Although the set of possible roots is larger than the
degree of any single minpoly $\Gamma_{i}$, in our case the set is always
bounded by the dimension of the matrix, so performing this precomputation once
is worthwhile.

###  5.1 Batched root finding

While fast root-finding procedures over many fields have already been
developed, we present a simple but efficient root-finding procedure for our
specific case that has the advantage of running in softly-linear time over any
coefficient field. The general idea is to first compute a coprime basis (also
called a gcd-free basis) for the minpolys $\Gamma_{i}$, then perform multi-
point evaluation with the precomputed list of possible roots.

The details of procedure FindRoots are in Algorithm 2 below. In the algorithm,
we use a _product tree_ constructed from a list of $t$ polynomials. This is a
binary tree of height $O\mathopen{}\left(\log t\right)\mathclose{}$, with the
original polynomials at the leaves, and where every internal node is the
product of its two children. In particular, the root of the tree is the
product of all $t$ polynomials. The total size of a product tree, and the cost
of computing it, is softly-linear in the total size of the input polynomials;
see Borodin and Munro [1975] for more details.

Input: Set of possible roots $\mathcal{Q}\subseteq\mathsf{F}$ and $r$
polynomials $\Gamma_{1},\ldots,\Gamma_{r}\in\mathsf{F}[z]$

Output: A list of sets
$[\mathcal{R}_{1},\ldots,\mathcal{R}_{r}]\in\mathbb{P}(\mathcal{Q})^{r}$ such
that $g_{i}(\alpha_{j})=0$ for every $1\leq i\leq r$ and
$\alpha_{j}\in\mathcal{R}_{i}$

1

2$\\{g_{0},\ldots,g_{t-1}\\}\leftarrow$ coprime basis of
$\\{\Gamma_{i}\\}_{1\leq i\leq r}$ using Algorithm 18.1 of Bernstein [2005]

3

/* Compute product tree from bottom-up */

4 for  _$j\leftarrow 0,1,\ldots,t-1$_ do $M_{0,j}\leftarrow g_{j}$

5 for  _$i\leftarrow 1,2,\ldots,\left\lceil\log_{2}t\right\rceil$_ do

6 for  _$j\leftarrow 0,1,\ldots,\left\lceil t/2^{i}\right\rceil-1$_ do

7 $M_{i,j}\leftarrow M_{i-1,2j}\cdot M_{i-1,2j+1}$

8

/* Find roots of basis elements */

9 $S_{\left\lceil\log_{2}t\right\rceil+1,0}\leftarrow\mathcal{Q}$

10 for
_$i\leftarrow\left\lceil\log_{2}t\right\rceil,\left\lceil\log_{2}t\right\rceil-1,\ldots,0$_
do

11 for  _$j\leftarrow 0,1,\ldots,\left\lceil t/2^{i}\right\rceil-1$_ do

12 Evaluate $M_{i,j}$ at each point in $S_{i+1,\left\lfloor j/2\right\rfloor}$
using fast multipoint evaluation

13 $S_{i,j}\leftarrow$ roots of $M_{i,j}$ found on previous step

14

/* Determine roots of original polynomials */

15 Compute exponents $(e_{i,j})_{1\leq i\leq r,0\leq j<t}$ such that
$\Gamma_{i}=g_{0}^{e_{i,0}}\cdots g_{t-1}^{e_{i,t-1}}$ for all $i$, using
algorithm 21.2 of Bernstein [2005]

16

17for  _$i\leftarrow 1,2,\ldots,r$_ do
$\mathcal{R}_{i}\leftarrow\bigcup_{0\leq j<t,e_{i,j}\geq 1}S_{0,j}$

Algorithm 2 FindRoots$(\mathcal{Q},[\Gamma_{1},\ldots,\Gamma_{r}])$

######  Lemma 5.1.

Given a set $\mathcal{Q}\subseteq\mathsf{F}$ with of size $\\#\mathcal{Q}=c$
and a list of $r$ polynomials $\Gamma_{1},\ldots,\Gamma_{r}\in\mathsf{F}[z]$
each with degree at most $s$, Algorithm 2 FindRoots correctly determines all
roots of all $\Gamma_{i}$’s which are in $\mathcal{Q}$ using
$\widetilde{O}\mathopen{}\left(rs+c\right)\mathclose{}$ field operations

###### Proof.

Algorithms 18.1 and 21.2 of Bernstein [2005] are deterministic and compute
(respectively) a coprime basis and then a factorization of the $\Gamma_{i}$’s
according to the basis elements.

The polynomials $M_{i,j}$ form a standard product tree from the polynomials in
the coprime basis. For any element of the product tree, its roots must be a
subset of the roots of its parent node. So the multi-point evaluations on
Algorithm 2 determine, for every polynomial in the tree, which elements of
$\mathcal{Q}$ are roots of that polynomial.

The cost of computing the coprime basis using Algorithm 18.1 in [Bernstein,
2005] is $\widetilde{O}\mathopen{}\left(rs\right)\mathclose{}$ with at least a
$\log^{5}(rs)$ factor. This gives the first term in the complexity (and
explains our use of soft-oh notation throughout). This cost dominates the cost
of constructing the product tree and factoring the $\Gamma_{i}$’s on Algorithm
2.

Using the standard product tree algorithm [Borodin and Munro, 1975], the cost
of multi-point evaluation of a degree-$n$ polynomial at each of $n$ points is
$\widetilde{O}\mathopen{}\left(n\right)\mathclose{}$ field operations. In our
algorithm, this happens at each level down the product tree. At each of
$\left\lceil\log_{2}t\right\rceil\in O(\log(rs))$ levels there are at most $c$
points (since the polynomials at each level are relatively prime). The total
degree at each level is at most $\deg\prod_{i}\Gamma_{i}\leq rs$. If $rs>c$,
then this cost is already bounded by that of determining the coprime basis.
Otherwise, the multi-point evaluations occur in $O(c/(rs))$ batches of $rs$
points each, giving the second term in the complexity statement. ∎

###  5.2 Batched sparse interpolation algorithm

We are now ready to present the full algorithm for performing simultaneous
sparse interpolation on $r$ polynomials based on $2s$ evaluations each.

Input: Bounds $r,n,s\in\mathbb{N}$, set
$\mathcal{S}\subseteq\\{0,\ldots,n-1\\}$ of possible exponents, high-order
element $\theta\in\mathsf{F}$, and evaluations $Y_{i,j}\in\mathsf{F}$ such
that $f_{i}(\theta^{j})=Y_{i,j}$ for all $1\leq i\leq r$ and $0\leq j<2s$

Output: Nonzero coefficients and corresponding exponents of
$f_{i}\in\mathsf{F}[x]$ or $0$ indicating failure, for each $1\leq i\leq r$

1 for  _$i=1,2,\ldots,r$_ do

2 $\Gamma_{i}\leftarrow$MinPoly$(Y_{i,0},\ldots,Y_{i,2s-1})$ using fast
Berlekamp-Massey algorithm

3

4$\mathcal{R}_{1},\ldots,\mathcal{R}_{r}\leftarrow\hyperlink{rootfind}{\texttt{FindRoots}{}}(\\{\,\theta^{d}\mid
d\in\mathcal{S}\,\\},\Gamma_{1},\ldots,\Gamma_{r})$

5 for  _$i=1,2,\ldots,r$_ do

6 if  _$\\#\mathcal{R}_{i}\neq\deg\Gamma_{i}$_ then $f_{i}\leftarrow 0$

7 else

8 $t_{i}\leftarrow\\#\mathcal{R}_{i}$

9
$\\{e_{i,1},\ldots,e_{i,t_{i}}\\}\leftarrow\\{\,d\in\mathcal{S}\mid\theta^{d}\in\mathcal{R}_{i}\,\\}$

10 $[a_{i,1},\ldots,a_{i,t_{i}}]\leftarrow$ solution of transposed Vandermonde
system from roots $[\theta^{e_{i,1}},\ldots,\theta^{e_{i,t_{t}}}]$ and right-
hand vector $[Y_{i,0},\ldots,Y_{i,t_{i}-1}]$

11 $f_{i}\leftarrow a_{i,1}x^{e_{i,1}}+\cdots+a_{i,t_{i}}x^{e_{i,t_{i}}}$

12

Algorithm 3 MultiSparseInterp$(r,n,s,\mathcal{S},\theta,Y)$

######  Theorem 5.2.

Let $\mathsf{F}$ be a field with an element $\theta\in\mathsf{F}$ whose
multiplicative order is at least $n$, and suppose
$f_{1},\ldots,f_{r}\in\mathsf{F}[x]$ are unknown univariate polynomials. Given
bounds $n,s\in\mathbb{N}$, a set $\mathcal{S}\subseteq\mathbb{N}$ of possible
exponents, and evaluations $f_{i}(\theta^{j})$ for all $1\leq i\leq r$ and
$0\leq j<2s$, Algorithm 3 MultiSparseInterp requires
$\widetilde{O}\mathopen{}\left(rs+\\#\mathcal{S}\log n\right)\mathclose{}$
field operations in the worst case.

The algorithm correctly recovers every polynomial $f_{i}$ such that
$\\#\mathsf{supp}f_{i}\leq s$, $\mathsf{supp}f_{i}\subseteq\mathcal{S}$, and
$\deg f_{i}<n$.

###### Proof.

From Lemma 5.1, the cost of FindRoots is
$\widetilde{O}\mathopen{}\left(rs+\\#\mathcal{S}\right)\mathclose{}$. This
dominates the cost of the minpoly and transposed Vandermonde computations,
which are both $\widetilde{O}\mathopen{}\left(rs\right)\mathclose{}$ [Kaltofen
and Yagati, 1989]. Using binary powering, computing the set
$\\{\,\theta^{d}\mid d\in\mathcal{S}\,\\}$ on Algorithm 3 requires
$O(\\#\mathcal{S}\log n)$ field operations.

If $f_{i}$ has at most $t_{i}\leq s$ nonzero terms, then the minpoly
$\Gamma_{i}$ computed on Algorithm 3 is a degree-$t_{i}$ polynomial with
exactly $t_{i}$ distinct roots. If $\mathsf{supp}(f_{i})\subseteq\mathcal{S}$,
then by Lemma 5.1, FindRoots correctly finds all these roots. Because these
are the only terms in $f_{i}$, solving the transposed Vandermonde system
correctly determines all corresponding coefficients; see [Kaltofen and Yagati,
1989] or [van der Hoeven and Lecerf, 2013] for details. ∎

##  6 Performing evaluations

Consider the matrix multiplication with errors problem, recovering an sparse
error matrix $E\in\mathsf{F}^{r\times n}$ according to the equation $E=C-AB$.
In this section, we show how to treat each row of $E$ as a sparse polynomial
and evaluate each of these at consecutive powers of a high-order element
$\theta$, as needed for the sparse interpolation algorithm from Section 5. The
same techniques will also be used in the matrix inverse with errors algorithm.

Consider a column vector of powers of an indeterminate,
$\mathbf{x}=[1,x,x^{2},\ldots,x^{n-1}]^{T}$. The matrix-vector product
$E\mathbf{x}$ then consists of $m$ polynomials, each degree less than $n$, and
with a 1-1 correspondence between nonzero terms of the polynomials and nonzero
entries in the corresponding row of $E$.

Evaluating every polynomial in $E\mathbf{x}$ at a point $\theta\in\mathsf{F}$
is the same as multiplying $E$ by a vector
$\mathbf{v}=[1,\theta,\theta^{2},\ldots,\theta^{n-1}]^{T}$. Evaluating each
polynomial in $E\mathbf{x}$ at the first $2s$ powers
$1,\theta,\theta^{2},\ldots,\theta^{2s-1}$ allows for the unique recovery of
all $t$-sparse rows, according to Theorem 5.2. This means multiplying $E$
times an _evaluation matrix_ $V\in\mathsf{F}^{n\times 2s}$ such that
$V_{i,j}=\theta^{ij}$.

Consider a single row vector $\mathbf{u}\in\mathsf{F}^{1\times n}$ with $c$
nonzero entries. Multiplying $\mathbf{u}V$ means evaluating a $c$-sparse
polynomial at the first $2s$ powers of $\theta$; we can view it as removing
the $n-c$ zero entries from $\mathbf{u}$ to give
$\mathbf{u}^{\prime}\in\mathsf{F}^{1\times c}$, and removing the same $n-c$
rows from $V$ to get $V^{\prime}\in\mathsf{F}^{c\times 2s}$, and then
computing the product $\mathbf{u}^{\prime}V^{\prime}$. This evaluation matrix
with removed rows $V^{\prime}$ is actually a _transposed Vandermonde matrix_
built from the entries $\theta^{i}$ for each index $i$ of a nonzero entry in
$\mathbf{u}$.

An explicit algorithm for this transposed Vandermonde matrix application is
given in Section 6.2 of Bostan et al. [2003], and they show that the
complexity of multiplying $\mathbf{u}^{\prime}V^{\prime}$ is
$\widetilde{O}\mathopen{}\left(c+s\right)\mathclose{}$ field operations, and
essentially amounts to a power series inversion and product tree traversal.
(See also Section 5.1 of van der Hoeven and Lecerf [2013] for a nice
description of this approach to fast evaluation of a sparse polynomial at
consecutive powers.) We repeat this for each row in a given matrix $M$ to
compute $AV$.

Now we are ready to show how to compute $EV=(C-AB)V$ efficiently, using the
approach just discussed to compute $CV$ and $BV$ before multiplying $A$ times
$BV$ and performing a subtraction. As we will show, the rectangular
multiplication of $A$ times $BV$ is the only step which is not softly-linear
time, and this dominates the complexity.

Input: Matrices $A\in\mathsf{F}^{r\times\ell}$, $B\in\mathsf{F}^{\ell\times
n}$, $C\in\mathsf{F}^{r\times n}$, high-order element $\theta\in\mathsf{F}$,
and integer $s\in\mathbb{N}$

Output: Matrix $Y=(C-AB)V\in\mathsf{F}^{r\times 2s}$, where $V$ is the
$r\times 2s$ evaluation matrix given by $(\theta^{ij})_{0\leq i<n,0\leq
j<2s}$.

1 Pre-compute $\theta^{j}$ for all indices $j$ of nonzero columns in $B$ or
$C$

2 $Y_{C}\leftarrow CV$ using the pre-computed $\theta^{j}$’s and fast
transposed Vandermonde applications row by row

3 $Y_{B}\leftarrow BV$ similarly

4 $Y_{AB}\leftarrow AY_{B}$ using dense or sparse multiplication, whichever is
faster

return _$Y_{C}-Y_{AB}$_

Algorithm 4 DiffEval$(A,B,C,\theta,s)$

######  Lemma 6.1.

Algorithm 4 DiffEval always returns the correct matrix of evaluations
$(C-AB)V$ and uses

$\widetilde{O}\mathopen{}\left(\\#B+\\#C+n+s\cdot\min\left(\\#A+r,\frac{r\ell}{\min(r,s,\ell)^{3-\omega}}\right)\right)\mathclose{}$

field operations.

###### Proof.

Correctness comes from the discussion above and the correctness of the
algorithm in section 6.2 of Bostan et al. [2003] which is used on Algorithms 4
and 4.

The cost of pre-computing all powers $\theta^{j}$ using binary powering is
$O(n\log n)$ field operations, since the number of nonzero columns in $B$ or
$C$ is at most $n$. Using these precomputed values, each transposed
Vandermonde apply costs
$\widetilde{O}\mathopen{}\left(k+s\right)\mathclose{}$, where $k$ is the
number of nonzero entries in the current row. Summing over all rows of $C$ and
$B$ gives $\widetilde{O}\mathopen{}\left(\\#B+\\#C+rs+\ell
s\right)\mathclose{}$ for these steps.

The most significant step is the rectangular multiplication of
$A\in\mathsf{F}^{r\times\ell}$ times $Y_{B}\in\mathsf{F}^{\ell\times s}$.
Using sparse multiplication, the cost is $O\mathopen{}\left(\\#A\cdot
s\right)\mathclose{}$, by 3.3. Using dense multiplication, the cost is
$O\mathopen{}\left(r\ell s/\min(r,\ell,s)^{\omega-3}\right)\mathclose{}$
according to 3.1.

Because all the relevant parameters $\\#A,s,r,\ell$ are known before the
multiplication, we assume that the underlying software makes the best choice
among these two algorithms.

Note that $rs$ and $\ell s$ are definitely dominated by the cost of the dense
multiplication. In the sparse multiplication case, the situation is more
subtle when $\\#A<\ell$. But in this case, we just suppose that rows of $BV$
which are not used in the sparse multiplication $A(BV)$ are never computed, so
the complexity is as stated. ∎

##  7 Multiplication with error correction algorithm

We now have all the necessary components to present the complete
multiplication algorithm, Algorithm 5 MultiplyEC. The general idea is to
repeatedly use FindNonzeroRows to determine the locations on nonzeros in $E$,
then DiffEval and MultiSparseInterp to recover half of the nonzero rows at
each step.

Setting the target sparsity of each recovered row to $s=\left\lceil
2k/r\right\rceil$, where $r$ is the number of nonzero rows, ensures that at
least half of the rows are recovered at each step. To get around the fact that
the number of errors $k$ is initially unknown, we start with an initial guess
of $k=1$, and double this guess whenever fewer than half of the remaining rows
are recovered.

The procedure is probabilistic of the Monte Carlo type _only_ because of the
Monte Carlo subroutine FindNonzeroRows to determine which rows are erroneous.
The other parts of the algorithm — computing evaluations and recovering
nonzero entries — are deterministic based on the bounds they are given.

Input: Matrices $A\in\mathsf{F}^{m\times\ell}$, $B\in\mathsf{F}^{\ell\times
n}$, $C\in\mathsf{F}^{m\times n}$, high-order element $\theta\in\mathsf{F}$,
and error bound $0<\epsilon<1$

Output: Matrix $E\in\mathsf{F}^{m\times n}$ such that, with probability at
least $1-\epsilon$, $AB=C-E$

1 $k\leftarrow 1$

2 $E\leftarrow\mathbf{0}^{m\times n}$

3
$\epsilon^{\prime}\leftarrow\epsilon/(4\left\lceil\log_{2}(mn)\right\rceil+1)$

4
$\mathcal{J}\leftarrow\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto
CV-A(BV),\epsilon^{\prime})$

5 while  _$\\#\mathcal{J}\geq 1$_ do

6 if
_$\\#\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto(V^{T}(C-E)-(V^{T}A)B)^{T},\epsilon^{\prime})
>\\#\mathcal{J}$_ then

7 Transpose $C$ and $E$, swap and transpose $A$ and $B$, replace $\mathcal{J}$

8

9 $r\leftarrow\\#\mathcal{J}$

10 $A^{\prime}\leftarrow$ submatrix of $A$ from rows in $\mathcal{J}$

11 $C^{\prime}\leftarrow$ submatrix of $(C-E)$ from rows in $\mathcal{J}$

12 $s\leftarrow\left\lceil 2(k-\\#E)/r\right\rceil$

13
$Y\leftarrow\hyperlink{diffeval}{\texttt{DiffEval}{}}(A^{\prime},B,C^{\prime},\theta,s)$

14
$f_{1},\ldots,f_{r}\leftarrow\hyperlink{spinterp}{\texttt{MultiSparseInterp}{}}(r,n,s,\theta,Y)$

15 for  _$i\leftarrow 1,2,\ldots,r$_ do

16 Set $(\mathcal{J}_{i},e)$th entry of $E$ to $c$ for each term $cx^{e}$ of
$f_{i}$

17

18
$\mathcal{J}\leftarrow\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto(C-E)V-A(BV),\epsilon^{\prime})$

19 if  _$\\#\mathcal{J} >r/2$_ then

20 $k\leftarrow 2k$

21 if  _$k\geq 2n\\#J$_ then return _$C-AB$_

22

23 foreach  _$i\in\mathcal{J}$_ do

24 Clear entries from row $i$ of $E$ added on this iteration

25

return _$E$_

Algorithm 5 MultiplyEC$(A,B,C,\theta,\epsilon)$

######  Theorem 7.1.

With probability at least $1-\epsilon$, Algorithm 5 MultiplyEC finds all
errors in $C$ and uses

$\widetilde{O}\mathopen{}\left(m+\left\lceil\log_{\\#\mathsf{F}}\tfrac{1}{\epsilon}\right\rceil(n+t)+k\cdot\min\left(\left\lceil\frac{t}{r}\right\rceil,\,\ell/\min(\ell,r,\tfrac{k}{r})^{3-\omega}\right)\right)\mathclose{}$

field operations, where $k$ is the _actual_ number of errors in the given
product $C$. Otherwise, it uses $\widetilde{O}\mathopen{}\left(m\ell
n/\min(m,\ell,n)^{3-\omega}\right)\mathclose{}$ field operations and may
return an incorrect result.

###### Proof.

The only probabilistic parts of the algorithm are the calls to
FindNonzeroRows, which can only fail by incorrectly returning too few rows. If
this never happens, then each iteration through the while loop either (a)
discovers that the value of $k$ was too small and increases it, or (b)
recovers half of the rows or columns of $E$.

So the total number of iterations of the while loop if the calls to
FindNonzeroRows are never incorrect is at most
$\left\lceil\log_{2}k\right\rceil+\left\lceil\log_{2}m\right\rceil+\left\lceil\log_{2}n\right\rceil.$
By the union bound, from the way that $\epsilon^{\prime}$ is computed and the
fact that there are two calls to FindNonzeroRows on every iteration, the
probability that FindNonzeroRows is never incorrect is at least $1-\epsilon$.

In this case, the rest of the correctness and the running time follow directly
from Lemmas 4.1, 6.1 and 5.2.

If a call to FindNonzeroRows returns an incorrect result, two bad things can
happen. First, if $\mathcal{J}$ is too small during some iteration, all of the
evaluations may be incorrect, leading to $k$ being increased. In extremely
unlucky circumstances, this may happen $O(\log(rn))$ times until the product
is naïvely recomputed on Algorithm 5, leading to the worst-case running time.

The second bad thing that can occur when a call to FindNonzeroRows fails is
that it may incorrectly report all errors have been found, leading to the
small $1-\epsilon$ chance that the algorithm returns an incorrect result. ∎

##  8 Inverse error correction algorithm

As discussed in Section 2, our algorithm for correcting errors in a matrix
inverse follows mostly the same outline as that for correcting a matrix
product, with two important changes. These are the basis for the next two
lemmas.

######  Lemma 8.1.

For any matrices $E,A\in\mathsf{F}^{n\times n}$ where $A$ is invertible, a
call to $\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto
EAV,\epsilon)$ correctly returns the nonzero rows of $E$ with probability at
least $1-\epsilon$.

###### Proof.

Recall from the proof of Lemma 4.1 that the correctness of FindNonzeroRows
depends on applying a matrix $V\in\mathsf{F}^{n\times\ell}$ of random entries
to the unknown matrix $E$.

In the current formulation, instead of applying the random matrix $V$ directly
to $E$, instead the product $AV$ is applied. But because $A$ is nonsingular,
there is a 1-1 correspondence between the set of all matrices
$V\in\mathsf{F}^{n\times\ell}$ and matrices the set $\\{\,AV\mid
V\in\mathsf{F}^{n\times\ell}\,\\}$. That is, applying a random matrix to $EA$
is the same as applying a random matrix to $E$ itself, and therefore the same
results hold. ∎

######  Lemma 8.2.

Given any rank-$r$ matrix $A\in\mathsf{F}^{n\times r}$, it is possible to
compute a matrix $X\in\mathsf{F}^{r\times r}$ formed from a subset of the rows
of $A$, and its inverse $X^{-1}$, using
$\widetilde{O}\mathopen{}\left(\\#A+r^{\omega}\right)\mathclose{}$ field
operations.

###### Proof.

This is a direct consequence of [Cheung et al., 2013, Theorem 2.11], along
with the classic algorithm of [Bunch and Hopcroft, 1974] for fast matrix
inversion. Alternatively, one could use the approach of [Storjohann and Yang,
2015] to compute the lexicographically minimal set of linearly independent
rows in $M$, as well as a representation of the inverse, in the same running
time. ∎

The resulting algorithm for matrix inversion with errors is presented in
Algorithm 6 InverseEC.

Input: Matrices $A,B\in\mathsf{F}^{n\times n}$, high-order element
$\theta\in\mathsf{F}$, and error bound $0<\epsilon<1$

Output: Matrix $E\in\mathsf{F}^{n\times n}$ such that, with probability at
least $1-\epsilon$, $A^{-1}=B+E$

1 $k\leftarrow 1$

2 $E\leftarrow\mathbf{0}^{n\times n}$

3 $\epsilon^{\prime}\leftarrow\epsilon/(8\left\lceil\log_{2}n\right\rceil+1)$

4
$\mathcal{J}\leftarrow\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto
V-B(AV),\epsilon^{\prime})$

5 while  _$\\#\mathcal{J}\geq 1$_ do

6 if  _$\\#\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto
V-A((B+E)V),\epsilon^{\prime}) >\\#\mathcal{J}$_ then

7 Transpose $A$, $B$, and $E$, and replace $\mathcal{J}$

8

9 $r\leftarrow\\#\mathcal{J}$

10 $A^{\prime}\leftarrow$ submatrix of $A$ from columns in $\mathcal{J}$

11 $\mathcal{J}^{\prime}\leftarrow$ set of $r$ linearly independent rows in
$A^{\prime}$ according to Lemma 8.2

12 $I^{\prime},A^{\prime\prime}\leftarrow$ submatrices of $I,A$ from the rows
chosen in $\mathcal{J}^{\prime}$

13 $X^{-1}\leftarrow$ inverse of submatrix of $A^{\prime}$ according to Lemma
8.2

14 $s\leftarrow\left\lceil 2(k-\\#E)/r\right\rceil$

15
$Y\leftarrow\hyperlink{diffeval}{\texttt{DiffEval}{}}(A^{\prime\prime},B,I^{\prime},\theta,s)$

16 $Y^{\prime}\leftarrow X^{-1}Y$

17
$f_{1},\ldots,f_{r}\leftarrow\hyperlink{spinterp}{\texttt{MultiSparseInterp}{}}(r,n,s,\theta,Y^{\prime})$

18 for  _$i\leftarrow 1,2,\ldots,r$_ do

19 Set $(\mathcal{J}_{i},e)$th entry of $E$ to $c$ for each term $cx^{e}$ of
$f_{i}$

20

21
$\mathcal{J}\leftarrow\hyperlink{nzrows}{\texttt{FindNonzeroRows}{}}(V\mapsto
V-(B+E)(AV),\epsilon^{\prime})$

22 if  _$\\#\mathcal{J} >r/2$_ then

23 $k\leftarrow 2k$

24 if  _$k\geq 2n\\#J$_ then return _$A^{-1}-B$_

25

26 foreach  _$i\in\mathcal{J}$_ do

27 Clear entries from row $i$ of $E$ added on this iteration

28

return _$E$_

Algorithm 6 InverseEC$(A,B,\theta,\epsilon)$

######  Theorem 8.3.

With probability at least $1-\epsilon$, Algorithm 6 InverseEC finds all errors
in $B$ and uses

$\widetilde{O}\mathopen{}\left(\left\lceil\log_{\\#\mathsf{F}}\tfrac{1}{\epsilon}\right\rceil
t+nk/\min(r,\tfrac{k}{r})^{3-\omega}+r^{\omega}\right)\mathclose{}$

field operations, where $k$ is the _actual_ number of errors in the given
product $C$. Otherwise, it uses
$\widetilde{O}\mathopen{}\left(n^{\omega}\right)\mathclose{}$ field operations
and may return an incorrect result.

###### Proof.

In this algorithm, we make use of two different formulas for $E$:

$\displaystyle EA$ $\displaystyle=I-BA$ (8.1) $\displaystyle AE$
$\displaystyle=I-AB$ (8.2)

The first formula (8.1) is used to determine the nonzero rows of $E$ on
Algorithms 6 and 6, and the second formula (8.2) is used to evaluate the rows
of $XE$ on Algorithm 6.

The main difference in this algorithm compared to MultiplyEC is the need to
find $r$ nonzero rows of $A^{\prime}$ to constitute the matrix $X$ and compute
its inverse $X^{-1}$. According to Lemma 8.2 these both cost
$\widetilde{O}\mathopen{}\left(r^{\omega}\right)\mathclose{}$, which gives the
additional term in the complexity statement. Note that this also eliminates
the utility of sparse multiplication in all cases, simplifying the complexity
statement somewhat compared to that of MultiplyEC.

The rest of the proof is identical to that of Theorem 7.1. ∎

## References

  * Amossen and Pagh [2009] Rasmus Resen Amossen and Rasmus Pagh.  Faster join-projects and sparse matrix multiplications.  In _Proceedings of the 12th International Conference on Database Theory_ , ICDT ’09, pages 121–126, New York, NY, USA, 2009. ACM.  doi: 10.1145/1514894.1514909.
  * Arnold et al. [2015] Andrew Arnold, Mark Giesbrecht, and Daniel S. Roche.  Faster sparse multivariate polynomial interpolation of straight-line programs.  _Journal of Symbolic Computation_ , 2015.  doi: 10.1016/j.jsc.2015.11.005.
  * Ben-Or and Tiwari [1988] Michael Ben-Or and Prasoon Tiwari.  A deterministic algorithm for sparse multivariate polynomial interpolation.  In _Proceedings of the twentieth annual ACM symposium on Theory of computing_ , STOC ’88, pages 301–309, New York, NY, USA, 1988. ACM.  doi: 10.1145/62212.62241.
  * Bernstein [2005] Daniel J. Bernstein.  Factoring into coprimes in essentially linear time.  _Journal of Algorithms_ , 54(1):1–30, 2005.  doi: 10.1016/j.jalgor.2004.04.009.
  * Böhm et al. [2015] Janko Böhm, Wolfram Decker, Claus Fieker, and Gerhard Pfister.  The use of bad primes in rational reconstruction.  _Math. Comp._ , 84(296):3013–3027, 2015.  doi: 10.1090/mcom/2951.
  * Borodin and Munro [1975] A. Borodin and I. Munro.  _The computational complexity of algebraic and numeric problems_.  Number 1 in Elsevier Computer Science Library; Theory of Computation Series. American Elsevier Pub. Co., New York, 1975.
  * Bostan et al. [2003] A. Bostan, G. Lecerf, and É. Schost.  Tellegen’s principle into practice.  In _Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’03, pages 37–44. ACM, 2003.  doi: 10.1145/860854.860870.
  * Boyer and Kaltofen [2014] Brice Boyer and Erich L. Kaltofen.  Numerical linear system solving with parametric entries by error correction.  In _Proceedings of the 2014 Symposium on Symbolic-Numeric Computation_ , SNC ’14, pages 33–38. ACM, 2014.  doi: 10.1145/2631948.2631956.
  * Bunch and Hopcroft [1974] James R. Bunch and John E. Hopcroft.  Triangular factorization and inversion by fast matrix multiplication.  _Mathematics of Computation_ , 28(125):231–236, 1974.  URL http://www.jstor.org/stable/2005828.
  * Cantor and Kaltofen [1991] David G. Cantor and Erich Kaltofen.  On fast multiplication of polynomials over arbitrary algebras.  _Acta Informatica_ , 28:693–701, 1991.  doi: 10.1007/BF01178683.
  * Cheung et al. [2013] Ho Yee Cheung, Tsz Chiu Kwok, and Lap Chi Lau.  Fast matrix rank algorithms and applications.  _J. ACM_ , 60(5):31:1–31:25, October 2013.  doi: 10.1145/2528404.
  * Comer et al. [2012] Matthew T. Comer, Erich L. Kaltofen, and Clément Pernet.  Sparse polynomial interpolation and Berlekamp/Massey algorithms that correct outlier errors in input values.  In _Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’12, pages 138–145, New York, NY, USA, 2012\. ACM.  doi: 10.1145/2442829.2442852.
  * Dumas and Kaltofen [2014] Jean-Guillaume Dumas and Erich Kaltofen.  Essentially optimal interactive certificates in linear algebra.  In _Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’14, pages 146–153. ACM, 2014.  doi: 10.1145/2608628.2608644.
  * Dumas et al. [2008] Jean-Guillaume Dumas, Pascal Giorgi, and Clément Pernet.  Dense linear algebra over word-size prime fields: the FFLAS and FFPACK packages.  _ACM Trans. Math. Softw._ , 35:19:1–19:42, October 2008\.  doi: 10.1145/1391989.1391992.
  * Dumas et al. [2016] Jean-Guillaume Dumas, Erich Kaltofen, Emmanuel Thomé, and Gilles Villard.  Linear time interactive certificates for the minimal polynomial and the determinant of a sparse matrix.  In _Proceedings of the ACM on International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’16, pages 199–206. ACM, 2016.  doi: 10.1145/2930889.2930908.
  * Dumas et al. [2017] Jean-Guillaume Dumas, David Lucas, and Clément Pernet.  Certificates for triangular equivalence and rank profiles.  In _Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’17, pages 133–140. ACM, 2017.  doi: 10.1145/3087604.3087609.
  * Freivalds [1979] Rūsiņš Freivalds.  Fast probabilistic algorithms.  In Jiří Bečvář, editor, _Mathematical Foundations of Computer Science 1979_ , pages 57–69. Springer Berlin Heidelberg, 1979.
  * Gasieniec et al. [2017] Leszek Gasieniec, Christos Levcopoulos, Andrzej Lingas, Rasmus Pagh, and Takeshi Tokuyama.  Efficiently correcting matrix products.  _Algorithmica_ , 79(2):428–443, Oct 2017.  doi: 10.1007/s00453-016-0202-3.
  * Goldreich et al. [1999] Oded Goldreich, Dana Ron, and Madhu Sudan.  Chinese remaindering with errors.  In _Proceedings of the Thirty-first Annual ACM Symposium on Theory of Computing_ , STOC ’99, pages 225–234. ACM, 1999.  doi: 10.1145/301250.301309.
  * Harvey et al. [2017] David Harvey, Joris van der Hoeven, and Grégoire Lecerf.  Faster polynomial multiplication over finite fields.  _J. ACM_ , 63(6):52:1–52:23, January 2017.  doi: 10.1145/3005344.
  * van der Hoeven and Lecerf [2013] Joris van der Hoeven and Grégoire Lecerf.  On the bit-complexity of sparse polynomial and series multiplication.  _Journal of Symbolic Computation_ , 50:227–0254, 2013.  doi: 10.1016/j.jsc.2012.06.004.
  * Huang and Gao [2017] Qiao-Long Huang and Xiao-Shan Gao.  Faster deterministic sparse interpolation algorithms for straight-line program multivariate polynomials.  _CoRR_ , abs/1709.08979, 2017.  URL http://arxiv.org/abs/1709.08979.
  * Javadi and Monagan [2010] Seyed Mohammad Mahdi Javadi and Michael Monagan.  Parallel sparse polynomial interpolation over finite fields.  In _Proceedings of the 4th International Workshop on Parallel and Symbolic Computation_ , PASCO ’10, pages 160–168, New York, NY, USA, 2010\. ACM.  doi: 10.1145/1837210.1837233.
  * Kaltofen and Yagati [1989] Erich Kaltofen and Lakshman Yagati.  Improved sparse multivariate polynomial interpolation algorithms.  In P. Gianni, editor, _Symbolic and Algebraic Computation_ , volume 358 of _Lecture Notes in Computer Science_ , pages 467–474. Springer Berlin / Heidelberg, 1989.  doi: 10.1007/3-540-51084-2_44.
  * Kaltofen [2010] Erich L. Kaltofen.  Fifteen years after DSC and WLSS2: What parallel computations I do today [invited lecture at PASCO 2010].  In _Proceedings of the 4th International Workshop on Parallel and Symbolic Computation_ , PASCO ’10, pages 10–17, New York, NY, USA, 2010. ACM.  doi: 10.1145/1837210.1837213.
  * Kaltofen and Yang [2013] Erich L. Kaltofen and Zhengfeng Yang.  Sparse multivariate function recovery from values with noise and outlier errors.  In _Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’13, pages 219–226. ACM, 2013.  doi: 10.1145/2465506.2465524.
  * Kaltofen et al. [2011] Erich L. Kaltofen, Michael Nehring, and B. David Saunders.  Quadratic-time certificates in linear algebra.  In _Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’11, pages 171–176. ACM, 2011.  doi: 10.1145/1993886.1993915.
  * Kaltofen et al. [2017] Erich L. Kaltofen, Clément Pernet, Arne Storjohann, and Cleveland Waddell.  Early termination in parametric linear system solving and rational function vector recovery with error correction.  In _Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’17, pages 237–244. ACM, 2017.  doi: 10.1145/3087604.3087645.
  * Khonji et al. [2010] Majid Khonji, Clément Pernet, Jean-Louis Roch, Thomas Roche, and Thomas Stalinski.  Output-sensitive decoding for redundant residue systems.  In _Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’10, pages 265–272. ACM, 2010.  doi: 10.1145/1837934.1837985.
  * Le Gall [2012] François Le Gall.  Faster algorithms for rectangular matrix multiplication.  In _2012 IEEE 53rd Annual Symposium on Foundations of Computer Science_ , pages 514–523, Oct 2012.  doi: 10.1109/FOCS.2012.80.
  * Le Gall [2014] François Le Gall.  Powers of tensors and fast matrix multiplication.  In _Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’14, pages 296–303, New York, NY, USA, 2014\. ACM.  doi: 10.1145/2608628.2608664.
  * Lingas [2009] Andrzej Lingas.  A fast output-sensitive algorithm for boolean matrix multiplication.  In Amos Fiat and Peter Sanders, editors, _Algorithms - ESA 2009_ , pages 408–419. Springer Berlin Heidelberg, 2009.  doi: 10.1007/978-3-642-04128-0_37.
  * Pagh [2013] Rasmus Pagh.  Compressed matrix multiplication.  _ACM Trans. Comput. Theory_ , 5(3):9:1–9:17, August 2013.  doi: 10.1145/2493252.2493254.
  * Storjohann and Yang [2015] Arne Storjohann and Shiyun Yang.  A relaxed algorithm for online matrix inversion.  In _Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’15, pages 339–346, New York, NY, USA, 2015. ACM.  doi: 10.1145/2755996.2756672.
  * Wang et al. [2013] Qian Wang, Xianyi Zhang, Yunquan Zhang, and Qing Yi.  Augem: Automatically generate high performance dense linear algebra kernels on x86 cpus.  In _Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_ , SC ’13, pages 25:1–25:12. ACM, 2013.  doi: 10.1145/2503210.2503219.
  * Yuster and Zwick [2005] Raphael Yuster and Uri Zwick.  Fast sparse matrix multiplication.  _ACM Trans. Algorithms_ , 1(1):2–13, July 2005\.  doi: 10.1145/1077464.1077466.