---
author:
- 'Mohammad Fahim and Viveck R. Cadambe[^1][^2]'
bibliography:
- 'IEEEabrv.bib'
- 'sample.bib'
title: Numerically Stable Polynomially Coded Computing
---

Introduction
============

The recently emerging area of “coded computing” focuses on incorporating redundancy based on coding-theory-inspired strategies to tackle central challenges in distributed computing, including stragglers, failures, processing errors, communication bottlenecks and security issues. Such ideas have been applied to different large scale distributed computations such as matrix multiplication [@dutta2016short; @polynomialcodes; @allerton17; @arxiv_allerton17; @genPolyDot], gradient methods [@tandon2016gradient; @tan17; @abbegrad], linear solvers [@YangGK17; @crit; @maity2018robust] and multi-variate polynomial evaluation [@yu2018lagrange]. An important idea that has emerged from this body of the work is the use of novel, Reed-Solomon like *polynomial* based methods for encoding data. In polynomial based methods, each computation node stores a linearly encoded combination of the data partitions, where data stored at different worker nodes can be interpreted as evaluation of an appropriate polynomial at different points. The nodes then perform computation on these encoded versions of the data, and a central master/fusion node aggregates the outputs of these computations to recover the overall result via a decoding process that inevitably involves polynomial interpolation. Much like Reed Solomon Codes, if the number of nodes performing the computation is higher than the number of evaluation points required for accurate interpolation, the overall computation is tolerant to faults and stragglers.

Perhaps the most striking application of polynomial based methods comes in the context of matrix multiplication. To multiply two $N \times N$ matrices $\mathbf{A},\mathbf{B},$ assuming that each node stores $1/m$ of each matrix, classical work in algorithm based fault tolerance [@Huang_TC_84] outlines a coding based method which has been analyzed in [@ProductCodes]. Reference [@polynomialcodes] showed through *polynomial* based encoding methods that the result of just $m^2$ nodes can be used by the master node to recover the matrix-product. Remarkably, this means that polynomial based codes ensure that *the recovery threshold* - the worst case number of nodes whose computation suffices to recover the overall matrix-product - does not grow with $P$, the number of the distributed system’s worker nodes, unlike the approaches of [@Huang_TC_84; @ProductCodes]. The recovery threshold for matrix multiplication has been improved to $2m-1$ via a code construction called MatDot Codes in [@allerton17], albeit at a higher communication/computation cost than codes in [@polynomialcodes]. A second prominent application of polynomial based methods is the idea of *Lagrange coded computing* [@yu2018lagrange], where coding is applied for multi-variate polynomial computing with guarantees of straggler resilience, security and privacy. In addition, polynomial-based methods are also useful for communication-efficient approaches for inverse problems and gradient methods [@abbegrad; @li2018polynomially; @crit].

Despite the enormous success, the scalability of polynomial based methods in practice are limited by an “inconvenient truth”, their numerical instability. The decoding methods for polynomial based methods require interpolating a degree $K-1$ polynomial using $K$ evaluation points. While this is numerically stable for classical error correcting codes for communication and storage which are implemented over finite fields, we are concerned here for data processing applications where the operations are typically real-valued. The main reason for the instability is that either implicitly or explicitly, interpolation effectively solves a linear system whose transform is characterized by a Vandermonde matrix. It is well known that the condition number of Vandermonde matrices with real-valued nodes grows exponentially in the dimension of the matrix [@gautschi1987lower; @gautschi1990stable; @gautschi1974norm; @reichel1991chebyshev]. The large condition number means that small perturbations of the Vandermonde matrix due to numerical precision errors can result in singular matrices [@quarteroni2010numerical; @trefethen2013approximation]. In practice, this can translate to large numerical errors even when the coded computation is distributed among few tens of nodes[^3].

Summary of Contributions
========================

In this paper, we develop a new, numerically stable, approach for polynomially coded computing. A significant difference from previous polynomial coding approaches is that we depart from the monomial basis, which allows us to circumvent the inherently ill-conditioned Vandermonde-matrices. We demonstrate our approach through two important applications of polynomially coded computing: matrix multiplication, and Lagrange coded computing.

To illustrate our results, consider the coded matrix multiplication problem, where the goal is to multiply two matrices $\mathbf{A},\mathbf{B}$ over $P$ computation nodes where each node stores $1/m$ of each of the two matrices. A master node encodes $\mathbf{A},\mathbf{B}$ into $P$ matrices each, and sends these matrices respectively to each worker node. Each worker node multiplies the received encoded matrices, and sends the product back to the fusion node[^4], which aims to recover $\mathbf{A}\mathbf{B}$ from a subset of the worker nodes. The recovery threshold is defined as a number $K$ such that the computation of any set of $K$ worker nodes suffices to recover the product $\mathbf{A}\mathbf{B}.$ The MatDot scheme of [@allerton17] achieves the best known recovery threshold of $2m-1$. We begin with an example of MatDot Codes for $m=2.$

**Example 1: MatDot Codes [@allerton17], recovery threshold = 3:** *Consider two $N \times N$ matrices* $$\mathbf{A}=\begin{bmatrix}\mathbf{{A}}_1 & \mathbf{{A}}_{2}\end{bmatrix},~~ \mathbf{B} = \begin{bmatrix}\mathbf{{B}}_1 \\ \mathbf{{B}}_{2}\end{bmatrix},$$ where $\mathbf{A}_{1},\mathbf{A}_{2}$ are $N \times N/2$ matrices and $\mathbf{B}_{1},\mathbf{B}_{2}$ are $N/2 \times N$ matrices. *Define $p_{\mathbf{A}}(x) = \mathbf{A}_{1} + \mathbf{A}_{2}x$ and $p_{\mathbf{B}}(x) = \mathbf{B}_{1}x+\mathbf{B}_{2},$ and let $x_1, \cdots, x_P$ be distinct real values. Notice that $\mathbf{AB}=\mathbf{A}_{1}\mathbf{B}_1+\mathbf{A}_2\mathbf{B}_2$ is the coefficient of $x$ in polynomial $p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$. In MatDot Codes, as illustrated in Fig. \[fig:Ex1\], worker node $i$ computes $p_{\mathbf{A}}(x_i)p_{\mathbf{B}}(x_i),~ i=1,2, \ldots P,$ so that from any $3$ of the $P$ nodes, the polynomial $p(x) = \mathbf{{A}}_{1}\mathbf{{B}}_2+(\mathbf{{A}}_{1}\mathbf{{B}}_{1} + \mathbf{{A}}_{2}\mathbf{{B}}_{2}) x + \mathbf{{A}}_2\mathbf{{B}}_1 x^{2}$ can be interpolated. Having interpolated the polynomial, the product $\mathbf{A}\mathbf{B}$ is simply the coefficient of $x$.*

![Example of MatDot Codes [@allerton17], with a recovery threshold of $3$. The matrix product $\mathbf{A}\mathbf{B}$ is the coefficient of $x$ in $p_\mathbf{A}(x)p_\mathbf{B}(x)$, and can be recovered at the fusion node upon receiving the output of any $3$ worker nodes and interpolating $p_\mathbf{A}(x)p_\mathbf{B}(x)$. []{data-label="fig:Ex1"}](Ex1.pdf)

A generalization of the above example leads to a recovery threshold of $2m-1$, with a decoding process that involves effectively inverting a $2m-1 \times 2m-1$ Vandermonde matrix. It has been shown that the condition number of the $n \times n$ Vandermonde matrix grows exponentially in $n$ with both $\ell_{\infty}$ and $\ell_{2}$ norms [@gautschi1987lower; @gautschi1990stable]. The intuition behind the inherent poor conditioning of the monomial basis $\{1,x,x^{2},\ldots, x^{2m-1}\}$ is demonstrated in Fig. \[fig:monomialplot\] and Fig. \[fig:monomialvec\].

![image](polys_mon-eps-converted-to.pdf){width="\textwidth"}

![image](Angle_poly.pdf){width="\textwidth"}

![image](polys_cheb-eps-converted-to.pdf){width="\textwidth"}

Motivated by Fig.\[fig:monomialvec\], we aim, in this paper, to choose polynomials that are orthonormal. However, it is not immediately clear whether orthonormal polynomials are applicable for matrix multiplications. We demonstrate the applicability of orthonormal codes for matrix multiplication. For the example below, let $q_0(x), q_1(x)$ denote two orthonormal polynomials such that $$\int_{-1}^{1} q_i(x) q_j(x) dx = \left\{\begin{array}{cc} 0 & \textrm{if }i=j \\ 1 & \textrm{otherwise} \end{array}\right.\label{eq:ortho}$$where $q_i(x),i=0,1$ has degree $i$.

**Example 2 : , recovery threshold = 3:** *For two $N \times N$ matrices $\mathbf{A}=\begin{bmatrix}\mathbf{{A}}_1 & \mathbf{{A}}_{2}\end{bmatrix}, \mathbf{B} = \begin{bmatrix}\mathbf{{B}}_1 \\ \mathbf{{B}}_{2}\end{bmatrix},$ let $p_{\mathbf{A}}(x) = \mathbf{A}_{1} q_0(x) + \mathbf{A}_{2}q_1(x)$ and $p_{\mathbf{B}}(x) = \mathbf{B}_{1}q_0(x)+\mathbf{B}_{2} q_1(x).$ Notice that because of (\[eq:ortho\]), we have $$\mathbf{A}\mathbf{B} = \int_{-1}^{1} p_{\mathbf{A}}(x) p_{\mathbf{B}}(x) dx.$$ This leads to the following coded computing scheme: worker node $i$ computes $p_{\mathbf{A}}(x_i)p_{\mathbf{B}}(x_i),~ i=1,2, \ldots P,$ where $x_1,\cdots, x_P$ are distinct real values, so that from any $3$ of the $P$ nodes, the fusion node can interpolate $p(x) = p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$. Having interpolated the polynomial, the fusion node obtains the product $\mathbf{A}\mathbf{B}$ by performing $\int_{-1}^{1}p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)dx$. This example is illustrated in Fig. \[fig:Ex2\]*.

![Example of our proposed orthonormal polynomials based codes, with a recovery threshold of $3$. The matrix product $\mathbf{A}\mathbf{B}$ is $\int_{-1}^{1} p_\mathbf{A}(x)p_\mathbf{B}(x)dx $, and can be recovered at the fusion node upon receiving the output of any $3$ worker nodes, then interpolating $p_\mathbf{A}(x)p_\mathbf{B}(x)$, and computing the integral $\int_{-1}^{1} p_\mathbf{A}(x)p_\mathbf{B}(x)dx.$[]{data-label="fig:Ex2"}](Ex2.pdf)

A simple generalization of the above example, described in Construction \[con:orthpoly\] in Section \[sec:ortho\], leads to a class of codes, with recovery threshold of $2m-1$, the same recovery threshold as MatDot Codes. In general, orthonormal polynomials are defined over arbitrary weight measure $\int_{-1}^{1} {\bf \cdot}~ w(x) dx;$ some well known classes of polynomials corresponding to different weight measures $w(x)$ include Legendre, Chebyshev, Jacobi and Laguerre Polynomials [@quarteroni2010numerical; @trefethen2013approximation] (See Section \[sec:prem\] for definitions). Our Codes in Section \[sec:ortho\] can use any weight measure, and therefore can be used with different classes of orthonormal polynomials. Of particular interest to our paper are the Chebyshev polynomials (Fig. \[fig:chebplot\]).

With our basic template, the task of developing numerically stable codes boils down to (A) interpolating $p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$ in a numerically stable manner, and (B) integrating this polynomial in a numerically stable manner. For task (B), we use a decoding procedure via *Gauss Quadrature* [@quarteroni2010numerical; @num97; @trefethen2013approximation] to recover the integral. Task (A) is particularly challenging in the coding setting, because our goal is to interpolate the coefficients of $p_\mathbf{A}(x)p_{\mathbf{B}}(x)$ - expanded over a series of orthonormal polynomials - from any $2m-1$ points among a set of $P$ points.

In Section \[sec:Chebyshev\], we provide a , a numerically stable matrix multiplication code construction . The construction specializes via the use of Chebyshev polynomials, which are a class of orthogonal polynomials that are ubiquitous in numerical methods and approximation theory [@trefethen2013approximation]. Construction \[con:chebpoly\] also specifies the choice of evaluation points $x_1, x_2, \ldots, x_P.$

The decoding procedure outlined involves the effective inversion of some $2m-1 \times 2m-1$ sub-matrix of a $2m-1 \times P$ *Chebyshev-Vandermonde* matrix [@reichel1991chebyshev], where each of the $i$-th column contains evaluations of the first $2m-1$ Chebyshev polynomials at $x_i,i=1,2,\ldots,P$. A key technical result of our paper shows that, with our choice of evaluation points $x_1, x_2, \ldots, x_P,$ every $2m-1 \times 2m-1$ square sub-matrix of the $2m-1 \times P$ Chebyshev-Vandermonde matrix is well-conditioned. More precisely, we show that, with our choice of $x_1, x_2, \ldots, x_P$, the condition number of *any* $2m-1 \times 2m-1$ sub-matrix of the Chebyshev-Vandermonde matrix grows at most polynomially in $P$ when the number of redundant parity nodes $\Delta=P-(2m-1)$ is fixed. Our condition number bound may be viewed as result of independent interest in the area of numerical methods, and requires non-trivial use of techniques from numerical approximation theory. This result is in contrast with the well known exponential growth for Vandermonde systems. We also show the significant improvement in stability via numerical experiments in Section \[sec:numerical\]. We also provide a preview of the results here in Table \[table:errors\], whose results demonstrate that remarkably, our Chebyhev-Vandermonde construction with even $P=150$ nodes has a smaller relative error than the Vandermonde-based MatDot Codes[^5] with $P=30$ nodes.

  ------------ ---------------------- ----------------------- ------------------------ ------------------------
     Number          [MatDot]{}                                   [[MatDot ]{}]{}      
   of Workers        worst case             worst case                average                  average
     $(P)$         relative error         relative error           relative error           relative error
       30       $1.54\times 10^{-6}$   $5.14\times 10^{-11}$   $1.36 \times 10^{-7}$    $1.36 \times 10^{-13}$
       50       $8.6 \times 10^{3}$    $1.27 \times 10^{-9}$     $2.00 \times 10^2$     $2.04 \times 10^{-13}$
       80       $2.45 \times 10^{6}$   $1.98 \times 10^{-8}$   $ 2.19 \times 10^{2} $   $3.08 \times 10^{-12}$
      150       $3.87 \times 10^{7}$   $7.84 \times 10^{-7}$    $8.73 \times 10^{2}$    $2.03 \times 10^{-11}$
  ------------ ---------------------- ----------------------- ------------------------ ------------------------

While MatDot Codes [@allerton17] have an optimal recovery threshold of $2m-1$, they have relatively higher computation cost per worker ($O(N^3/m)$) and worker node to fusion node communication cost ($O(N^2)$) as compared to Polynomial Codes [@polynomialcodes] which have a computation cost per worker of $O(N^3/m^2)$ and worker node to fusion node communication cost of $O(N^2/m^2)$. In particular, each worker in MatDot Codes performs an “outer” product of an $N \times N/m$ matrix with a $N/m \times N$ matrix, whereas each worker in Polynomial Codes performs an “inner” product of a $N/m \times N$ matrix with a $N \times N/m$ matrix. The reduced computation/communication comes at the cost of weaker fault-tolerance - Polynomial Codes have a higher recovery threshold of $m^2$ as compared with MatDot Codes ($2m-1$). In Section \[sec:chebpolycod\], we develop numerically stable codes for matrix multiplication, again via orthogonal polynomials, that achieve the same low computation/communication costs as Polynomial Codes as well as the same recovery threshold, .

The trade-off between computation/communication cost and recovery threshold imposed by MatDot Codes and Polynomial Codes has motivated general code constructions that interpolates both of them [@allerton17; @genPolyDot; @entPoly], albeit using the monomial basis. In Section \[sec:polydot-chebyshev\], we extend our approach to a general matrix multiplication code construction, that offers a computation/communication cost vs recovery threshold trade-off, following the research thread for the monomial basis [@allerton17; @genPolyDot; @entPoly], however we also target numerical stability in our proposed construction. While our , i.e., they achieve the same optimal recovery threshold as Codes when allowing for the same computation/communication cost as Codes, they do not specialize to Codes. Specifically, codes have higher recovery threshold than Codes when allowing for the same computation/communication cost as OrthoPoly Codes. In Section \[sec:lcc\], we exploit the result obtained in Theorem \[thm:bound\] on the condition number of the square $K \times K$ sub-matrices of the $K\times P$ Chebyshev-Vandermonde matrices to propose a numerically stable algorithm for Lagrange coded computing. In Section \[sec:conc\], we conclude with a discussion on other related problems such as matrix-vector multiplication [@Huang_TC_84; @lee2016speeding], and describe some related open questions.

Preliminaries on Numerical Analysis and Notations {#sec:prem}
=================================================

[ We discuss, in this section, the problem of finite precision in representing real numbers on digital machines and how it may horribly affect the output of computation problems performed on these machines. In addition, we also introduce some basic definitions and results from the area of numerical approximation theory that will be used in this paper[@num97], [@shortApp]. At the end of this section, we provide most of the common notations that will be used in this paper. ]{}

Preliminaries on Numerical Analysis
-----------------------------------

Since digital machines have finite memory, real numbers are digitally stored using a finite number of bits, i.e., finite precision. However, storing real numbers using a finite number of bits leads to inevitable errors since a finite number of bits can only represent a finite number of real numbers with no errors. On the other hand, real numbers that cannot be directly represented using the specified finite number of bits have to be either truncated or rounded-off in order to fit in the memory. Although such perturbation (e.g., truncation/round-off error) of real numbers due to the finite precision of digital machines can be negligibly small, the perturbation of the output of any computation that uses such “small" perturbed stored real numbers as input does not necessarily be small as well. In fact, a very small perturbation to the input of some computation may lead to an output that is totally wrong and irrelevant to the correct output. The condition number of a computation problem captures/measures this observation.

Let $f$ be a function representing a computation problem with input $x$, and let $\delta x$ be a small perturbation of $x$, and define $\delta f=f(x+\delta x)-f(x)$ to be the perturbation of $f$ due to $\delta x$, the condition number of the problem at $x$ with respect to some norm $||\cdot||$ is $$\begin{aligned}
\kappa(x)=\sup_{\delta x} \left(\frac{||\delta f||}{||f(x)||}\middle/\frac{||\delta x||}{||x||}\right).\end{aligned}$$

Given the above definition of condition number, a problem is said to be “ill-conditioned" if small perturbations in the input lead to large perturbation in the output (i.e., the condition number is large). On the other hand, a problem is said to be “well-conditioned" if small perturbations in the input lead to small perturbations in the output (i.e., the condition number is small).

[ In what follows, we discuss the condition number of two computation problems: the matrix-vector multiplication and solving a system of linear equations. For both problems, consider the system of linear equations represented in the matrix form $\mathbf{A}\mathbf{x}=\mathbf{y}$, where $\mathbf{A}\in \mathbb{R}^{n,n}$ and non-singular, and $\mathbf{x},\mathbf{y}\in\mathbb{R}^n$, and let $||\cdot||$ be some matrix norm. Then, let $\mathbf{A}$ be fixed, the condition number of this matrix-vector multiplication problem with $\mathbf{y}$ as its output given small perturbations in the input $\mathbf{x}$ is $\kappa(\mathbf{x})\leq ||\mathbf{A}||||\mathbf{A}^{-1}||$, for any $\mathbf{x} \in \mathbb{R}^n$. Also, for the problem of solving the system of linear equations $\mathbf{A}\mathbf{x}=\mathbf{y}$, with $\mathbf{A}$ still fixed, the condition number of the problem of solving this system of linear equations, given small perturbations in the input $\mathbf{y}$, where $\mathbf{x}$ is the output, is $\kappa(\mathbf{y}) \leq ||\mathbf{A}||||\mathbf{A}^{-1}||$, for any $\mathbf{y}  \in \mathbb{R}^n$. ]{}

Since we focus on polynomially coded computing, next, we introduce some basic tools of numerical approximation theory that will be used throughout this paper. Notice that, in the following, $C[a,b]$ denotes the vector space of continuous integrable functions defined on the interval $[a,b]$.

For any $f,g \in C[a,b]$, and given a non-negative integrable weight function $w$, $$\begin{aligned}
\langle f,g \rangle=\int_{a}^{b} f(x)g(x)w(x) dx\notag\end{aligned}$$ defines an inner product on $C[a,b]$ relative to $w$.

Consider a non-negative integrable weight function $w$, the polynomials $\{q_i\}_{i\geq 0}$ in $C[a,b]$ and $$\begin{aligned}
    \langle q_i, q_j \rangle= \left\{ \begin{array}{ll}
       c_{i} & \text{if $i=j$,} \\
        0 &\text{otherwise,}
        \end{array}\right.\end{aligned}$$ for some non-zero values $c_{i}$, where the inner product is relative to $w$, are called orthogonal polynomials relative to $w$, .

Consider a non-negative integrable weight function $w$, the polynomials $\{q_i\}_{i\geq 0}$, , in $C[a,b]$ such that $$\begin{aligned}
    \langle q_i, q_j \rangle= \left\{ \begin{array}{ll}
       1 & \text{if $i=j$,} \\
        0 &\text{otherwise,}
        \end{array}\right.\end{aligned}$$ where the inner product is relative to $w$, are called orthonormal polynomials relative to $w$.

\[ex:cheb\] The following recurrence relation defines the Chebyshev polynomials of the first kind: $$T_n(x)=2xT_{n-1}(x)-T_{n-2}(x),$$ where, $T_0(x)=1, T_1(x)=x$. These Chebyshev polynomials are the corner stone of modern numerical approximation theory and practice with applications to numerical integration, and least-square approximations of continuous functions [@num97],[@shortApp]. $\frac{1}{\sqrt{2}} T_0, T_1, T_2, \cdots$ are orthonormal relative to the weight function $\frac{2}{\pi \sqrt{1-x^2}}$. In general, Chebyshev polynomials are defined over $x\in \mathbb{R}$. However, for $x\in [-1,1]$, $T_n(x) = \cos(n\arccos(x))$, for any $n \in \mathbb{N}$. For the rest of this paper, unless otherwise is stated, whenever Chebyshev polynomials are used, they are restricted only to the range $[-1,1]$.

Let $w$ be a weight function on the range $[a,b]$, , and let $x_1, \cdots, x_n$ be distinct real numbers such that $a < x_1 < \cdots < x_n < b$, there exist unique weights $a_1, \cdots, a_n$ such that $$\int_a^b f(x) w(x) dx =\sum_{i=1}^n a_i f(x_i),$$ for all polynomials $f$ with degree less than $n$. \[thm:integration\]

\[thm:GQ\] Fix a weight function $w$, and let $\{q_i\}_{i\geq 0}$ be a set of orthonormal polynomials in $C[a,b]$ relative to $w$. Given $n$, let $\eta_1, \cdots, \eta_n$ be the roots of $q_n$ such that $a\leq \eta_1< \eta_2 < \cdots <\eta_n\leq b$, and choose real values $a_1, \cdots, a_n$ such that $\sum_{i=1}^{n} a_i f(\eta_i) = \int_{a}^{b} f(x) w(x) dx$, for any $f \in C[a,b]$ with degree less than $n$. Then, $\sum_{i=1}^{n} a_i f(\eta_i) = \int_{a}^{b} f(x) w(x) dx$, for any polynomial $f$ with degree less than $2n$.

\[rmk:orth\]

1.  Consider any orthonormal polynomials $\{q_i\}_{i>0}$. For any $n \in \mathbb{N}$, the set $\{q_0, q_1, \cdots, q_{n-1}\}$ forms a basis for the vector space of polynomials with degree less than $n$.

2.  In Theorem \[thm:GQ\], $a_1, \cdots, a_n$ can be chosen as $$\begin{aligned}
    \label{eq:ais}
    a_i= \int_{a}^{b} \bigg(\prod_{j\in[n]-i} \frac{x-\eta_j}{\eta_i-\eta_j}\bigg) w(x)dx ,~ i\in[n].\end{aligned}$$

3.  In Theorem \[thm:GQ\], the roots of $q_n$, i.e., $\eta_1, \cdots ,\eta_n$ are, in fact, real and distinct. Moreover, the Chebyshev polynomial of the first kind $T_n$ has the following roots $$\begin{aligned}
    \label{eq:Chebnodes}
     \rho^{(n)}_i =  \cos\left(\frac{2i-1}{2n} \pi\right),~ i\in [n].\end{aligned}$$ The set $\{\rho^{(n)}_1, \cdots, \rho^{(n)}_n\}$ is often called the $n$-point Chebyshev grid, and its elements $\rho^{(n)}_1, \cdots, \rho^{(n)}_n$ are called “Chebyshev nodes” of degree $n$. We here discard the term “node” and use the term “Chebyshev points” to avoid confusion with computation nodes. We also denote by $\boldsymbol{\rho}^{(n)}$ the vector $(\rho_1^{(n)}, \cdots, \rho_n^{(n)})$. It is useful to note that $T_n(x)$ can be written as $$\begin{aligned}
    T_n(x)=2^{n-1}\prod_{i=1}^{n} (x-\rho_i^{(n)}),    \end{aligned}$$ and for $T_n(x)$, the weights $a_i$ in (\[eq:ais\]) are all equal to ${2}/{n}$ when $w(x) =\frac{2}{\pi \sqrt{1-x^2}}.$

Notations {#sec:notation}
---------

Throughout this paper, we use lowercase bold letters to denote vectors and uppercase bold letters to denote matrices. In addition, for any positive integers $k,n$, and given a set of orthogonal polynomials $q_0,q_1, \cdots, q_{k-1}$ on the interval $[a,b]$, let $\mathbf{x}=(x_1, \cdots, x_n)$ be a vector with entries in $[a,b]$, we define the $k\times n$ matrix $\mathbf{Q}^{(k,n)}(\mathbf{x})$ as: $$\begin{aligned}
  \mathbf{Q}^{(k,n)}(\mathbf{x})= \left(\hspace{-2mm}\begin{array}{ccc} q_0(x_1)&  \cdots & q_{0}(x_n)\\\vdots & \ddots &\vdots \\ 
   q_{k-1}(x_1)&\cdots& q_{k-1}(x_n)\end{array}\hspace{-2mm}\right).\end{aligned}$$ For any subset $\mathcal{S}=\{s_1, \cdots, s_r\}\subset [n]$, we denote by $\mathbf{Q}^{(k,n)}_\mathcal{S}(\mathbf{x})$ the sub-matrix of $\mathbf{Q}^{(k,n)}(\mathbf{x})$ formed by concatenating columns with indices in $\mathcal{S}$, i.e., $$\begin{aligned}
  \mathbf{Q}_\mathcal{S}^{(k,n)}(\mathbf{x})= \left(\hspace{-2mm}\begin{array}{ccc} q_0(x_{s_1})&  \cdots & q_{0}(x_{s_r})\\\vdots & \ddots &\vdots \\ 
   q_{k-1}(x_{s_1})&\cdots& q_{k-1}(x_{s_r})\end{array}\hspace{-2mm}\right).\end{aligned}$$

For the special case where the orthogonal polynomials are the Chebyshev polynomials of the first kind $T_0, T_1, \cdots, T_{k-1}$, we define the $k\times n$ matrix $\mathbf{G}^{(k,n)}(\mathbf{x})$ as: $$\begin{aligned}
\label{eq:GTT}
  \mathbf{G}^{(k,n)}(\mathbf{x})= \left(\hspace{-2mm}\begin{array}{ccc} T_0(x_1)&  \cdots & T_{0}(x_n) \\ 
  \vdots & \ddots &\vdots \\ 
   T_{k-1}(x_1)&\cdots& T_{k-1}(x_n)\end{array}\hspace{-2mm}\right),\end{aligned}$$ we denote by $\mathbf{G}^{(k,n)}_\mathcal{S}(\mathbf{x})$ the sub-matrix of $\mathbf{G}^{(k,n)}(\mathbf{x})$ formed by concatenating columns with indices in $\mathcal{S}$, i.e., $$\begin{aligned}
\mathbf{G}_\mathcal{S}^{(k,n)}(\mathbf{x})= \left(\hspace{-2mm}\begin{array}{ccc}
  T_0(x_{s_1})&  \cdots & T_{0}(x_{s_r})\\\vdots & \ddots &\vdots \\ 
   T_{k-1}(x_{s_1})&\cdots& T_{k-1}(x_{s_r})\end{array}\hspace{-2mm}\right).\end{aligned}$$

Also, for the case where the orthogonal polynomials are the “orthonormal” Chebyshev polynomials $\frac{1}{\sqrt{2}}T_0, T_1, \cdots, T_{k-1}$, we define the $k\times n$ matrix $\tilde{\mathbf{G}}^{(k,n)}(\mathbf{x})$ as: $$\begin{aligned}
\label{eq:GT}
  \tilde{\mathbf{G}}^{(k,n)}(\mathbf{x})= \left(\hspace{-2mm}\begin{array}{ccc} T_0(x_1)/\sqrt{2}&  \cdots & T_{0}(x_n)/\sqrt{2} \\ 
   T_1(x_1)&  \cdots & T_{1}(x_n)\\\vdots & \ddots &\vdots \\ 
   T_{k-1}(x_1)&\cdots& T_{k-1}(x_n)\end{array}\hspace{-2mm}\right),\end{aligned}$$ and we denote by $\tilde{\mathbf{G}}^{(k,n)}_\mathcal{S}(\mathbf{x})$ the sub-matrix of $\tilde{\mathbf{G}}^{(k,n)}(\mathbf{x})$ formed by concatenating columns with indices in $\mathcal{S}$, i.e., $$\begin{aligned}
  \tilde{\mathbf{G}}_\mathcal{S}^{(k,n)}(\mathbf{x})= \left(\hspace{-2mm}\begin{array}{ccc} T_0(x_{s_1})/\sqrt{2}&  \cdots & T_{0}(x_{s_r})/\sqrt{2}\\
  T_1(x_{s_1})&  \cdots & T_{1}(x_{s_r})\\\vdots & \ddots &\vdots \\ 
   T_{k-1}(x_{s_1})&\cdots& T_{k-1}(x_{s_r})\end{array}\hspace{-2mm}\right).\end{aligned}$$ Wherever there is no ambiguity on $\mathbf{x}$, it may be dropped from the notation.

In the next section, we show that orthonormal polynomials can be used for designing codes for the distributed large scale matrix multiplication problem.

 Orthonormal Polynomials based Codes for Distributed Matrix Multiplication {#sec:ortho}
==========================================================================

In this section, we present a new orthonormal polynomials based class of codes for matrix-multiplication . These codes achieve the same recovery threshold as MatDot Codes, and have similar computational complexity as MatDot. The main advantage of the proposed codes is that they avoid dealing with the ill-conditioned monomial basis used in previous work (e.g., in [@allerton17; @polynomialcodes; @genPolyDot; @entPoly]). In Section \[sec:Chebyshev\], Codes will be specialized and demonstrated to have higher numerical stability as compared with state of the art. We begin with a formal problem formulation in Section \[sec:sysprob\], and describe our codes in Section \[sec:orthosub\].

System Model and Problem Formulation {#sec:sysprob}
------------------------------------

### System Model {#sec:sysmod}

![The distributed system framework[]{data-label="fig:Model"}](Model.pdf)

We consider the distributed framework depicted in Fig. \[fig:Model\] that consists of a master node, $P$ worker nodes, and a fusion node where the only communication allowed is from the master node to the different worker nodes and from the worker nodes to the fusion node. It can happen that the fusion node and the master node be represented by the same node. In this case, the only communication allowed is the communication between the master node and every worker node.

### Problem Formulation {#sec:probfor}

The master node possesses two real-valued input matrices $\mathbf{A}$, $\mathbf{B}$ with dimensions $N_1\times N_2$, $N_2\times N_3$, respectively. Every worker node receives from the master node an encoded matrix of $\mathbf{A}$ of dimension $N_1 \times N_2/m$ and an encoded matrix of $\mathbf{B}$ of dimension $N_2/m \times N_3$, and performs matrix multiplication of these two received inputs. Upon performing the matrix multiplication, each worker node sends the result to the fusion node. The fusion node needs to recover the matrix multiplication $\mathbf{A}\mathbf{B}$ once it receives the results of any $K$ worker nodes, where $K \leq P$. In this case, $K$ is denoted by *the recovery threshold* of the distributed computing scheme.

 Code Construction {#sec:orthosub}
------------------

Our result regarding the existence of achievable codes solving the distributed matrix multiplication problem using orthonormal polynomials is stated in the following theorem.

\[thm:orthpoly\] For the matrix multiplication problem described in Section \[sec:probfor\] computed on the system defined in Section \[sec:sysmod\], a recovery threshold of $2m-1$ is achievable using any set of orthonormal polynomials $\{q_i\}_{i\geq 0}$ relative to some weight polynomial $w$ and defined on a range $[a,b]$.

Before proving this theorem, we first present , a code construction that achieves the recovery threshold of $2m-1$ given any set $\{q_i\}_{i\geq 0}$ of orthonormal polynomials relative to a weight polynomial $w(x)$ and defined on a range $[a,b]$. In our code construction, we assume that matrix $\mathbf{A}$ is split vertically into $m$ equal sub-matrices, of dimension $N_1\times N_2/m$ each, and matrix $\mathbf{B}$ is split horizontally into $m$ equal sub-matrices, of dimension $N_2/m \times N_3 $ each, as follows: $$\label{eq:spltA}
\mathbf{A} = \left(\mathbf{A}_0 \ \mathbf{A}_1 \ \ldots \ \mathbf{A}_{m-1}\right),\;\;\; \mathbf{B}=\left(\begin{array}c \mathbf{B}_0\\\mathbf{B}_1\\\vdots \\\mathbf{B}_{m-1}\end{array}\right),$$ we also define a set of $P$ distinct real numbers $x_1, \cdots, x_P$ in the range $[a,b]$, and define two encoding polynomials $p_\mathbf{A}(x)=\sum_{i=0}^{m-1} \mathbf{A}_i q_{i}(x)$ and $p_\mathbf{B}(x)=\sum_{i=0}^{m-1} \mathbf{B}_i q_{i}(x),$ and let $p_\mathbf{C}(x)=p_{\mathbf{A}}(x)p_\mathbf{B}(x)$.

In the following, we briefly describe the construction. First, for every $r \in[P]$, the master node sends to the $r$-th worker node evaluations of $p_\mathbf{A}(x),p_\mathbf{B}(x)$ at $x=x_r$, that is, it sends $p_\mathbf{A}(x_r)$ and $p_\mathbf{B}(x_r)$ to the $r$-th worker node. Next, for every $r \in [P]$, the $r$-th worker node computes the matrix product $p_{\mathbf{C}}(x_r)=p_\mathbf{A}(x_r) p_\mathbf{B}(x_r)$ and sends the result to the fusion node. Once the fusion node receives the output of any $2m-1$ worker nodes, it interpolates the polynomial $p_{\mathbf{C}}(x)=p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$, and evaluates $p_{\mathbf{C}}(x)$ at $\eta_1, \cdots, \eta_{m}$, where $\eta_1, \cdots, \eta_m$ are the roots of $q_m$. Then, it performs the summation $\sum_{r=1}^{m} a_r\hspace{0.5mm} p_{\mathbf{C}}(\eta_r)$, where $a_1, \cdots, a_m$ are as in (\[eq:ais\]).

We formally present in Construction \[con:orthpoly\]. Construction \[con:orthpoly\] uses the following notation. The output of the algorithm is the $N_1\times N_3$ matrix $\hat{\mathbf{C}}.$ The $(i,j)$-th entries of the matrix polynomial $p_{\mathbf{C}}(x)$ and the matrix $\hat{\mathbf{C}}$ are respectively denoted as $p^{(i,j)}_{\mathbf{C}}(x)$ and $\hat{C}(i,j).$ The reader may also recall the definition of matrices $\mathbf{Q}^{(2m-1,P)}(\mathbf{x})$ and $\mathbf{Q}_{\mathcal{R}}^{(2m-1,P)}(\mathbf{x}),$ for any subset $\mathcal{R}=\{{r_1}, \cdots, {r_{2m-1}}\}  \subset [P]$. $\boldsymbol{\eta}=(\eta_1, \cdots, \eta_m)$ is the vector of the roots of $q_m$. Based on Construction \[con:orthpoly\], we state the following claim.

$r \gets 1$ $p_\mathbf{A}(x_r) \gets \sum_{i=0}^{m-1} \mathbf{A}_i q_{i}(x_r)$ $p_\mathbf{B}(x_r) \gets \sum_{i=0}^{m-1} \mathbf{B}_i q_{i}(x_r)$ **send** $p_{\mathbf{A}}(x_r), p_{\mathbf{B}}(x_r)$ **to worker node** $r$ $r \gets r+1$ $p_\mathbf{C}(x_r) \gets p_{\mathbf{A}}(x_r) p_{\mathbf{B}}(x_r)$ **send** $p_{\mathbf{C}}(x_r)$ **to the fusion node** $\mathbf{Q}_{\operatorname{inv}}\gets\left(\mathbf{Q}_{\mathcal{R}}^{(2m-1,P)}\right)^{-1}$ $(c^{(i,j)}_0, \cdots, c^{(i,j)}_{2m-2}) \gets (p_{\mathbf{C}}^{(i,j)}(x_{r_1}), \cdots,p_{\mathbf{C}}^{(i,j)}(x_{r_{2m-1}}))\mathbf{Q}_{\operatorname{inv}}$ $(p_\mathbf{C}^{(i,j)}(\eta_1), \cdots, p_{\mathbf{C}}^{(i,j)}(\eta_m))\gets (c^{(i,j)}_0, \cdots, c^{(i,j)}_{2m-2}) \mathbf{Q}^{(2m-1,m)}(\boldsymbol{\eta})$ $\hat{C}(i,j)\gets (p_\mathbf{C}^{(i,j)}(\eta_1), \cdots, p_{\mathbf{C}}^{(i,j)}(\eta_m)) (a_1, \cdots, a_m)^T $ **return** $\hat{\mathbf{C}}$

\[cl:orthAB\] $\mathbf{A}\mathbf{B}=\sum_{r=1}^{m} a_r\hspace{0.5mm} p_{\mathbf{C}}(\eta_r).$

The proof of Claim \[cl:orthAB\] is provided in Appendix \[app:ortho\].

Now, we can prove Theorem \[thm:orthpoly\].

In order to prove the theorem, it suffices to show that Construction \[con:orthpoly\] is a valid construction with a recovery threshold of $2m-1$. Therefore, in the following, we prove that Construction \[con:orthpoly\] can recover $\mathbf{A}\mathbf{B}$ after the fusion node receives the output of at most $2m-1$ worker nodes. Assume that the fusion node has already received the results of any $2m-1$ worker nodes. Now, because the polynomial $p_{\mathbf{C}}(x)$ has degree $2m-2$, the evaluations of $p_{\mathbf{C}}(x)$ at any $2m-1$ distinct points is sufficient to interpolate the polynomial, and since $x_1, \cdots, x_P$ are distinct, the fusion node can interpolate $p_{\mathbf{C}}(x)$ once it receives the output of any $2m-1$ worker nodes. Afterwards, given that $\mathbf{A}\mathbf{B}=\sum_{r=1}^{m} a_r\hspace{0.5mm} p_{\mathbf{C}}(\eta_r)$ (Claim \[cl:orthAB\]), the fusion node can evaluate $p_{\mathbf{C}}(\eta_1), \cdots, p_{\mathbf{C}}(\eta_m)$ and perform the scaled summation $\sum_{r=1}^{m} a_r\hspace{0.5mm} p_{\mathbf{C}}(\eta_r)$ to recover $\mathbf{A}\mathbf{B}$.

In Construction \[con:orthpoly\], setting $x_1, \cdots, x_m$ to be the roots of $q_m$ leads to a faster decoding for the scenarios in which the first $m$ worker nodes send their results but only less than $2m-1$ workers succeed to send their outputs. For such scenarios, we have $\sum_{r=1}^{m} a_r\hspace{0.5mm} p_{\mathbf{C}}(x_r)=$ $\sum_{r=1}^{m} a_r\hspace{0.5mm} p_{\mathbf{C}}(\eta_r)=\mathbf{A}\mathbf{B}$, where the last equality follows from Claim \[cl:orthAB\].

Next, we study the computational and communication costs of .

### Complexity Analyses of  {#sec:con1complexity}

**Encoding Complexity**: Encoding for each worker requires performing two additions, each adding $m$ scaled matrices of size $N_1 N_2/m$ and $N_2 N_3 /m$, for an overall encoding complexity for each worker of $O( N_1 N_2+N_2N_3)$. Therefore, the overall computational complexity of encoding for $P$ workers is $O(N_1N_2P+N_2N_3P)$.

**Computational Cost per Worker**: Each worker multiplies two matrices of dimensions $N_1\times N_2/m$ and $N_2/m\times N_3$, requiring $O(N_1N_2N_3/m)$ operations.

**Decoding Complexity**: Since $p_\mathbf{C}(x)$ has degree $2m-2$, the interpolation of $p_\mathbf{C}(x)$ requires the inversion of a $2m-1 \times 2m-1$ matrix, with complexity $O(m^3)$, and performing $N_1N_3$ matrix-vector multiplications, each of them is between the inverted matrix and a column vector of length $2m-1$ of the received evaluations of the matrix polynomial $p_\mathbf{C}(x)$ at some position $(i,j) \in [N_1] \times [N_3]$, with complexity $O(N_1N_3m^2)$. Next, the evaluation of the polynomial $p_\mathbf{C}(x)$ at $\eta_1, \cdots, \eta_m$ requires a complexity of $O(N_1N_3m^2)$. Finally, performing the summation $\sum_{r=1}^m a_r p_\mathbf{C}(\eta_r)$ requires a complexity of $O(N_1N_3m)$. Thus, assuming that $m \ll N_1,N_3$, the overall decoding complexity is $O(m^3+2N_1N_3m^2+N_1N_3m)=O(N_1N_3m^2)$.

**Communication Cost**: The master node sends $O( N_1N_2P/m+N_2N_3P/m)$ symbols, and the fusion node receives $O(N_1 N_3 m)$ symbols from the successful worker nodes.

With the reasonable assumption that the dimensions of the input matrices $\mathbf{A},\mathbf{B}$ are large enough such that $N_1,N_2,N_3 \gg m,P$, we can conclude that the encoding and decoding costs at the master and fusion nodes, respectively, are negligible compared to the computation cost at each worker node.

Numerically Stable Codes for Matrix Multiplication via  {#sec:Chebyshev}
=======================================================

In this section, we specialize Codes by restricting the orthonormal polynomials to be Chebyshev polynomials of the first kind $\{T_i\}_{i\geq 0}$ with the evaluation points chosen to be the $P$-dimensional Chebyshev grid, i.e., $x_i=\rho_i^{(P)}, i\in [P]$. Our , described in Construction \[con:chebpoly\] in Section \[subsec:Chebyshev\], develops a decoding that involves inversion of a $2m-1 \times 2m-1$ sub-matrix of a $2m-1 \times P$ Chebyshev-Vandermonde matrix. One of the main technical results of this section (and paper), presented in Theorem \[thm:bound\] in Section \[subsec:condbound\], is an upper bound to the worst case condition number over all possible $2m-1 \times 2m-1$ sub-matrices of the $2m-1 \times P$ Chebeshev-Vandermonde matrix for the case where the distinct evaluation points $x_1, \cdots, x_P$ are chosen as the Chebyshev points of degree $P$, i.e., $x_i=\rho_i^{(P)}, i\in [P]$. In fact, the derived bound shows that the worst case condition number grows at most polynomially in $P$ at a fixed number of straggler/parity worker nodes. This is in contrast with the monomial basis codes where the condition number grows exponentially in $P$, even when there is no redundancy [@gautschi1987lower; @gautschi1990stable; @gautschi1974norm; @reichel1991chebyshev]. We show through numerical experiments in Section \[sec:numerical\] that our proposed codes provide significantly lower numerical errors as compared to MatDot Codes in [@allerton17].

Chebyshev Polynomials based Code Construction {#subsec:Chebyshev}
---------------------------------------------

Recalling from Example \[ex:cheb\] that $\frac{1}{\sqrt{2}} T_0, T_1, T_2, \cdots$ form an orthonormal polynomial set relative to the weight function $w(x)= \frac{2}{\pi \sqrt{1-x^2}}$, in Construction \[con:chebpoly\], we explain the application of Chebyshev polynomials of the first kind to Construction \[con:orthpoly\]. Note that, in Construction \[con:chebpoly\], we assume that the input matrices $\mathbf{A}$ and $\mathbf{B}$ are also split as in (\[eq:spltA\]), and let $x_1, x_2, \ldots, x_P$ be distinct real numbers in the range $[-1,1]$, and define the encoding functions $p_\mathbf{A}(x), p_{\mathbf{B}}(x)$ as $p_\mathbf{A}(x)=\frac{1}{\sqrt{2}}\mathbf{A}_0 T_0(x)+\sum_{i=1}^{m-1} \mathbf{A}_i T_{i}(x)$ and $p_\mathbf{B}(x)=\frac{1}{\sqrt{2}}\mathbf{B}_0 T_0(x)+\sum_{i=1}^{m-1} \mathbf{B}_i T_{i}(x),$ and let $p_\mathbf{C}(x)=p_\mathbf{A}(x)p_\mathbf{B}(x)$.

The idea of our Chebyshev polynomials based code is as follows: First, for every $r \in [P]$, the master node sends to the $r$-th worker node $p_\mathbf{A}(\rho^{(P)}_r)$ and $p_\mathbf{B}(\rho^{(P)}_r)$. Next, for every $r \in [P]$, the $r$-th worker node computes the matrix product $p_\mathbf{C}(\rho^{(P)}_r)=p_\mathbf{A}(\rho^{(P)}_r)p_\mathbf{B}(\rho^{(P)}_r)$ and sends the result to the fusion node. Once the fusion node receives the output of any $2m-1$ worker nodes, it interpolates $p_{\mathbf{C}}(x)$. Then, it evaluates $p_{\mathbf{C}}(x)$ at $\rho^{(m)}_1, \cdots, \rho^{(m)}_m,$ and computes $\sum_{i=1}^{m} a_i\hspace{0.5mm} p_{\mathbf{C}}(\rho^{(m)}_i)$,

A formal description of our Chebyshev polynomials based code is provided in Construction \[con:chebpoly\]. Finally, we assume that our construction returns an $N_1\times N_3$ matrix $\hat{\mathbf{C}}$ representing the result of the product $\mathbf{A}\mathbf{B}$, where the $(i,j)$-th entry of $\hat{\mathbf{C}}$ is $\hat{C}(i,j)$.

$r \gets 1$ $p_\mathbf{A}(\rho^{(P)}_r) \gets \frac{1}{\sqrt{2}}\mathbf{A}_0+ \sum_{i=1}^{m-1} \mathbf{A}_i T_{i}(\rho^{(P)}_r)$ $p_\mathbf{B}(\rho^{(P)}_r) \gets \frac{1}{\sqrt{2}}\mathbf{B}_0+\sum_{i=1}^{m-1} \mathbf{B}_i T_{i}(\rho^{(P)}_r)$ **send** $p_{\mathbf{A}}(\rho^{(P)}_r), p_{\mathbf{B}}(\rho^{(P)}_r)$ **to worker node** $r$ $r \gets r+1$ $p_\mathbf{C}(\rho^{(P)}_r) \gets p_{\mathbf{A}}(\rho^{(P)}_r) p_{\mathbf{B}}(\rho^{(P)}_r)$ **send** $p_{\mathbf{C}}(\rho^{(P)}_r)$ **to the fusion node** $\mathbf{G}_{\operatorname{inv}}\gets\left(\tilde{\mathbf{G}}_{\mathcal{R}}^{(2m-1,P)}\right)^{-1}$ $(c^{(i,j)}_0, \cdots, c^{(i,j)}_{2m-2}) \gets (p_{\mathbf{C}}^{(i,j)}(\rho^{(P)}_{r_1}), \cdots,p_{\mathbf{C}}^{(i,j)}(\rho^{(P)}_{r_{2m-1}}))\mathbf{G}_{\operatorname{inv}}$ $(p_\mathbf{C}^{(i,j)}(\rho_1^{(m)}),\hspace{-1mm}\cdots\hspace{-1mm}, p_{\mathbf{C}}^{(i,j)}(\rho_m^{(m)}))\gets (c^{(i,j)}_0, \cdots, c^{(i,j)}_{2m-2}) \tilde{\mathbf{G}}^{(2m-1,m)}(\boldsymbol{\rho}^{(m)})$ $\hat{C}(i,j)\gets \frac{2}{m} (p_\mathbf{C}^{(i,j)}(\rho^{(m)}_1), \cdots, p_{\mathbf{C}}^{(i,j)}(\rho^{(m)}_m)) (1, \cdots, 1)^T $ **return** $\hat{\mathbf{C}}$

![image](Pall-eps-converted-to.pdf)

------------------------------------------------------------------------

### Complexity Analyses:

The different encoding complexity, computational complexity per worker, decoding complexity and communication cost for stated in Section \[sec:con1complexity\].

Evaluation Points and Condition Number Bound {#subsec:condbound}
--------------------------------------------

\[thm:bound\] For any $s \in [n-1]$, $$\kappa^{max}_F(\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})) = O\left( (n-s)\sqrt{ns(n-s)}\left({2}n^2\right)^{s-1}\right),$$ where $\kappa^{max}_F$ denotes the worst case condition number over all possible $n-s \times n-s$ sub-matrices of $\mathbf{G}^{(n-s,n)}(\mathbf{x})$ with respect to the Frobenius norm, $\boldsymbol{\rho}^{(n)}=(\rho_1^{(n)},\rho_2^{(n)},\ldots,\rho_n^{(n)})$ are the roots of the Chebyshev polynomial $T_n$, i.e., $\rho^{(n)}_i =  \cos\left(\frac{2i-1}{2n} \pi\right), i\in [n]$.

Since $||.||_2 \leq ||.||_F,$ the above bound applies to the standard $\ell_2$ matrix norm as well. The proof uses techniques from numerical methods, and is provided in Appendix \[app:proofgen\].

Although the bound in Theorem \[thm:bound\] is derived for $\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})$, the theorem also applies for $\tilde{\mathbf{G}}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})$. This is because it can be shown using simple matrix operations that for any $\tilde{\mathbf{G}}^{(n-s,n)}_{\mathcal{R}}$, for a subset $\mathcal{R}\subset [n]$ such that $|\mathcal{R}|=n-s$, $\kappa_F(\tilde{\mathbf{G}}^{(n-s,n)}_{\mathcal{R}})<\sqrt{2}\hspace{1mm}\kappa_F(\mathbf{G}^{(n-s,n)}_{\mathcal{R}}).$

Numerical Results {#sec:numerical}
-----------------

![image](fixed_parity-eps-converted-to.pdf)

------------------------------------------------------------------------

$$\begin{aligned}
  \mathbf{M}= \left(\hspace{-2mm}\begin{array}{cccc} 1& 1 & \cdots & 1\\ x_1 & x_2 & \cdots & x_P\\\vdots & \vdots &\ddots &\vdots \\ 
  x_1^{2m-2}& x_2^{2m-2}&\cdots& x_P^{2m-2}\end{array}\hspace{-2mm}\right).\end{aligned}$$

As a reflection to the significant higher stability of compared to MatDot Codes, Fig. \[fig:fixparerr\] shows that provides much more accurate outputs compared to MatDot Codes. For the experiments whose results are shown in Fig. \[fig:fixparerr\], the entries of the input matrices $\mathbf{A},\mathbf{B}$ are chosen independently according to the standard Gaussian distribution $\mathcal{N}(0,1)$. In addition, for any two input matrices $\mathbf{A}, \mathbf{B}$, let $\hat{\mathbf{C}}$ be the output of the distributed system (which is not necessarily equal to the correct answer $\mathbf{A}\mathbf{B}$), we define the relative error between $\mathbf{A}\mathbf{B}$ and $\hat{\mathbf{C}}$ to be $$\begin{aligned}
    E_r(\mathbf{A}\mathbf{B},\hat{\mathbf{C}})=\frac{||\mathbf{A}\mathbf{B}-\hat{\mathbf{C}}||_F}{||\mathbf{A}\mathbf{B}||_F}.\end{aligned}$$

![image](fixed_parity_error-eps-converted-to.pdf)

------------------------------------------------------------------------

Fig. \[fig:fixparerr\] shows how the maximum relative error (the worst case relative error given a fixed number of parity workers $s$ among all the $P-s$ successful nodes scenarios) grows with the size of the distributed system. In Fig. \[fig:fixparerr\], we plot the average result of five different realizations of the system at each system size $P$. The figure shows that MatDot Codes crushes after the size of the system exceeds $50$ workers, providing a relative error of around $10^5$. On the other hand, can support systems with sizes up to $150$ worker nodes only allowing for a relative error $< 10^{-5}$. It is also worth mentioning that in our experiments, we use the MATLAB command $inv()$ [@matlab] for matrix inversion. We have also tried matrices inversion through the Bjork-Pereyra algorithm [@bjorck1970solution], however, its results were much less accurate than $inv()$, especially for large systems with a number of worker nodes $> 50$.

A main challenge in this work is that we assume operations over the real field. For finite fields, one can always perform arithmetic operations with no errors. Although this fact may motivate a simple solution to the numerical stability of real-valued computations by rounding the computation’s inputs to a finite field’s elements and performing computations over this finite field, such solution has limited applicability, especially for inputs with wide range, due to the following reason. Since performing arithmetic operations over a finite field $\mathbb{F}_{2^n}$ requires representing each element of $\mathbb{F}_{2^n}$ as an element in $\mathbb{F}_{2}^{n}$ through a bit representation, this solution is applicable in machines with fixed point operations and word sizes of at least $n$. However, the solution is not applicable in machines with floating point operations since in floating point representation not all the intermediate values between the minimum and the maximum representable values can be represented, this is a drawback of the floating point representation over the fixed point representation, though floating point representation can represent a wider range of values than fixed point representation for the same word size.

 {#sec:chebpolycod}

While MatDot Codes [@allerton17] have an optimal recovery threshold of $2m-1$, they have relatively higher computation cost per worker and worker node to fusion node communication cost as compared to Polynomial Codes [@polynomialcodes]. In this section, motivated by the condition number bound in Theorem \[thm:bound\], we use the idea of using Chebyshev polynomials to provide a numerically stable code construction for matrix multiplication that has the same low communication/computation costs as Polynomial Codes, as well as the same recovery threshold. However, as will be shown in this section, our proposed codes, , provides lower numerical errors than Polynomial Codes. In this section, we follow the same system model as in Section \[sec:sysmod\], and solve the problem statement formulated in Section \[sec:chebprobfor\]. We provide a motivating example in Section \[sec:chebpolyex\], then we provide the general code construction in Section \[sec:chebpolygen\]. Finally, in Section \[sec:chebpolyexp\], we show experimentally that achieve lower numerical errors as compared to Polynomial Codes.

Problem Formulation {#sec:chebprobfor}
-------------------

The master node possesses two real-valued input matrices $\mathbf{A}$, $\mathbf{B}$ with dimensions $N_1\times N_2$, $N_2\times N_3$, respectively. Every worker node receives from the master node an encoded matrix of $\mathbf{A}$ of dimension $N_1/m \times N_2$ and an encoded matrix of $\mathbf{B}$ of dimension $N_2 \times N_3/n$, and performs matrix multiplication of these two received inputs. Upon performing the matrix multiplication, each worker node sends the result to the fusion node. The fusion node needs to recover the matrix multiplication $\mathbf{A}\mathbf{B}$ once it receives the results of any $mn$ worker nodes.

Example $(m=n=3)$ {#sec:chebpolyex}
-----------------

Consider computing the matrix multiplication $\mathbf{A}\mathbf{B}$, for some two real matrices $\mathbf{A}, \mathbf{B}$ of dimensions $N_1\times N_2$ and $N_2 \times N_3$, respectively, over a distributed system of $P \geq 9$ workers such that:

1.  Each worker receives an encoded matrix of $\mathbf{A}$ of dimension $N_1 /3\times N_2$, and an encoded matrix of $\mathbf{B}$ of dimension $N_2 \times N_3/3$.

2.  The product $\mathbf{A}\mathbf{B}$ can be recovered by the fusion node given the results of any $9$ worker nodes.

A solution can be as follows: First, matrices $\mathbf{A}, \mathbf{B}$ can be partitioned as $$\begin{aligned}
\label{eq:polsplt}
 \mathbf{A}&=\left(\begin{array}{ccc}\mathbf{A}_{0}\\\mathbf{A}_{1}\\\mathbf{A}_{2}\end{array}\right),~~
 \mathbf{B}=\left(\begin{array}{ccc}\mathbf{B}_{0}~~  \mathbf{B}_{1}~~\mathbf{B}_{2}\end{array}\right), \end{aligned}$$ where, for any $i\in\{0,1,2\}$, $\mathbf{A}_{i}$ has dimension $N_1/3 \times N_2$, and $\mathbf{B}_{i}$ has dimension $N_2 \times N_3/3$. Next, let $$\begin{aligned}
p_\mathbf{A}(x)=\mathbf{A}_{0}T_{0}(x)+\mathbf{A}_{1}T_{1}(x)+\mathbf{A}_{2}T_{2}(x),\notag\\
p_\mathbf{B}(x)=\mathbf{B}_{0}T_{0}(x)+\mathbf{B}_{1} T_{3}(x)+\mathbf{B}_{2}T_{6}(x).\notag\end{aligned}$$ Now, $p_\mathbf{A}(x)p_\mathbf{B}(x)$ can be written as $$\begin{aligned}
p_\mathbf{A}(x)p_\mathbf{B}(x)&=\big(\mathbf{A}_{0}T_{0}(x)+\mathbf{A}_{1}T_{1}(x)+\mathbf{A}_{2}T_{2}(x)\big)\big(\mathbf{B}_{0}T_{0}(x)+\mathbf{B}_{1} T_{3}(x)+\mathbf{B}_{2}T_{6}(x)\big)\notag\\
&= \mathbf{A}_0\mathbf{B}_0+\big(\mathbf{A}_1\mathbf{B}_0+\frac{1}{2}\mathbf{A}_2\mathbf{B}_1\big)T_1(x)+\big(\mathbf{A}_2\mathbf{B}_0+\frac{1}{2}\mathbf{A}_1\mathbf{B}_1\big)T_2(x)\notag\\
&+  \mathbf{A}_0\mathbf{B}_1 T_3(x)+\frac{1}{2}\big(\mathbf{A}_1\mathbf{B}_1+\mathbf{A}_2\mathbf{B}_2\big)T_4(x)+\frac{1}{2}\big(\mathbf{A}_1\mathbf{B}_2+\mathbf{A}_2\mathbf{B}_1\big)T_5(x)\notag\\
&+\mathbf{A}_0\mathbf{B}_2 T_6(x)+\frac{1}{2}\mathbf{A}_1\mathbf{B}_2T_7(x)+\frac{1}{2}\mathbf{A}_2\mathbf{B}_2T_8(x)\end{aligned}$$ Since $p_\mathbf{A}(x)p_\mathbf{B}(x)$ is a degree $8$ polynomial, once the fusion node receives the output of any $9$ workers, it can interpolate $p_\mathbf{A}(x)p_\mathbf{B}(x)$, i.e., obtain its matrix coefficients, let such matrix coefficients be $\mathbf{C}_{T_0}, \cdots, \mathbf{C}_{T_8}$. Specifically, for any $i\in \{0, \cdots, 8\}$, let $\mathbf{C}_{T_i}$ be the matrix coefficient of $T_i$ in $p_\mathbf{A}(x)p_\mathbf{B}(x)$. Now, recalling (\[eq:polsplt\]), the product $\mathbf{A}\mathbf{B}$ can be written as $$\begin{aligned}
\mathbf{A}\mathbf{B}=    \left(\begin{array}{ccc}\mathbf{A}_0\mathbf{B}_0&  \mathbf{A}_0\mathbf{B}_1& \mathbf{A}_0\mathbf{B}_2\\\mathbf{A}_1\mathbf{B}_0& \mathbf{A}_1\mathbf{B}_1 &\mathbf{A}_1\mathbf{B}_2\\\mathbf{A}_2\mathbf{B}_0&\mathbf{A}_2\mathbf{B}_1 &\mathbf{A}_2\mathbf{B}_2\end{array}\right).\end{aligned}$$ While the obtained set of matrix coefficients $\{\mathbf{C}_{T_i}: i \in \{0, \cdots, 8\}\}$ is not equal to $\{\mathbf{A}_i\mathbf{B}_j : i,j\in\{0,1,2\}\}$, $\mathbf{C}_{T_{i}}$’s are linear combinations of $\mathbf{A}_i\mathbf{B}_j$’s. Specifically, for any $\mathbf{C}_{T_i}$, $i \in \{0, \cdots, 8\}$, let $\mathbf{C}_{T_i}^{(k,l)}$ be its $(k,l)$-th entry, and, for any $i,j \in\{0,1,2\}$, let $(\mathbf{A}_i\mathbf{B}_j)^{(k,l)}$ be the $(k,l)$-th entry of the product $\mathbf{A}_i\mathbf{B}_j$, we can write $$\begin{aligned}
\label{eq:tcoef}
\left(\begin{array}{c}\mathbf{C}_{T_0}^{(k,l)}\\  \mathbf{C}_{T_1}^{(k,l)}\\ \mathbf{C}_{T_2}^{(k,l)}\\\mathbf{C}_{T_3}^{(k,l)}\\ \mathbf{C}_{T_4}^{(k,l)}\\\mathbf{C}_{T_5}^{(k,l)}\\\mathbf{C}_{T_6}^{(k,l)}\\\mathbf{C}_{T_7}^{(k,l)}\\\mathbf{C}_{T_8}^{(k,l)}\end{array}\right)= \left(\begin{array}{ccccccccc}
1&0&0&0&0&0&0&0&0\\0&1&0&0&0&1/2&0&0&0\\0&0&1&0&1/2&0&0&0&0
\\0&0&0&1&0&0&0&0&0\\0&0&0&0&1/2&0&0&0&1/2\\0&0&0&0&0&1/2&0&1/2&0\\0&0&0&0&0&0&1&0&0\\0&0&0&0&0&0&0&1/2&0\\0&0&0&0&0&0&0&0&1/2
\end{array}\right)\left(\begin{array}{c}(\mathbf{A}_0\mathbf{B}_0)^{(k,l)}\\  (\mathbf{A}_1\mathbf{B}_0)^{(k,l)}\\ (\mathbf{A}_2\mathbf{B}_0)^{(k,l)}\\(\mathbf{A}_0\mathbf{B}_1)^{(k,l)}\\ (\mathbf{A}_1\mathbf{B}_1)^{(k,l)}\\(\mathbf{A}_2\mathbf{B}_1)^{(k,l)}\\(\mathbf{A}_0\mathbf{B}_2)^{(k,l)}\\(\mathbf{A}_1\mathbf{B}_2)^{(k,l)}\\(\mathbf{A}_2\mathbf{B}_2)^{(k,l)}\end{array}\right),\end{aligned}$$ for any $(k,l) \in [N_1/3]\times[N_3/3]$. Thus, the products $\mathbf{A}_i\mathbf{B}_j,  i,j\in\{0,1,2\}$ can be obtained by computing $$\begin{aligned}
\label{eq:ABcoef}
\left(\begin{array}{c}(\mathbf{A}_0\mathbf{B}_0)^{(k,l)}\\  (\mathbf{A}_1\mathbf{B}_0)^{(k,l)}\\ (\mathbf{A}_2\mathbf{B}_0)^{(k,l)}\\(\mathbf{A}_0\mathbf{B}_1)^{(k,l)}\\ (\mathbf{A}_1\mathbf{B}_1)^{(k,l)}\\(\mathbf{A}_2\mathbf{B}_1)^{(k,l)}\\(\mathbf{A}_0\mathbf{B}_2)^{(k,l)}\\(\mathbf{A}_1\mathbf{B}_2)^{(k,l)}\\(\mathbf{A}_2\mathbf{B}_2)^{(k,l)}\end{array}\right)=\left(\begin{array}{ccccccccc}
1&0&0&0&0&0&0&0&0\\0&1&0&0&0&1/2&0&0&0\\0&0&1&0&1/2&0&0&0&0
\\0&0&0&1&0&0&0&0&0\\0&0&0&0&1/2&0&0&0&1/2\\0&0&0&0&0&1/2&0&1/2&0\\0&0&0&0&0&0&1&0&0\\0&0&0&0&0&0&0&1/2&0\\0&0&0&0&0&0&0&0&1/2
\end{array}\right)^{-1} \left(\begin{array}{c}\mathbf{C}_{T_0}^{(k,l)}\\  \mathbf{C}_{T_1}^{(k,l)}\\ \mathbf{C}_{T_2}^{(k,l)}\\\mathbf{C}_{T_3}^{(k,l)}\\ \mathbf{C}_{T_4}^{(k,l)}\\\mathbf{C}_{T_5}^{(k,l)}\\\mathbf{C}_{T_6}^{(k,l)}\\\mathbf{C}_{T_7}^{(k,l)}\\\mathbf{C}_{T_8}^{(k,l)}\end{array}\right),\end{aligned}$$ for all $(k,l) \in [N_1/3]\times[N_3/3]$. In the following, we provide the general code construction.

 Code Construction {#sec:chebpolygen}
------------------

We assume that matrix $\mathbf{A}$ is split horizontally into $m$ equal sub-matrices, of dimension $N_1/m\times N_2$ each, and matrix $\mathbf{B}$ is split vertically into $n$ equal sub-matrices, of dimension $N_2 \times N_3/n $ each, as follows: $$\label{eq:spltApol}
\mathbf{A}=\left(\begin{array}c \mathbf{A}_0\\\mathbf{A}_1\\\vdots \\\mathbf{A}_{m-1}\end{array}\right),\;\;\;\mathbf{B} = \left(\mathbf{B}_0 \ \mathbf{B}_1 \ \ldots \ \mathbf{B}_{n-1}\right),$$ and define two encoding polynomials $p_\mathbf{A}(x)=\sum_{i=0}^{m-1} \mathbf{A}_i T_{i}(x)$ and $p_\mathbf{B}(x)=\sum_{i=0}^{n-1} \mathbf{B}_i T_{im}(x),$ and let $p_\mathbf{C}(x)=p_{\mathbf{A}}(x)p_\mathbf{B}(x)$. We describe, next, the idea of the general code construction. First, for all $r \in [P]$, the master node sends to the $r$-th worker evaluations of $p_\mathbf{A}(x)$ and $p_\mathbf{B}(x)$ at $x=\rho^{(P)}_r$, that is, it sends $p_\mathbf{A}(\rho^{(P)}_r)$ and $p_\mathbf{B}(\rho^{(P)}_r)$ to the $r$-th worker. Next, for every $r \in [P]$, the $r$-th worker node computes the matrix product $p_\mathbf{C}(\rho^{(P)}_r)=p_{\mathbf{A}}(\rho^{(P)}_r)p_{\mathbf{B}}(\rho^{(P)}_r)$ and sends the result to the fusion node. Once the fusion node receives the output of any $mn$ worker nodes, it interpolates $p_\mathbf{C}(x)$. Next, the fusion node recovers the products $\mathbf{A}_i\mathbf{B}_j, i \in \{0, \cdots, m-1\}, j \in \{0, \cdots, n-1\}$, from the matrix coefficients of $p_\mathbf{C}(x)$ using a low complexity matrix-vector multiplication, specified later in Construction \[con:polycheb2\]. We formally present our Codes in Construction \[con:polycheb2\]. In the following, we explain the notation used in Construction \[con:polycheb2\]. The output of the algorithm is the $N_1\times N_3$ matrix $\hat{\mathbf{C}},$ where the $(k,l)$-th block of $\hat{\mathbf{C}}$ is the $N_1/m \times N_3/n$ matrix $\hat{\mathbf{C}}_{k,l}$, and the $(i,j)$-th entry of any matrix $\hat{\mathbf{C}}_{k,l}$ is $\hat{{c}}_{k,l}^{(i,j)}$. The $(i,j)$-th entry of the matrix polynomial $p_{\mathbf{C}}(x)$ is denoted as $p^{(i,j)}_{\mathbf{C}}(x)$, and Section \[sec:notation\] defines matrices $\mathbf{G}^{(mn,P)}(\boldsymbol{\rho}^{(P)})$ and $\mathbf{G}^{(mn,P)}_{\mathcal{R}}(\boldsymbol{\rho}^{(P)})$, for any subset $\mathcal{R}=\{{r_1}, \cdots, {r_{mn}}\}  \subset [P]$. In addition, $\mathbf{H}$ is an $mn \times mn$ matrix of the following form $\mathbf{H}=\left(\begin{array}{ccc}\mathbf{H}_{0}~~  \mathbf{H}_{1}~~\cdots~~\mathbf{H}_{n-1}\end{array}\right),$ where $\mathbf{H}_0$ is an $mn \times m$ matrix with ones on the main diagonal and zeros elsewhere, and for any $i\in \{1, \cdots, n-1\}$, $\mathbf{H}_i$ is an $mn\times m$ matrix of the following structure $$\mathbf{H}_i=\left(\begin{array}{ccccccccc}
0&0&0&\cdots&0\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&0&\cdots&0\\\vdots&\vdots&\vdots&\iddots&1/2\\0&0&0&\iddots&0\\0&0&1/2&\iddots&\vdots\\0&1/2&0&\cdots&0\\1&0&0&\cdots&0\\0&1/2&0&\cdots&0\\0&0&1/2&\ddots&\vdots\\0&0&0&\ddots&0\\\vdots&\vdots&\vdots&\ddots&1/2\\0&0&0&\cdots&0\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&0&\cdots&0
\end{array}\right),$$ where the value $1$ in the first column is at the $(im+1)$-th row of $\mathbf{H}_i$.

$r \gets 1$ $p_\mathbf{A}(\rho^{(P)}_r) \gets \sum_{i=0}^{m-1}\mathbf{A}_{i} T_{i}(\rho^{(P)}_r)$ $p_\mathbf{B}(\rho^{(P)}_r)\gets\sum_{i=0}^{n-1} \mathbf{B}_{i} T_{im}(\rho^{(P)}_r)$ **send** $p_{\mathbf{A}}(\rho^{(P)}_r), p_{\mathbf{B}}(\rho^{(P)}_r)$ **to worker node** $r$ $r \gets r+1$ $p_\mathbf{C}(\rho^{(P)}_r) \gets p_{\mathbf{A}}(\rho^{(P)}_r) p_{\mathbf{B}}(\rho^{(P)}_r)$ **send** $p_{\mathbf{C}}(\rho^{(P)}_r)$ **to the fusion node** $\mathbf{G}_{\operatorname{inv}}\gets \left(\mathbf{G}^{(mn,P)}_{\mathcal{R}}\right)^{-1}$ $(c^{(i,j)}_0, \cdots, c^{(i,j)}_{mn-1}) \gets (p_{\mathbf{C}}^{(i,j)}(\rho^{(P)}_{r_1}), \cdots,p_{\mathbf{C}}^{(i,j)}(\rho^{(P)}_{r_{mn}}))\mathbf{G}_{\operatorname{inv}}$ $(\hat{c}_{0,0}^{(i,j)} \cdots \hat{c}_{m-1,0}^{(i,j)} \cdots\cdots \hat{c}_{0,n-1}^{(i,j)} \cdots \hat{c}_{m-1,n-1}^{(i,j)}) \gets (c^{(i,j)}_0, \cdots, c^{(i,j)}_{mn-1}) (\mathbf{H}^{-1})^T$ **return** $\hat{\mathbf{C}}$

### Complexity Analyses of  {#sec:con4complexity}

**Encoding Complexity**: Encoding for each worker requires performing two additions, the first one adds $m$ scaled matrices of size $N_1 N_2/ m$ and the other adds $n$ scaled matrices of size $N_2 N_3 /n$, for an overall encoding complexity for each worker of $O(N_1 N_2+N_2N_3)$. Therefore, the overall computational complexity of encoding for $P$ workers is $O(N_1 N_2P+N_2N_3P)$.

**Computational Cost per Worker**: Each worker multiplies two matrices of dimensions $N_1  /m\times N_2$ and $N_2\times N_3/n$, requiring $O(N_1N_2N_3/mn)$ operations.

**Decoding Complexity**: Since $p_\mathbf{A}(x)p_\mathbf{B}(x)$ has degree $mn-1$, the interpolation of $p_\mathbf{C}(x)$ requires the inversion of a $mn \times mn$ matrix, with complexity $O(m^3n^3)$, and performing $N_1N_3/mn$ matrix-vector multiplications, each of them is between the inverted matrix and a column vector of length $mn$ of the received evaluations of the matrix polynomial $p_\mathbf{C}(x)$ at some position $(i,j) \in [N_1/m] \times [N_3/n]$, with complexity $O(N_1N_3m^2n^2/(mn))=O(N_1N_3mn)$. Thus, assuming that $mn \ll N_1,N_3$, the overall decoding complexity is $O(N_1N_3 m n)$.

**Communication Cost**: The master node sends $O( N_1N_2P/m+N_2N_3P/n)$ symbols, and the fusion node receives $O(N_1 N_3)$ symbols from the successful worker nodes.

With the reasonable assumption that the dimensions of the input matrices $\mathbf{A},\mathbf{B}$ are large enough such that $N_1,N_2,N_3 \gg m,n,P$, we can conclude that the encoding and decoding costs at the master and fusion nodes, respectively, are negligible compared to the computation cost at each worker node.

Numerical Results {#sec:chebpolyexp}
-----------------

In our experiments, the entries of the input matrices $\mathbf{A},\mathbf{B}$ are chosen independently according to the standard Gaussian distribution $\mathcal{N}(0,1)$. In addition, for any two input matrices $\mathbf{A}, \mathbf{B}$, let $\hat{\mathbf{C}}$ be the output of the distributed system, we define the relative error between $\mathbf{A}\mathbf{B}$ and $\hat{\mathbf{C}}$ to be $$\begin{aligned}
    E_r(\mathbf{A}\mathbf{B},\hat{\mathbf{C}})=\frac{||\mathbf{A}\mathbf{B}-\hat{\mathbf{C}}||_F}{||\mathbf{A}\mathbf{B}||_F}.\end{aligned}$$ Fig. \[fig:fixparerrpoly\] shows how the maximum relative error (the worst case relative error given a fixed number of parity workers $s$ among all the $P-s$ successful nodes scenarios) grows with the size of the distributed system for both Construction \[con:polycheb2\] and Polynomial Codes. In Fig. \[fig:fixparerrpoly\], we plot the average result of five different realizations of the system at each system size $P$. The figure shows that Polynomial Codes have unacceptable relative errors after the size of the system exceeds $50$ workers, providing a relative error of around $10^5$. On the other hand, can support systems with sizes up to $170$ worker nodes only allowing for a relative error $< 10^{-5}$.

![image](fixed_parity_error_polycod-eps-converted-to.pdf)

------------------------------------------------------------------------

 {#sec:polydot-chebyshev}

Although MatDot Codes [@allerton17] have a low recovery threshold of $2m-1$ as compared with Polynomial Codes [@polynomialcodes] which have a recovery threshold of $mn$, MatDot Codes’ worker to fusion nodes communication cost and computation cost per worker are higher than Polynomial Codes. Codes proposed in [@arxiv_allerton17; @genPolyDot; @entPoly] offer a trade-off between the communication/computation cost and the recovery threshold. However, all of these codes are based on the “ill-conditioned” monomial basis. In this section, we offer a numerically stable code construction, that offers a trade-off between communication/computation costs and recovery threshold. Our construction incurs a higher recovery threshold than the codes of [@genPolyDot; @entPoly] by a factor of at most $4$ for the same communication/computation cost. We provide in Section \[sec:probfor2\] the formal problem statement considered in this section. We describe an example of our construction in Section \[sec:polydot\_example\], provide the general code construction in Section \[sec:Chebyshev\_polydot\], and describe our numerical experiments in Section \[sec:numerical\_polydot\].

System Model and Problem Formulation {#sec:probfor2}
------------------------------------

We consider the same system model and problem formulation as in Section \[sec:sysprob\] with the following change: We assume that the master node is allowed to send an encoded ${1}/{m}$ fraction of matrix $\mathbf{A}$, and an encoded ${1}/{n}$ fraction of matrix $\mathbf{B}$, where $m$ and $n$ are not necessarily equal, and $\mathbf{A}$ and $\mathbf{B}$ are split as follows $$\begin{aligned}
\label{eq:spltA2}
 \mathbf{A}&=\left(\begin{array}{ccc}\mathbf{A}_{0,0}&  \cdots& \mathbf{A}_{0,m_2-1}\\\vdots& \ddots &\vdots\\\mathbf{A}_{m_1-1,0} &\cdots &\mathbf{A}_{m_1-1,m_2-1}\end{array}\right),\notag\\
 \mathbf{B}&=\left(\begin{array}{ccc}\mathbf{B}_{0,0}&  \cdots& \mathbf{B}_{0,m_3-1}\\\vdots& \ddots &\vdots\\\mathbf{B}_{m_2-1,0} &\cdots &\mathbf{B}_{m_2-1,m_3-1}\end{array}\right), \end{aligned}$$ where $m_1,m_2,m_3$ divide $N_1,N_2,N_3$, respectively, and $m=m_1m_2, n=m_2m_3$. In addition, we assume that each worker node receives a linear combination of sub-matrices $\mathbf{A}_{i,j}$, and another linear combination of sub-matrices $\mathbf{B}_{i,j}$.

Although, in this section, we offer a code construction with lower condition numbers than codes in [@genPolyDot; @entPoly], the recovery threshold of our codes are higher by a factor of at most $4$ than the codes of these references. Specifically, codes have a recovery threshold of $4m_1m_2m_3-2(m_1m_2+m_2m_3+m_3m_1)$ $+m_1+2m_2+m_3-1$ while both codes in [@genPolyDot; @entPoly] have a recovery threshold of $m_1m_2m_3+m_2-1$. This increased recovery threshold is due to the fact that Codes are based on Chebyshev polynomials which have the following property: For any $i,j \in \mathbb{N}$, $T_{i}(x)T_{j}(x)=1/2~(T_{i+j}(x)+T_{|i-j|}(x))$. This property allows for a higher number of *undesired* terms in the multiplication of the encoding polynomials $p_{\mathbf{A}}(x), p_{\mathbf{B}}(x)$. In order to avoid combining undesired and desired terms at the same degree, higher degree Chebyshev polynomials have to be used in $p_{\mathbf{B}}(x)$, yielding a higher recovery threshold. It is still an open question whether the recovery threshold in [@genPolyDot; @entPoly] can be achieved using orthonormal polynomials.

Example $(m_1=m_2=m_3=2)$ {#sec:polydot_example}
-------------------------

Consider computing the matrix multiplication $\mathbf{A}\mathbf{B}$, for some two real matrices $\mathbf{A}, \mathbf{B}$ of dimensions $N_1\times N_2$ and $N_2 \times N_3$, respectively, over a distributed system of $P \geq 15$ workers such that:

1.  Each worker receives an encoded matrix of $\mathbf{A}$ of dimension $N_1/2 \times N_2/2$, and an encoded matrix of $\mathbf{B}$ of dimension $N_2/2 \times N_3/2$.

2.  The product $\mathbf{A}\mathbf{B}$ can be recovered by the fusion node given the results of any $15$ worker nodes.

A solution can be as follows: First, matrices $\mathbf{A}, \mathbf{B}$ can be partitioned as $$\begin{aligned}
\label{eq:exsplt}
 \mathbf{A}&=\left(\begin{array}{ccc}\mathbf{A}_{0,0}&  \mathbf{A}_{0,1}\\\mathbf{A}_{1,0}  &\mathbf{A}_{1,1}\end{array}\right),
 \mathbf{B}=\left(\begin{array}{ccc}\mathbf{B}_{0,0}& \mathbf{B}_{0,1}\\\mathbf{B}_{1,0} &\mathbf{A}_{1,1}\end{array}\right), \end{aligned}$$ where, for $i,j \in\{0,1\}$, $\mathbf{A}_{i,j}$ has dimension $N_1/2 \times N_2/2$, and $\mathbf{B}_{i,j}$ has dimension $N_2/2 \times N_3/2$. Next, let $$\begin{aligned}
p_\mathbf{A}(x)=\mathbf{A}_{0,0}T_{1}(x)+\frac{1}{2}\mathbf{A}_{0,1}T_{0}(x)+\mathbf{A}_{1,0}T_{\alpha+1}(x)+\mathbf{A}_{1,1}T_{\alpha}(x),\notag\\
p_\mathbf{B}(x)=\frac{1}{2}\mathbf{B}_{0,0}T_{0}(x)+\mathbf{B}_{1,0} T_{1}(x)+\mathbf{B}_{0,1}T_{\beta}(x)+\mathbf{B}_{1,1}T_{\beta+1}(x),\notag\end{aligned}$$ where $\alpha, \beta$ to be specified next, and define $P$ distinct real numbers $x_1,x_2, \cdots, x_P$ in the range $[-1,1]$. For each worker node $r \in [P]$, the master node sends $p_{\mathbf{A}}(x_r)p_{\mathbf{B}}(x_r)$.

Now, in order to specify the best values for $\alpha,\beta$, we expand the polynomial $p_\mathbf{A}(x)p_\mathbf{B}(x)$ in the Chebyshev basis, and then point out some observations. $$\begin{aligned}
\label{eq:expandab}
    p_\mathbf{A}(x)p_\mathbf{B}(x)&=\frac{1}{4} \mathbf{A}_{0,1}\mathbf{B}_{0,0}+\frac{1}{2} \mathbf{A}_{0,0}\mathbf{B}_{0,0}T_{1}(x)+\frac{1}{2}\mathbf{A}_{0,1}\mathbf{B}_{1,0}T_{1}(x)+ \mathbf{A}_{0,0}\mathbf{B}_{1,0}T_1(x)T_{1}(x)\notag\\
    &+ \frac{1}{2} \mathbf{A}_{1,1}\mathbf{B}_{0,0}T_{\alpha}(x)+  \mathbf{A}_{1,1}\mathbf{B}_{1,0}T_{1}(x)T_{\alpha}(x)+ \frac{1}{2} \mathbf{A}_{1,0}\mathbf{B}_{0,0}T_{\alpha+1}(x)+  \mathbf{A}_{1,0}\mathbf{B}_{1,0}T_{1}(x)T_{\alpha+1}(x)\notag\\
    &+ \frac{1}{2} \mathbf{A}_{0,1}\mathbf{B}_{0,1}T_{\beta}(x)+  \mathbf{A}_{0,0}\mathbf{B}_{0,1}T_{1}(x)T_{\beta}(x)+ \frac{1}{2} \mathbf{A}_{0,1}\mathbf{B}_{1,1}T_{\beta+1}(x)+  \mathbf{A}_{0,0}\mathbf{B}_{1,1}T_{1}(x)T_{\beta+1}(x)\notag\\
     &+ \mathbf{A}_{0,1}\mathbf{B}_{1,1}T_{1}(x)T_{\beta+1}(x)+  \mathbf{A}_{1,0}\mathbf{B}_{0,1}T_{\alpha+1}(x)T_{\beta}(x)+\mathbf{A}_{1,1}\mathbf{B}_{1,1}T_{\alpha}(x)T_{\beta+1}(x) \notag\\
     &+\mathbf{A}_{1,0}\mathbf{B}_{1,1}T_{\alpha+1}(x)T_{\beta+1}(x).\end{aligned}$$ Using the property of the Chebyshev polynomials that for any $i,j \in \mathbb{N}$, $T_{i}(x)T_{j}(x)=1/2~(T_{i+j}(x)+T_{|i-j|}(x))$, (\[eq:expandab\]) can be rewritten as $$\begin{aligned}
\label{eq:expandab2}
    p_\mathbf{A}(x)p_\mathbf{B}(x)&=\frac{1}{4} \mathbf{A}_{0,1}\mathbf{B}_{0,0}+\frac{1}{2}\mathbf{A}_{0,0}\mathbf{B}_{1,0}+\frac{1}{2} \left(\mathbf{A}_{0,0}\mathbf{B}_{0,0}+\mathbf{A}_{0,1}\mathbf{B}_{1,0}\right)T_{1}(x)+\frac{1}{2} \mathbf{A}_{0,0}\mathbf{B}_{1,0}T_{2}(x)\notag\\
    &+\frac{1}{2}\mathbf{A}_{1,1}\mathbf{B}_{1,0}T_{\alpha-1}(x)+\frac{1}{2} \left(\mathbf{A}_{1,1}\mathbf{B}_{0,0}+\mathbf{A}_{1,0}\mathbf{B}_{1,0}\right) T_{\alpha}(x)+\frac{1}{2} \left(\mathbf{A}_{1,0}\mathbf{B}_{0,0}+\mathbf{A}_{1,1}\mathbf{B}_{1,0}\right) T_{\alpha+1}(x)\notag\\
    &+\frac{1}{2}\mathbf{A}_{1,0}\mathbf{B}_{1,0}T_{\alpha+2}(x)+\frac{1}{2}\mathbf{A}_{1,0}\mathbf{B}_{0,1}T_{\beta-\alpha-1}(x)+
    \frac{1}{2}\left(\mathbf{A}_{1,1}\mathbf{B}_{0,1}+\mathbf{A}_{1,0}\mathbf{B}_{1,1}\right)T_{\beta-\alpha}(x)\notag\\
    &+\frac{1}{2}\mathbf{A}_{1,1}\mathbf{B}_{1,1}T_{\beta-\alpha+1}(x)+\frac{1}{2}\mathbf{A}_{0,0}\mathbf{B}_{0,1}T_{\beta-1}(x)+\frac{1}{2}\left(\mathbf{A}_{0,1}\mathbf{B}_{0,1}+\mathbf{A}_{0,0}\mathbf{B}_{1,1}\right)T_{\beta}(x)\notag\\
    &+ \frac{1}{2}\left(\mathbf{A}_{0,0}\mathbf{B}_{0,1}+\mathbf{A}_{0,1}\mathbf{B}_{1,1}\right)T_{\beta+1}(x)+\frac{1}{2}\mathbf{A}_{0,0}\mathbf{B}_{1,1}T_{\beta+2}(x)+\frac{1}{2}\mathbf{A}_{1,1}\mathbf{B}_{0,1}T_{\beta+\alpha}(x)\notag\\
    &+ \frac{1}{2} \left(\mathbf{A}_{1,0}\mathbf{B}_{0,1}+\mathbf{A}_{1,1}\mathbf{B}_{1,1} \right)T_{\beta+\alpha+1}(x)+\frac{1}{2}\mathbf{A}_{1,0}\mathbf{B}_{1,1}T_{\beta+\alpha+2}(x).\end{aligned}$$ Now, note the following regrading $p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$ in (\[eq:expandab2\]):

(i) $\frac{1}{2}\left(\mathbf{A}_{0,0}\mathbf{B}_{0,0}+\mathbf{A}_{0,1}\mathbf{B}_{1,0}\right)$ is the coefficient of $T_{1}(x)$,

(ii) $\frac{1}{2}\left(\mathbf{A}_{1,0}\mathbf{B}_{0,0}+\mathbf{A}_{1,1}\mathbf{B}_{1,0}\right)$ is the coefficient of $T_{\alpha+1}(x)$,

(iii) $\frac{1}{2}\left(\mathbf{A}_{0,0}\mathbf{B}_{0,1}+\mathbf{A}_{0,1}\mathbf{B}_{1,1}\right)$ is the coefficient of $T_{\beta+1}(x)$,

(iv) $\frac{1}{2}\left(\mathbf{A}_{1,0}\mathbf{B}_{0,1}+\mathbf{A}_{1,1}\mathbf{B}_{1,1}\right)$ is the coefficient of $T_{\beta+\alpha+1}(x)$.

[ Since $p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$ has degree $\beta+\alpha+2$, and this polynomial is evaluated at distinct value at each worker node, once the fusion node receives the output of any $\beta+\alpha+3$ worker nodes, it can interpolate $p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$ and extract the product $\mathbf{A}\mathbf{B}$ (i.e., the matrix coefficients of $T_1(x),T_{\alpha+1}(x)$, $ T_{\beta+1}(x), T_{\beta+\alpha+1}(x)$). Now, we aim for picking values for $\alpha,\beta$ such that the degree of $p_\mathbf{A}(x)p_\mathbf{B}(x)$ is minimal; and hence, the recovery threshold is minimal as well. These minimal values for $\alpha,\beta$ must be chosen such that the desired coefficients in (i)-(iv) are separate. That is, each of them is neither combined with another desired nor undesired term. This constraint leads to the following two inequalities: ]{} $$\alpha-1>1,\text{~and }
\alpha+1<\beta-\alpha-1,$$ which implies that $\alpha=3,\beta=9$. Next, we provide our general code construction for the

 {#sec:Chebyshev_polydot}

\[thm:genCheb\] For the matrix multiplication problem described in Section \[sec:probfor2\] computed on the system defined in Section \[sec:sysmod\], there exists a coding strategy with recovery threshold $$\begin{aligned}
\label{eq:recgen}
\hspace{-0mm}4m_1m_2&m_3-2(m_1m_2+m_2m_3+m_3m_1)\notag\\
&+m_1+2m_2+m_3-1.\end{aligned}$$

Notice that the problem specified in Section \[sec:probfor2\] restricts the output matrix of each worker node to be of dimension $N_1/m_1 \times N_3/m_3$, for some positive integers $m_1,m_3$ that divide $N_1,N_3$, respectively. This is smaller than the dimensions of the output matrix of each worker node according to the problem specified in Section \[sec:probfor\] (i.e., $N_1\times N_3$) by a factor of $m_1m_3$. However, according to Theorem \[thm:genCheb\],

For ease of exposition in the remaining of this section, we use $T^{'}_0,T^{'}_1,T^{'}_2, \cdots$ to denote $\frac{1}{2}T_0,T_1,T_2, \cdots$, respectively.

In order to prove Theorem \[thm:genCheb\], we first present a code construction that achieves the recovery threshold in (\[eq:recgen\]), then we prove that the presented code construction is valid. First, note that in code construction, we assume that the two input matrices $\mathbf{A}, \mathbf{B}$ are split as in (\[eq:spltA2\]). Also, note that given this partitioning of input matrices, we can write $\mathbf{C}=\mathbf{A}\mathbf{B}$, where $\mathbf{C}$ is written as $$\begin{aligned}
\mathbf{C}=    \left(\begin{array}{ccc}\mathbf{C}_{0,0}&  \cdots& \mathbf{C}_{0,m_3-1}\\\vdots& \ddots &\vdots\\\mathbf{C}_{m_1-1,0} &\cdots &\mathbf{C}_{m_1-1,m_3-1}\end{array}\right),\end{aligned}$$ and each of $\mathbf{C}_{i,l}$ has dimension ${N_1}/{m_1}\times{N_3}/{m_3}$ and can be expressed as $\mathbf{C}_{i,l}=\sum_{j=0}^{m_2-1} \mathbf{A}_{i,j}\mathbf{B}_{j,l},$ for any $i\in \{0, 1, \cdots, m_1-1\},$ and $l\in \{0,1, \cdots, m_3-1\}$. Also, let $x_1, \cdots, x_P$ be distinct real numbers in the range $[-1,1]$, and define encoding polynomials $$\begin{aligned}
\label{eq:genpolys}
p_\mathbf{A}(x)&=\sum_{i=0}^{m_1-1}\sum_{j=0}^{m_2-1}\mathbf{A}_{i,j} T^{'}_{m_2-1-j+i(2m_2-1)}(x), \notag\\ p_\mathbf{B}(x)&=\sum_{k=0}^{m_2-1}\sum_{l=0}^{m_3-1} \mathbf{B}_{k,l} T^{'}_{k+l(2m_1-1)(2m_2-1)}(x),\end{aligned}$$ and let $p_\mathbf{C}(x)=p_\mathbf{A}(x)p_\mathbf{B}(x)$. Notice that $p_{\mathbf{C}}(x)$ is a polynomial matrix of degree equals $\deg_\mathbf{C}:=4m_1m_2m_3-2(m_1m_2+m_2m_3+m_3m_1)$ $+m_1+2m_2+m_3-2$.

\[cl:chebcoef\] For any $i\in \{0, 1, \cdots, m_1-1\}$ and $l\in \{0,1, \cdots, m_3-1\}$, $\frac{1}{2}\mathbf{C}_{i,l}$ is the matrix coefficient of $T_{m_2-1+i(2m_2-1)+l(2m_1-1)(2m_2-1)}$ in $p_{\mathbf{C}}(x)$,

The proof of this claim is in Appendix \[app:proof\].

We describe, next, the idea of our proposed . First, for all $r \in [P]$, the master node sends to the $r$-th worker evaluations of $p_\mathbf{A}(x)$ and $p_\mathbf{B}(x)$ at $x=\rho^{(P)}_r$, that is, it sends $p_\mathbf{A}(\rho^{(P)}_r)$ and $p_\mathbf{B}(\rho^{(P)}_r)$ to the $r$-th worker. Next, for every $r \in [P]$, the $r$-th worker node computes the matrix product $p_\mathbf{C}(\rho^{(P)}_r)=p_{\mathbf{A}}(\rho^{(P)}_r)p_{\mathbf{B}}(\rho^{(P)}_r)$ and sends the result to the fusion node. Once the fusion node receives the output of any $\deg_\mathbf{C}+1$ worker nodes, it interpolates $p_\mathbf{C}(x)$.

We formally present our in Construction \[con:gencheb\].

$r \gets 1$ $p_\mathbf{A}(\rho^{(P)}_r) \gets \sum_{i=0}^{m_1-1}\sum_{j=0}^{m_2-1}\mathbf{A}_{i,j} T^{'}_{m_2-1-j+i(2m_2-1)}(\rho^{(P)}_r)$ $p_\mathbf{B}(\rho^{(P)}_r)\gets\sum_{k=0}^{m_2-1}\sum_{l=0}^{m_3-1} \mathbf{B}_{k,l} T^{'}_{k+l(2m_1-1)(2m_2-1)}(\rho^{(P)}_r)$ **send** $p_{\mathbf{A}}(\rho^{(P)}_r), p_{\mathbf{B}}(\rho^{(P)}_r)$ **to worker node** $r$ $r \gets r+1$ $p_\mathbf{C}(\rho^{(P)}_r) \gets p_{\mathbf{A}}(\rho^{(P)}_r) p_{\mathbf{B}}(\rho^{(P)}_r)$ **send** $p_{\mathbf{C}}(\rho^{(P)}_r)$ **to the fusion node**

$\mathbf{G}_{\operatorname{inv}}\gets \left(\mathbf{G}^{(\deg_\mathbf{C}+1,P)}_{\mathcal{R}}\right)^{-1}$ $(c^{(i,j)}_0, \cdots, c^{(i,j)}_{\deg_\mathbf{C}}) \gets (p_{\mathbf{C}}^{(i,j)}(\rho^{(P)}_{r_1}), \cdots,p_{\mathbf{C}}^{(i,j)}(\rho^{(P)}_{r_{\deg_\mathbf{C}+1}}))\mathbf{G}_{\operatorname{inv}}$ $\hat{c}_{k,l}^{(i,j)} \gets 2 c^{(i,j)}_{m_2-1+(k-1)(2m_2-1)+(l-1)(2m_1-1)(2m_2-1)}$ **return** $\hat{\mathbf{C}}$

Now, we prove Theorem \[thm:genCheb\].

To prove the theorem, it suffices to prove that Construction \[con:gencheb\] is valid. Noting that $p_\mathbf{A}(x)p_\mathbf{B}(x)$ has degree $4m_1m_2m_3-2(m_1m_2+m_2m_3+m_3m_1)+m_1+2m_2+m_3-2$ and every worker node sends an evaluation of $p_\mathbf{A}(x)p_\mathbf{B}(x)$ at a distinct point, once the fusion node receives the output of any $4m_1m_2m_3-2(m_1m_2+m_2m_3+m_3m_1)+m_1+2m_2+m_3-1$ worker node, it can interpolate $p_\mathbf{A}(x)p_\mathbf{B}(x)$ (i.e., obtain all its matrix coefficients). This includes the coefficients of $T_{m_2-1+i(2m_2-1)+l(2m_1-1)(2m_2-1)}$ for all $i\in \{0, 1, \cdots, m_1-1\},$ and $l\in \{0,1, \cdots, m_3-1\}$, i.e., $\mathbf{C}_{i,l}$, for all $i\in \{0, 1, \cdots, m_1-1\},$ and $l\in \{0,1, \cdots, m_3-1\}$ (Claim \[cl:chebcoef\]), which completes the proof.

Next, we provide the different complexity analyses of

### Complexity Analyses of  {#sec:con3complexity}

**Encoding Complexity**: Encoding for each worker requires performing two additions, the first one adds $m_1m_2$ scaled matrices of size $N_1 N_2/ (m_1m_2)$ and the other adds $m_2 m_3$ scaled matrices of size $N_2 N_3 /(m_2m_3)$, for an overall encoding complexity for each worker of $O(N_1 N_2+N_2N_3)$. Therefore, the overall computational complexity of encoding for $P$ workers is $O(N_1 N_2P+N_2N_3P)$.

**Computational Cost per Worker**: Each worker multiplies two matrices of dimensions $N_1  /m_1\times N_2/m_2$ and $N_2/m_2\times N_3/m_3$, requiring $O(N_1N_2N_3/(m_1m_2m_3))$ operations.

**Decoding Complexity**: Since $p_\mathbf{A}(x)p_\mathbf{B}(x)$ has degree $k-1:=4m_1m_2m_3- 2(m_1m_2+m_2m_3+m_3m_1)+m_1+2m_2+m_3-2$, the interpolation of $p_\mathbf{C}(x)$ requires the inversion of a $k \times k$ matrix, with complexity $O(k^3)=O(m_1^3m_2^3m_3^3)$, and performing $N_1N_3/(m_1m_3)$ matrix-vector multiplications, each of them is between the inverted matrix and a column vector of length $k$ of the received evaluations of the matrix polynomial $p_\mathbf{C}(x)$ at some position $(i,j) \in [N_1/m_1] \times [N_3/m_3]$, with complexity $O(N_1N_3k^2/(m_1m_3))=O(N_1N_3m_1m_2^2m_3)$. Thus, assuming that $m_1,m_3 \ll N_1,N_3$, the overall decoding complexity is $O(N_1N_3m_1m_2^2m_3)=O(N_1N_3 m n)$.

**Communication Cost**: The master node sends $O( N_1N_2P/(m_1m_2)+N_2N_3P/(m_2m_3))$ symbols, and the fusion node receives $O(N_1 N_3 m_2)$ symbols from the successful worker nodes.

With the reasonable assumption that the dimensions of the input matrices $\mathbf{A},\mathbf{B}$ are large enough such that $N_1,N_2,N_3 \gg m_1,m_2,m_3,P$, we can conclude that the encoding and decoding costs at the master and fusion nodes, respectively, are negligible compared to the computation cost at each worker node.

![image](Pall_gen-eps-converted-to.pdf)

------------------------------------------------------------------------

Numerical Results {#sec:numerical_polydot}
-----------------

In our experiments on Construction \[con:gencheb\], we considered distributed systems with $P=16,25$ worker nodes. Fig. \[fig:Pallgen\] shows that, for every examined system, the condition number of the interpolation matrix using the Codes is less than its counterpart codes in [@genPolyDot; @entPoly]. The results in Fig. \[fig:Pallgen\] also show that, for the same system, as the partitioning factor $m_1$ decreases (i.e., as the redundancy in worker nodes increases), the stability of decreases; however, it is still better than the monomial-basis based codes in any cases.

Numerically Stable Lagrange Coded Computing {#sec:lcc}
===========================================

In this section, we study the numerical stability of Lagrange coded computing [@yu2018lagrange] that lifts coded computing beyond matrix-vector and matrix-matrix multiplications, to multi-variate polynomial computations. As shown in [@yu2018lagrange], Lagrange coded computing has applications in gradient coding, privacy and secrecy. Our main contribution here is to develop a numerically stable approach towards Lagrange coded computing inspired by our result of Theorem \[thm:bound\]. In particular, our contribution involves (a) careful choice of evaluation points, and (b) a careful decoding algorithm that involves inversion of the appropriate Chebyshev Vandermonde matrix. We describe the system model in Section \[sec:lccsysprob\]. We overview the Lagrange coded computing technique of [@yu2018lagrange] in Section \[sec:Lagrange\]. We describe our numerically stable approach in Section \[sec:Lagrange\_stable\], and present the results of our numerical experiments in Section \[sec:Lagrange\_numerical\].

System Model and Problem Formulation {#sec:lccsysprob}
------------------------------------

We consider, for this section, the distributed computing framework depicted in Fig. \[fig:lccsys\], that is used in [@yu2018lagrange] and consists of a master node, $P$ worker nodes, and a fusion node where the only communication allowed is from the master node to the different worker nodes and from the worker nodes to the fusion node. The worker nodes have a prior knowledge of a polynomial function of interest $f: \mathbb{R}^d \rightarrow \mathbb{R}^v$ of degree $\operatorname{deg}(f)$, where $d,v\in\mathbb{N}^{+}$. In addition, the master node possesses a set of data points $\mathcal{X}=\{X_1, \cdots, X_m\}$, where $X_i \in \mathbb{R}^d$, $i \in [m]$. For every worker node $i \in [P]$, the master node is allowed to send some encoded vector $\tilde{X}_i(X_1,\cdots, X_m) \in \mathbb{R}^d$. Once a worker node receives the encoded vector on its input, it evaluates $f$ at this encoded vector and sends the evaluation to the fusion node. That is, for $i\in [P]$, worker node $i$ receives $\tilde{X}_i$ on its input, evaluates $f(\tilde{X}_i)$, then it sends the result to the fusion node. Finally, the fusion node is expected to numerically stably decode the set of evaluations $\mathcal{F}=\{f(X_1), \cdots, f(X_m)\}$ after it receives the output of any $K$ worker nodes.

![The Lagrange coded computing system framework[]{data-label="fig:lccsys"}](Model_LagCC.pdf)

Background on Lagrange Coded Computing {#sec:Lagrange}
--------------------------------------

In this section, we review the baseline Lagrange coded computing method introduced in [@yu2018lagrange] considering the framework in Section \[sec:lccsysprob\]. Notice that although the method in [@yu2018lagrange] is more general, here, for simplicity, we limit our discussion to the *systematic* Lagrange coded computing. That is, we assume that for $i\in [m]$, worker node $i$ receives the $i$-th data point from the master node. In other words, we assume that $\tilde{X}_i=X_i, i\in[m]$. Now, the encoding procedure goes as follows: First, let $x_1, \cdots, x_P$ be distinct real values, an encoding function $g(x)$ is defined as: $$\begin{aligned}
\label{eq:lagenc}
g(x)=\sum_{i=1}^{m} X_i \prod_{j\in[m]-i} \frac{x-x_j}{x_i-x_j}.\end{aligned}$$ Given this encoding function, the master node sends the encoded vector $\tilde{X}_i=g(x_i)$ to the worker node $i$, for every $i \in [P]$. Notice that the encoding function $g(x)$ indeed leads to a systematic encoding since $\tilde{X}_i=g(x_i)=X_i,$ for all $i\in [m]$. Every worker node $i$ computes $f(\tilde{X}_i)$ upon the reception of $\tilde{X}_i$, and sends the result to the fusion node. The fusion node waits till receiving the output of any $K:=(m-1)\deg(f)+1$.

Numerically Stable Lagrange Coded Computing {#sec:Lagrange_stable}
-------------------------------------------

Our numerically stable algorithm for Lagrange coded computing is formally described in Construction \[con:Lag\]. Finally, we assume that our construction returns as output the set of evaluations $\hat{\mathcal{F}}=$ $\{\hat{f}(X_1),$ $\cdots, \hat{f}(X_m)\}$, where for each $\hat{f}(X_i), i\in [m]$, we have $\hat{f}(X_i)= (\hat{f}^{(1)}(x_i), \cdots, \hat{f}^{(v)}(x_i))$, where for every $i\in[m], j \in [v], \hat{f}^{(j)}(x_i)$ and ${f}^{(j)}(x_i)$ would be the same if the machine had infinite precision.

In the following, we show through numerical experiments the stability of our proposed Construction \[con:Lag\].

$r \gets 1$ $\tilde{X}_r \gets X_i$ $\tilde{X}_r \gets\sum_{i=1}^{m} X_i \prod_{j\in[m]-i} \frac{\rho^{(P)}_r-\rho^{(P)}_j}{\rho^{(P)}_i-\rho^{(P)}_j}$ **send** $\tilde{X}_r$ **to worker node** $r$ $r \gets r+1$ ${\text{Out}}_r \gets f(\tilde{X}_r)$ **send** ${\text{Out}}_r$ **to the fusion node** $\mathbf{G}_{\text{inv}} \gets \left(\mathbf{G}_{\mathcal{R}}^{(K,P)}\right)^{-1}$ $(c^{(i)}_0, \cdots, c^{(i)}_{K-1}) \gets ({\text{Out}}^{(i)}_{{r_1}}, \cdots,{\text{Out}}^{(i)}_{{r_{K}}})\mathbf{G}_{\operatorname{inv}}$ $(\hat{f}^{(i)}(x_1), \cdots, \hat{f}^{(i)}(x_m))\gets (c^{(i)}_0, \cdots, c^{(i)}_{K-1}) \mathbf{G}^{(K,P)}_{[m]}$ **return** $\hat{\mathcal{F}}$

Numerical Results {#sec:Lagrange_numerical}
-----------------

In our experiments, we assume that we have a distributed system of $P$ worker nodes, $m=P-2$ data points/input vectors $X_1, \cdots, X_m$, each of them is of dimension $d=10$, where each entry of every input vector is picked independently, according to the standard Gaussian distribution $\mathcal{N}(0,1)$. The function of interest in this system is $f(X)=Y^TX$, where $Y$ is some $d$-dimensional vector with entries picked independently according to the standard Gaussian distribution $\mathcal{N}(0,1)$. In our experiments, we compare between Construction \[con:Lag\], where the Chebyshev basis is used for interpolation, and the case where the monomial basis is used for interpolation instead. Let $\hat{\mathbf{f}}=(\hat{f}(X_1) \cdots \hat{f}(X_m))$ be the system’s output vector, and ${\mathbf{f}}=({f}(X_1) \cdots {f}(X_m))$ be the correct output vector, we define the relative error between $\mathbf{f}$ and $\hat{\mathbf{f}}$ to be $$\begin{aligned}
    E_r(\mathbf{f},\hat{\mathbf{f}})=\frac{||\mathbf{f}-\hat{\mathbf{f}}||_2}{||\mathbf{f}||_2}.\end{aligned}$$

![image](Lag_fixed_par1.pdf)

------------------------------------------------------------------------

The results, shown in Fig. \[fig:lagfixpar\], illustrates that using the Chebyshev basis for interpolation provides less relative error/higher stability than the monomial basis at every system size. Fig. \[fig:lagfixpar\] also shows that under a certain relative error constraint, Construction \[con:Lag\] provides higher scalability than the monomial basis case. Specifically, let us assume that a relative error up to $0.1$ can be tolerated, Fig. \[fig:lagfixpar\] shows that the monomial-basis interpolation construction can support systems with a number of worker nodes only less than $40$. However, for the same relative error constraint, Construction \[con:Lag\] can support systems with a number of worker nodes up to $100$.

Concluding Remarks {#sec:conc}
==================

In this paper, we develop numerically stable codes for matrix-matrix multiplication and Lagrange coded computing. A distinctive character of our work is the infusion of principles of numerical approximation theory into coded computing towards the end goal of numerical stability. In particular, our work is marked by the use of orthogonal polynomials for encoding, Gauss quadrature techniques for decoding and new bounds on the condition number of Chebyshev Vandermonde matrices. Notably, our constructions obtain the same recovery threshold as MatDot Codes and Polynomial Codes for matrix multiplication as well as for Lagrange Coded Computing. However, our construction in Section \[sec:polydot-chebyshev\] obtains a weaker (higher) recovery threshold than previous constructions [@entPoly; @genPolyDot] for the problem of coded matrix multiplication when the computation/communication cost is constrained to be lower than that of MatDot Codes. The search of numerically stable codes for this application with the same recovery threshold as [@entPoly; @genPolyDot] remains open.

While our paper focuses on applications where polynomial based encoding are particularly useful, our results might be useful for other applications as well. For instance, for the simple matrix-multiplication problem $\mathbf{A}\mathbf{x}$ performed in a distributed setting over $P$ worker nodes, where the goal is to encode $\mathbf{A}$ such that each worker stores a partition $1/m$ of matrix $\mathbf{A},$ it is well known that MDS type codes can be used [@Huang_TC_84; @lee2016speeding]. Specifically, let $\mathbf{A}=\begin{bmatrix}\mathbf{A}_{1} \\ \mathbf{A}_{2} \\ \vdots \\ \mathbf{A}_{m}\end{bmatrix}$ and let $\mathbf{H}=(h_{ij})$ be an $m \times P$ matrix where every $m \times m$ submatrix of $\mathbf{H}$ has a full rank of $m$. Then the $p$-th worker for $p \in \{1,2,\ldots,P\}$ can compute $\left(\sum_{i=1}^{m} {h}_{ip} \mathbf{A}_{i}\right)\mathbf{x};$ the product $\mathbf{A}\mathbf{x}$ can be recovered from any $m$ of the $P$ nodes. The instinctual, Reed-Solomon inspired solution of choosing $\mathbf{H}$ to be a Vandermode matrix is ill-conditioned over real numbers. Note however that, unlike the matrix multiplication problem, the matrix $\mathbf{H}$ does not need to have a polynomial structure. Indeed, choosing $\mathbf{H}$ to be a random Gaussian matrix leads to well-conditioned solutions with high probability. In particular, the following result follows from elementary arguments that build on [@azais2004].

Let $\mathbf{H}$ be an $m \times P$ matrix, $P\geq m\geq 3$, and let the entries of $\mathbf{H}$ be independent and identically distributed standard Gaussian random variables. Then, $$\operatorname{Pr}\Big(\kappa_2^{max}(\mathbf{H}) > m P^{2(P-m)}\big)\Big) <\frac{5.6}{P^{(P-m)}}.$$ \[thm:iidG\]

The theorem which is proved in Appendix \[app:iidG\], formally demonstrates that for a fixed number of redundant workers $s=P-m,$ the worst case condition number grows as $O(mP^{2s})$ with high probability. However, the random Gaussian matrix approach has two drawbacks: (i) for a given realization of the random variables, it is difficult to verify whether it is well-conditioned, and (ii) the lack of structure could lead to more complex decoding. Our result of Theorem \[thm:bound\] also indicates that choosing $\mathbf{H}=\mathbf{G}^{(m,P)}(\boldsymbol{\rho}^{(P)}),$ i.e., to be a Chebyshev Vandermonde matrix, naturally provides a well-conditioned solution to this problem. Another solution for the matrix-vector multiplication problem is provided in [@ramamoorthy2019universally] via universally decodable matrices [@ganesan2007existence]; in this work numerical stability is demonstrated empirically.

It is, however, important to note that the problems resolved in our paper here are more restrictive since matrix multiplication codes - where both matrices are to be encoded so that the product can be recovered - require much more structure than matrix-multiplication where only one matrix is to be encoded. For instance, random Gaussian encoding does not naturally work for matrix multiplication to get a recovery threshold of $2m-1$, and it is not clear whether the solution of [@ramamoorthy2019universally] is applicable either. The utility of Chebyshev-Vandermonde matrices for a variety of coded computing problems including matrix-vector multiplication, matrix multiplication and Lagrange coded computing motivates the study of low-complexity decoding and error correction mechanisms for these systems.

Proof of Claim \[cl:orthAB\] {#app:ortho}
============================

We have, $$\begin{aligned}
\label{eq:clorth}
\int_{a}^{b} p_{\mathbf{A}}(x) p_{\mathbf{B}}(x) w(x) dx 
&=\int_{a}^{b} \left(\sum_{i=0}^{m-1} \mathbf{A}_i q_{i}(x) \right)\left(\sum_{j=0}^{m-1} \mathbf{B}_j q_{j}(x)\right) w(x) dx\notag\\
&=\int_{a}^{b} \sum_{i=0}^{m-1} \sum_{j=0}^{m-1} \mathbf{A}_i \mathbf{B}_j q_{i}(x)  q_{j}(x) w(x) dx\notag\\
&= \sum_{i=0}^{m-1} \sum_{j=0}^{m-1} \mathbf{A}_i \mathbf{B}_j\int_{a}^{b} q_{i}(x)  q_{j}(x) w(x) dx\notag\\
&=\sum_{i=0}^{m-1} \sum_{j=0}^{m-1} \mathbf{A}_i \mathbf{B}_j\ \langle q_{i},  q_{j}\rangle \notag\\
&=\sum_{i=0}^{m-1} \mathbf{A}_i \mathbf{B}_i\  \notag\\
&=\mathbf{A}\mathbf{B}.\end{aligned}$$ In addition, noting that $p_\mathbf{A}(x)p_\mathbf{B}(x)$ (i.e., $p_\mathbf{C}(x)$) is of degree $2m-2$ (less than $2m$), Theorem \[thm:GQ\] implies that $$\begin{aligned}
\label{eq:clorth2}
\int_{a}^{b}   p_\mathbf{A}(x)p_\mathbf{B}(x) w(x) dx &= \sum_{r=1}^{m} a_r p_{\mathbf{A}}(\eta_r) p_{\mathbf{B}}(\eta_r)    \notag\\
&=\sum_{r=1}^{m} a_r p_{\mathbf{C}}(\eta_r). \end{aligned}$$ Finally, combining (\[eq:clorth\]) and (\[eq:clorth2\]) completes the proof.[$\Box$]{}

Proof of Theorem \[thm:bound\] {#app:proofgen}
==============================

We use the following trigonometric identity in our proof.

For $n \geq 0$, let $x_i$ be chosen as (\[eq:Chebnodes\]). Then $\prod_{j \neq i} (x_i - x_j) = (-1)^{i-1} \frac{2^{1-n}n}{\sin(\frac{(2i-1) \pi}{2n})}$ \[lem:trig\]

Note that $2^{n-1}\prod_{i=1}^{n}(x-x_i) = T_n(x) =  \cos(n\cos^{-1}(x))$. Therefore, $$2^{n-1}\prod_{j \neq i} (x_i - x_j) = T_n'(x_i) = \frac{n}{ \sqrt{1-x_i^2}} \sin(n \cos^{-1}(x_i))$$ where $T_n'(x)$ denotes the derivative of $T_n(x).$ Using $x_i=\cos(\frac{(2i-1) \pi}{2n})$ above we get the desired result.

[We show that any square sub-matrix of $\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})$ formed by any $n-s$ columns of $\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})$ satisfies the bound stated in the theorem.]{} [Let $\mathcal{S}$ be a subset of $[n]$ such that $|\mathcal{S}|=s$, for some $s<n$, and define $\mathbf{G}^{(n-s,n)}_{[n]-\mathcal{S}}(\boldsymbol{\rho}^{(n)})$ to be the square $n-s \times n-s$ submatrix of $\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})$ after removing the columns with indices in $\mathcal{S}$. Recalling the structure of $\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})$ from (\[eq:GT\]), we can write it as $$\begin{aligned}
  \mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})= \left(\hspace{-2mm}\begin{array}{ccc} T_0(\rho^{(n)}_1)&  \cdots & T_{0}(\rho^{(n)}_n)\\\vdots & \ddots &\vdots \\ 
   T_{n-s-1}(\rho^{(n)}_{1})&\cdots& T_{n-s-1}(\rho^{(n)}_{n})\end{array}\hspace{-2mm}\right).\notag\end{aligned}$$ Moreover, for any $\mathcal{S}\subset[n]$ such that $|\mathcal{S}|=s$, we can write $$\begin{aligned}
 \mathbf{G}_{[n]-\mathcal{S}}^{(n-s,n)}(\boldsymbol{\rho}^{(n)}):= \mathbf{G}_{\Gamma}^{(n-s,n)}:= \left(\hspace{-2mm}\begin{array}{ccc} T_0(\gamma_1)&  \cdots & T_{0}(\gamma_{n-s})\\\vdots & \ddots &\vdots \\ 
   T_{n-s-1}(\gamma_{1})&\cdots& T_{n-s-1}(\gamma_{n-s})\end{array}\hspace{-2mm}\right),\notag\end{aligned}$$ where ${\Gamma}=(\gamma_1, \gamma_2, \cdots, \gamma_{n-s})=(\rho^{(n)}_{g_1}, \rho^{(n)}_{g_2}, \cdots, \rho^{(n)}_{g_{n-s}})$, where $\{g_i\}_{i\in[n-s]}=[n]-\mathcal{S}$ and $g_1 < g_2 < \cdots < g_{n-s}$.]{} [Now, notice that $||\mathbf{G}_{\Gamma}^{(n-s,n)}||^2_{F}=\sum_{i=1}^{n-s}\sum_{j=1}^{n-s} |T_{i-1}(\gamma_j)|^2$, and $|T_i(\gamma_j)| \leq 1$ for any $i,j \in [n-s]$. Therefore, we have $$\begin{aligned}
\label{eq:GT2}
||\mathbf{G}_{\Gamma}^{(n-s,n)}||^2_{F} \leq (n-s)^2.\end{aligned}$$ In the following, we obtain an upper bound on $||(\mathbf{G}_{\Gamma}^{(n-s,n)})^{-1}||_{F}$. Let $L_{\Gamma,k}$ be the $k$-th Lagrange polynomial associated with $\Gamma$, that is, $$\begin{aligned}
\label{eq:kLag}
L_{\Gamma,k}(x)= \prod_{i\in[n-s]-\{k\}} \frac{x-\gamma_i}{\gamma_k-\gamma_i}     \end{aligned}$$ Since $L_{\Gamma,k}(x)$ has a degree of $n-s-1$, it can be written in terms of the Chebyshev basis $T_{0}(x), \cdots, T_{n-s-1}(x)$ as $$\begin{aligned}
\label{eq:kLagCh}
L_{\Gamma,k}(x)=\sum_{i=0}^{n-s-1} a_{i,k} T_{i}(x),    \end{aligned}$$ for some real coefficients $a_{0,k}, \cdots, a_{n-s-1,k}$. Now, from (\[eq:kLag\]), note the following property regarding $L_{\Gamma,k}(x)$: $$\begin{aligned}
L_{\Gamma,k}(x)=\left\{\begin{array}{cl} 1, &\text{if } x=\gamma_k\\0, &\text{if } x\in\{\gamma_i\}_{i\in[n-s]-{k}}.\end{array}\right. \notag    \end{aligned}$$ Using this property and observing (\[eq:kLagCh\]), we conclude that, for any $j \in [n-s]$, $\sum_{i=0}^{n-s-1} a_{i,k} T_{i}(\gamma_j)=\delta(k-j)$. Therefore, $$\begin{aligned}
\left(
\begin{array}{ccc}
a_{0,1}    &\cdots & a_{n-s-1,1}  \\
\vdots    &\ddots & \vdots \\
a_{0,n-s}    & \cdots & a_{n-s-1,n-s}
\end{array}\right)    \mathbf{G}_{\Gamma}^{(n-s,n)}=\mathbf{I}_{n-s},\notag\end{aligned}$$ where $\mathbf{I}_{n-s}$ is the $n-s\times n-s$ identity matrix. That is, $$\begin{aligned}
\left(\mathbf{G}_{\Gamma}^{(n-s,n)}\right)^{-1}=\left(
\begin{array}{ccc}
a_{0,1}    &\cdots & a_{n-s-1,1}  \\
\vdots    &\ddots & \vdots \\
a_{0,n-s}    & \cdots & a_{n-s-1,n-s}
\end{array}\right),\end{aligned}$$ Therefore, $$\begin{aligned}
\label{eq:inv1}
\Big| \Big|\left(\mathbf{G}_{\Gamma}^{(n-s,n)}\right)^{-1}\Big|\Big|_F^2=\sum_{i=1}^{n-s}\sum_{j=1}^{n-s} |a_{i-1,j}|^2.\end{aligned}$$ In addition, we have that $$\begin{aligned}
\label{eq:inv2}
\sum_{k=1}^{n-s} \int_{-1}^{1} L^2_{\Gamma,k}(x) w(x )dx&=  \sum_{k=1}^{n-s} \int_{-1}^{1} \sum_{i=0}^{n-s-1}\sum_{j=0}^{n-s-1} a_{i,k}a_{j,k} T_{i}(x)T_{j}(x) w(x)dx\notag\\
&=\sum_{k=1}^{n-s}  \sum_{i=0}^{n-s-1}\sum_{j=0}^{n-s-1} a_{i,k}a_{j,k}\int_{-1}^{1} T_{i}(x)T_{j}(x) w(x)dx\notag\\
&=\sum_{k=1}^{n-s}  \sum_{i=0}^{n-s-1}\sum_{j=0}^{n-s-1} a_{i,k}a_{j,k} \langle T_{i}, T_{j}\rangle \notag\\
&=\sum_{k=1}^{n-s}  \sum_{i=0}^{n-s-1} |a_{i,k}|^2.\end{aligned}$$ From (\[eq:inv1\]) and (\[eq:inv2\]), we conclude that $||(\mathbf{G}_{\Gamma}^{(n-s,n)})^{-1}||_F^2=\sum_{k=1}^{n-s} \int_{-1}^{1} L^2_{\Gamma,k}(x) w(x )dx$.]{}

Now, we express the integral $\int_{-1}^{1} L^2_{\Gamma,k}(x) w(x )dx$ in the Gauss quadrature form using the $n$ roots of $T_n(x):$ $\rho_1^{(n)}, \cdots, \rho_n^{(n)}$. Note that this is a “trick” we use in the proof - it is possible to use the Gauss quadrature formula over $n-s$ nodes to express the integral of the degree $2(n-s-1)$ polynomial $L^2_{\Gamma,k}(x)$. However, the use of $n$ nodes instead of $n-s$ nodes leads to simple tractable bound for $||(\mathbf{G}_{\Gamma}^{(n-s,n)})^{-1}||^2_F.$ Now, we can write $$\begin{aligned}
    \int_{-1}^{1} L^2_{\Gamma,k}(x) w(x )dx=\sum_{i=1}^{n} c_i L_{\Gamma,k}^2(\rho_i^{(n)}),\end{aligned}$$ for some constants $c_1, \cdots, c_n$. Moreover, $c_1, \cdots, c_n$ for the Chebyshev polynomials of the first kind are, in fact, all equal to $\pi /n$. Therefore, we have $$\begin{aligned}
\label{eq:kthLSum}
    \int_{-1}^{1} L^2_{\Gamma,k}(x) w(x )dx=\frac{\pi}{n}\sum_{i=1}^{n}  L_{\Gamma,k}^2(\rho_i^{(n)}),\end{aligned}$$ and, consequently, $$\begin{aligned}
\label{eq:kthLSum2}
    \Big|\Big|\left(\mathbf{G}_{\Gamma}^{(n-s,n)}\right)^{-1}\Big|\Big|_F^2=\frac{\pi}{n}\sum_{k=1}^{n-s} \sum_{i=1}^{n} L^2_{\Gamma,k}(\rho_i^{(n)}).\end{aligned}$$ Now, from (\[eq:kLag\]), note that $L_{\Gamma,k}(x)$ has the following evaluations $$\begin{aligned}
L_{\Gamma,k}(\rho_{i}^{(n)})=\left\{\begin{array}{cl} 1, &\text{if } i=g_k\\0, &\text{if } i\in\{g_i\}_{i\in[n-s]-{k}}\\\prod_{j\in[n-s]-\{k\}} \frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}   , &\text{if } i\in \mathcal{S}.\end{array}\right. \end{aligned}$$ Therefore, (\[eq:kthLSum2\]) can be written as $$\begin{aligned}
\label{eq:kthLSum3}
    \Big|\Big|\left(\mathbf{G}_{\Gamma}^{(n-s,n)}\right)^{-1}\Big|\Big|_F^2&=\frac{\pi}{n}\sum_{k=1}^{n-s}\left( 1+\sum_{i\in \mathcal{S}} \prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}\right)^2\right)\notag\\
    &=\frac{\pi(n-s)}{n}+\frac{\pi}{n}\sum_{k=1}^{n-s}\sum_{i\in \mathcal{S}} \prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}\right)^2\end{aligned}$$ In order to obtain our upper bound on $||(\mathbf{G}_{\Gamma}^{(n-s,n)})^{-1}||_F^2$, in the following, we get an upper bound on the term $\prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}\right)^2$ in (\[eq:kthLSum3\]). Notice that $\prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}\right)^2$ can be written as $$\begin{aligned}
\label{eq:product}
    \prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}\right)^2&= \prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\rho_{g_j}^{(n)}}{\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}}\right)^2\notag\\
    &=\left[\prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\rho_{g_j}^{(n)}}{\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}}\right)^2\right] \frac{\prod_{j\in \mathcal{S}\cup \{g_k\}-\{i\}}\left(\rho_{i}^{(n)}-\rho_{j}^{(n)}\right)^2}{\prod_{j\in \mathcal{S}\cup \{g_k\}-\{i\}}\left(\rho_{i}^{(n)}-\rho^{(n)}_{j}\right)^2}\notag\\&=\frac{\prod_{j\in[n]-\{i\}}\left(\rho_{i}^{(n)}-\rho_{j}^{(n)}\right)^2}{\prod_{j\in[n-s]-\{k\}}\left(\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}\right)^2\prod_{j\in \mathcal{S}\cup\{g_k\}-\{i\}}\left(\rho_{i}^{(n)}-\rho^{(n)}_{j}\right)^2}\notag\\
    &\stackrel{}{=}\frac{ \left(2^{1-n}n / {\sin(\frac{(2i-1) \pi}{2n})}\right)^2}{\prod_{j\in[n-s]-\{k\}}\left(\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}\right)^2\prod_{j\in \mathcal{S}\cup\{g_k\}-\{i\}}\left(\rho_{i}^{(n)}-\rho^{(n)}_{j}\right)^2},\end{aligned}$$ where the last equality follows from Lemma \[lem:trig\]. Moreover, the product $\prod_{j\in[n-s]-\{k\}}\left(\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}\right)^2$ in (\[eq:product\]) can be written as $$\begin{aligned}
\label{eq:product2}
    \prod_{j\in[n-s]-\{k\}}\left(\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}\right)^2&=
   \frac{ \prod_{j\in[n-s]-\{k\}}\left(\rho_{g_k}^{(n)}-\rho^{(n)}_{g_j}\right)^2\prod_{j\in \mathcal{S}} \left(\rho^{(n)}_{g_k}-\rho^{(n)}_{j}\right)^2}{\prod_{j\in \mathcal{S}} \left(\rho^{(n)}_{g_k}-\rho^{(n)}_{j}\right)^2}\notag\\
   &= \frac{ \left(2^{1-n}n / {\sin(\frac{(2g_k-1) \pi}{2n})}\right)^2}{\prod_{j\in \mathcal{S}} \left(\rho^{(n)}_{g_k}-\rho^{(n)}_{j}\right)^2},\end{aligned}$$ where the last equality follows from Lemma \[lem:trig\]. Now, substituting from (\[eq:product2\]) in (\[eq:product\]) yields $$\begin{aligned}
\label{eq:product3}
    \prod_{j\in[n-s]-\{k\}} \left(\frac{\rho_{i}^{(n)}-\gamma_j}{\gamma_k-\gamma_j}\right)^2&= \frac{ \left({\sin(\frac{(2g_k-1) \pi}{2n})}\right)^2}{\left({\sin(\frac{(2i-1) \pi}{2n})}\right)^2}\frac{\prod_{j\in \mathcal{S}} \left(\rho^{(n)}_{g_k}-\rho^{(n)}_{j}\right)^2}{\prod_{j\in \mathcal{S}\cup\{g_k\}-\{i\}}\left(\rho_{i}^{(n)}-\rho^{(n)}_{j}\right)^2}\notag\\
    &= \frac{ \left({\sin(\frac{(2g_k-1) \pi}{2n})}\right)^2}{\left({\sin(\frac{(2i-1) \pi}{2n})}\right)^2}\frac{\prod_{j\in \mathcal{S}-\{i\}} \left(\rho^{(n)}_{g_k}-\rho^{(n)}_{j}\right)^2}{\prod_{j\in \mathcal{S}-\{i\}}\left(\rho_{i}^{(n)}-\rho^{(n)}_{j}\right)^2},\notag\\
    &\leq\frac{ \left({\sin(\frac{(2g_k-1) \pi}{2n})}\right)^2}{\left({\sin(\frac{(2i-1) \pi}{2n})}\right)^2}\left[\frac{\max_{j\in \mathcal{S}-\{i\}} \left(\rho^{(n)}_{g_k}-\rho^{(n)}_{j}\right)^2}{\min_{j\in \mathcal{S}-\{i\}}\left(\rho_{i}^{(n)}-\rho^{(n)}_{j}\right)^2}\right]^{s-1}\notag\\
    &= \frac{ 1}{\left({\sin(\frac{(2i-1) \pi}{2n})}\right)^2}\left[\frac{4}{ \left(\cos(\frac{\pi}{2n})-\cos(\frac{3\pi}{2n})\right)^2}\right]^{s-1}\notag\\
    &=\frac{4^{s-1}}{\left({\sin(\frac{(2i-1) \pi}{2n})}\right)^2 \left(\cos(\frac{\pi}{2n})-\cos(\frac{3\pi}{2n})\right)^{2(s-1)}}\notag\\
    &=O(4^{s-1}n^{2+4(s-1)}). \end{aligned}$$ Using (\[eq:product3\]) in (\[eq:kthLSum3\]), we conclude that $$\begin{aligned}
\label{eq:invG}
    \Big|\Big|\left(\mathbf{G}_{\Gamma}^{(n-s,n)}\right)^{-1}\Big|\Big|_F^2=O(4^{s-1}(n-s)sn^{1+4(s-1)}).\end{aligned}$$ Finally, combining (\[eq:GT2\]) and (\[eq:invG\]), we conclude that $$\kappa^{max}_F(\mathbf{G}^{(n-s,n)}(\boldsymbol{\rho}^{(n)})) = O\left( (n-s)\sqrt{ns(n-s)}\left({2}n^2\right)^{s-1}\right).$$

Proof of Claim \[cl:chebcoef\] {#app:proof}
==============================

Let $\alpha=2m_2-1, \gamma=\alpha (2m_1-1)$. $p_{\mathbf{A}}(x)$ in (\[eq:genpolys\]) can be written as $$\begin{aligned}
  p_\mathbf{A}(x)&=\sum_{i=0}^{m_1-1}\sum_{j=0}^{m_2-1}\mathbf{A}_{i,j} T^{'}_{m_2-1-j+i\alpha}(x) \notag\\
  &=\notag  \sum_{j=0}^{m_2-2}\mathbf{A}_{0,j} T_{m_2-1-j}(x)+ 1/2~ \mathbf{A}_{0,m_2-1} T_{0}(x)+ \sum_{i=1}^{m_1-1}\sum_{j=0}^{m_2-1}\mathbf{A}_{i,j} T_{m_2-1-j+i\alpha}(x)\end{aligned}$$ Similarly, $p_{\mathbf{B}}(x)$ in (\[eq:genpolys\]) can be written as $$\begin{aligned}
    p_\mathbf{B}(x)&=\sum_{k=0}^{m_2-1}\sum_{l=0}^{m_3-1} \mathbf{B}_{k,l} T^{'}_{k+l\gamma}(x)\notag\\
    &= 1/2~\mathbf{B}_{0,0} T_0(x) + \sum_{k=1}^{m_2-1} \mathbf{B}_{k,0} T_{k}(x)+\sum_{k=0}^{m_2-1}\sum_{l=1}^{m_3-1} \mathbf{B}_{k,l} T_{k+l\gamma}(x)\end{aligned}$$ Now, the product $p_{\mathbf{A}}(x)p_{\mathbf{B}}(x)$ can be written as $$\begin{aligned}
    p_\mathbf{A}&(x)p_\mathbf{B}(x)=\frac{1}{2}\big(p_1(x)+p_2(x)\big)\end{aligned}$$ where, $$\begin{aligned}
p_1(x)=& \sum_{j=0}^{m_2-2}\mathbf{A}_{0,j} \mathbf{B}_{0,0}T_{m_2-1-j}(x)+ \frac{1}{2} \mathbf{A}_{0,m_2-1}\mathbf{B}_{0,0} T_0(x)+\sum_{i=1}^{m_1-1}\sum_{j=0}^{m_2-1}\mathbf{A}_{i,j}\mathbf{B}_{0,0} T_{m_2-1-j+i\alpha}(x)\notag\\
    &+  \sum_{j=0}^{m_2-2}\sum_{k=1}^{m_2-1}\mathbf{A}_{0,j}  \mathbf{B}_{k,0} T_{m_2-1-j+k}(x)+ \sum_{k=1}^{m_2-1}\mathbf{A}_{0,m_2-1}  \mathbf{B}_{k,0} T_{k}(x)\notag\\
    &+\sum_{i=1}^{m_1-1}\sum_{j=0}^{m_2-1}\sum_{k=1}^{m_2-1}\mathbf{A}_{i,j}  \mathbf{B}_{k,0} T_{m_2-1-j+i\alpha+k}(x)+\sum_{j=0}^{m_2-2}\sum_{k=0}^{m_2-1}\sum_{l=1}^{m_3-1}\mathbf{A}_{0,j} \mathbf{B}_{k,l} T_{m_2-1-j+k+l\gamma}(x)\notag\\
    &+\ \sum_{k=0}^{m_2-1}\sum_{l=1}^{m_3-1}\mathbf{A}_{0,m_2-1} \mathbf{B}_{k,l} T_{k+l\gamma}(x)+\sum_{i=1}^{m_1-1}\sum_{j=0}^{m_2-1}\sum_{k=0}^{m_2-1}\sum_{l=1}^{m_3-1}\mathbf{A}_{i,j}  \mathbf{B}_{k,l} T_{m_2-1-j+i\alpha+k+l\gamma}(x)\notag\\\end{aligned}$$ and, $$\begin{aligned}
\label{eq:p2}
&p_2(x)=     \sum_{j=0}^{m_2-2}\sum_{k=1}^{m_2-1}\mathbf{A}_{0,j}\mathbf{B}_{k,0} T_{|m_2-1-j-k|}(x) + \sum_{i=1}^{m_1-1}\sum_{j=0}^{m_2-1}\sum_{k=1}^{m_2-1}\mathbf{A}_{i,j}\mathbf{B}_{k,0}  T_{|m_2-1-j+i\alpha-k|}(x)\notag\\
    &+ \sum_{j=0}^{m_2-2}\sum_{k=0}^{m_2-1}\sum_{l=1}^{m_3-1}\mathbf{A}_{0,j}\mathbf{B}_{k,l} T_{|m_2-1-j-k-l\gamma|}(x) + \sum_{i=1}^{m_1-1}\sum_{j=0}^{m_2-1}\sum_{k=0}^{m_2-1}\sum_{l=1}^{m_3-1} \mathbf{A}_{i,j}\mathbf{B}_{k,l} T_{|m_2-1-j+i\alpha-k-l\gamma|}(x).\end{aligned}$$ Now, in order to prove the claim, it suffices to prove the following two statements:

1.  For any $i\in\{0, \cdots, m_1-1\}, l\in \{0, \cdots, m_3-1\}$, $\mathbf{C}_{i,l}$ is the matrix coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_1(x)$.

2.  For any $i\in\{0, \cdots, m_1-1\}, l\in \{0, \cdots, m_3-1\}$, the matrix coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_{2}(x)$ is $\mathbf{0}_{{N_1}/{m_1}\times{N_3}/{m_3}}$, where $\mathbf{0}_{{N_1}/{m_1}\times{N_3}/{m_3}}$ is the ${N_1}/{m_1}\times{N_3}/{m_3}$ all zeros matrix.

In the following, we prove that statement 1) is true. In order to find the coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_1(x)$, we find the set $\mathcal{S}_1=\{(i',j',k',l'): m_2-1-j'+i'\alpha+k'+l'\gamma=m_2-1+i\alpha+l\gamma\}$. Rewriting $m_2-1-j'+i'\alpha+k'+l'\gamma=m_2-1+i\alpha+l\gamma$, we have $$\begin{aligned}
\label{eqn:coef}
    (k'-j')+(i'-i) \alpha+(l'-l)\gamma=0.\end{aligned}$$ (\[eqn:coef\]) implies that $l'=l$. Suppose $l'\neq l$, this means that $(k'-j')+(i'-i)\alpha= c \gamma$ for some integer $c$. However, this is a contradiction since $|(k'-j')+(i'-i)\alpha| < \gamma$, for any $i,i',j',k'$. Now, (\[eqn:coef\]) can be written as $$\begin{aligned}
\label{eqn:coef2}
    (k'-j')+(i'-i) \alpha=0.\end{aligned}$$ Again, (\[eqn:coef2\]) implies $i'=i$. Suppose $i'\neq i$, this means $k'-j'=c\alpha$, for some integer $c$. However, this is a contradiction since $|k'-j'|<\alpha$. Now, since $i'=i$, (\[eqn:coef2\]) implies $j'=k'$. Thus, $\mathcal{S}_1=\{(i,j',j',k): j'\in\{0, \cdots, m_2-1\}\}$. That is, for any $i\in\{0, \cdots, m_1-1\}, j\in \{0, \cdots, m_3-1\}$, the matrix coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_1(x)$ is $\sum_{j'=0}^{m_2-1} \mathbf{A}_{i,j'}\mathbf{B}_{j',l} =\mathbf{C}_{i,l}$.

Now, it remains to prove statement 2). That is, for any $i\in\{0, \cdots, m_1-1\}, l\in \{0, \cdots, m_3-1\}$, the matrix coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_{2}(x)$ is $\mathbf{0}_{{N_1}/{m_1}\times{N_3}/{m_3}}$. In order to find the coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_2(x)$, we find the sets $\mathcal{S}_2^{(1)}=\{(i',j',k',l'): m_2-1-j'+i'\alpha-k'-l'\gamma=m_2-1+i\alpha+l\gamma\}$, and $\mathcal{S}_2^{(2)}=\{(i',j',k',l'): -m_2+1+j'-i'\alpha+k'+l'\gamma=m_2-1+i\alpha+l\gamma\}$.

First, for the set $\mathcal{S}_2^{(1)}$, rewriting $m_2-1-j'+i'\alpha-k'-l'\gamma=m_2-1+i\alpha+l\gamma$, we get $$\begin{aligned}
\label{eq:coef3}
(-j'-k')+(i'-i)\alpha+(l+l')\gamma=0.    \end{aligned}$$ From (\[eq:coef3\]), we conclude that $l+l'=0$. Otherwise, $(-j'-k')+(i'-i)\alpha=c \gamma$, for some integer $c$, a contradiction since $|(-j'-k')+(i'-i)\alpha| < \gamma$. Since $l+l'=0$ and both $l,l'$ are non-negative, we conclude that $l'=l=0$. Moreover, now (\[eq:coef3\]) reduces to $$\begin{aligned}
\label{eq:coef4}
(-j'-k')+(i'-i)\alpha=0.    \end{aligned}$$ Again, since $|-j'-k'| < \alpha$, we conclude that $i'=i$, which implies that $j'+k'=0$. Since $j'+k'=0$ and both $j',k'$ are non-negative, we conclude that $j'=k'=0$. Thus, $\mathcal{S}_2^{(1)}=\{(i,0,0,0)\}$. Now, noticing from (\[eq:p2\]) that $\mathbf{A}_{i,0}\mathbf{B}_{0,0}$ does not contribute to any term in $p_2(x)$, we conclude that the matrix coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_{2}(x)$ is only due to the set $\mathcal{S}_2^{(2)}$. Recall that $\mathcal{S}_2^{(2)}=\{(i',j',k',l'): -m_2+1+j'-i'\alpha+k'+l'\gamma=m_2-1+i\alpha+l\gamma\}$, we rewrite $ -m_2+1+j'-i'\alpha+k'+l'\gamma=m_2-1+i\alpha+l\gamma$ as $$\begin{aligned}
\label{eq:coef5}
    (j'+k'-2m_2+2)-(i'+i)\alpha+(l'-l)\gamma=0\end{aligned}$$ From (\[eq:coef5\]), we conclude that $l=l'$. Otherwise, $(j'+k'-2m_2+2)-(i'+i)\alpha=c \gamma$, for some integer $c$, a contradiction since $|(j'+k'-2m_2+2)-(i'+i)\alpha| < \gamma$. Moreover, now (\[eq:coef5\]) reduces to $$\begin{aligned}
\label{eq:coef6}
(j'+k'-2m_2+2)+(i'+i)\alpha=0.    \end{aligned}$$ Again, since $|j'+k'-2m_2+2| < \alpha$, we conclude that $i'+i=0$. Since $i'+i=0$ and both $i,i'$ are non-negative, we conclude that $i'=i=0$, which implies that $j'+k'=2m_2-2$. Since $j'+k'=2m_2-2$ and both $j',k'\leq m_2-1$, we conclude that $j'=k'=m_2-1$. Thus, $\mathcal{S}_2^{(2)}=\{(0,m_2-1,m_2-1,l)\}$. Now, noticing from (\[eq:p2\]) that $\mathbf{A}_{0,m_2-1}\mathbf{B}_{m_2-1,l}$ does not contribute to any term in $p_2(x)$, we conclude that the matrix coefficient of $T_{m_2-1+i\alpha+l\gamma}$ in $p_{2}(x)$ is $\mathbf{0}_{{N_1}/{m_1}\times{N_3}/{m_3}}$.[$\Box$]{}

Upper Bound on the Condition Number of Gaussian Matrices {#app:iidG}
========================================================

We first introduce the following theorem from [@azais2004].

\[thm:azais\] Let $\mathbf{A}$ be an $m\times m$ matrix, $m\geq 3$, and let the entries of $\mathbf{A}$ be independent and identically distributed standard Gaussian random variables. Then, for all $\alpha > 1$, $$\operatorname{Pr}(\kappa_2(\mathbf{A})>m\alpha)<\frac{5.6}{\alpha},$$ where $\kappa_2(\mathbf{A})$ is the condition number of $\mathbf{A}$ with respect to the matrix norm induced by $\ell_2$.

As a consequence, in the following, we extend the result in Theorem \[thm:azais\] to bound the condition number of every $m\times m$ sub-matrix of a random $m\times P$ matrix with $i.i.d$ standard Gaussian entries, $P \geq m$.

For any subset $\mathcal{S} \subseteq \{1,2,\ldots,P\}$, let $\mathbf{H}_{\mathcal{S}}$ denote the $|\mathcal{S}| \times m$ sub-matrix of $\mathbf{H}$ containing the columns $\mathbf{H}$ corresponding to $\mathcal{S},$ and let $s=P-m$. Then we have $$\begin{aligned}
\operatorname{Pr}\Big(\kappa^{max}_2(\mathbf{H}) > mP^{2s}\Big) &=  \operatorname{Pr}\Big(\underset{\substack{\mathcal{S}'\subset [P], |\mathcal{S}'|=m}}{\bigcup}\big(\kappa_2(\mathbf{H}_{\mathcal{S}'})> m P^{2s}\big)\Big)\\ &\stackrel{(1)}{\leq}\underset{\substack{\mathcal{S}'\subset [P], |\mathcal{S}'|=m}}{\sum}\operatorname{Pr}\big(\kappa_2(\mathbf{H}_{\mathcal{S}'})> m P^{2s}\big)   \notag \\
&\stackrel{}{=}\left(\begin{array}{cc}
     P\\s \end{array}\right)  \operatorname{Pr}\big(\kappa_2(\mathbf{H}_{\mathcal{S}'})> m P^{2s}\big), \text{~for any } \mathcal{S}'\subset [P] \text{~such that  } |\mathcal{S}'|= m\notag\\
&\stackrel{(2)}{<} P^s  \frac{5.6}{P^{2s}}    \notag\\
&=\frac{5.6}{P^s},\notag\end{aligned}$$ where $(1)$ follows from the union bound, and $(2)$ follows from the fact that $\left(\hspace{-2mm}\begin{array}{cc}
     P\\s \end{array}\hspace{-2mm}\right) \leq P^s$ and Theorem \[thm:azais\].

[^1]: M. Fahim and V. Cadambe are with the Department of Electrical Engineering, Pennsylvania State University, University Park, PA 16802.

[^2]: This work will be presented in part at the IEEE International Symposium on Information Theory (ISIT), July 2019.

[^3]: For example, [@seth2018bigdata], reports that “In our experiments we observed large floating point errors when inverting high degree Vandermonde matrices for polynomial interpolation”.

[^4]: The master and fusion nodes are logical entities; in practice, they may be the same node, or may be emulated in a decentralized manner by the computation nodes.

[^5]: We note that the numerical error depends not only on the condition number of the matrix, but also the algorithm used for solving the linear system. However, we are not aware of any approach that can accurately solve, say, a $150 \times 150$ linear system with a Vandermonde matrix (See e.g., [@demmel2005accurate; @ramamoorthy2019universally])
