# Distributed Structured Matrix Multiplication

Derya Malak  Co-funded by the European Union (ERC, SENSIBILITÉ, 101077361).
Views and opinions expressed are however those of the author only and do not
necessarily reflect those of the European Union or the European Research
Council. Neither the European Union nor the granting authority can be held
responsible for them. This research was partially supported by a Huawei
France-funded Chair towards Future Wireless Networks, and supported by the
program “PEPR Networks of the Future" of France 2030. Communication Systems
Department, EURECOM, Biot Sophia Antipolis, FRANCE
<EMAIL_ADDRESS>
###### Abstract

We devise achievable encoding schemes for distributed source compression for
computing inner products, symmetric matrix products, and more generally,
square matrix products, which are a class of nonlinear transformations. To
that end, our approach relies on devising nonlinear mappings of distributed
sources, which are then followed by the structured linear encoding scheme,
introduced by Körner and Marton. For different computation scenarios, we
contrast our findings on the achievable sum rate with the state of the art to
demonstrate the possible savings in compression rate. When the sources have
special correlation structures, it is possible to achieve unbounded gains, as
demonstrated by the analysis and numerical simulations.

###### Index Terms:

Distributed computation, inner product, structured codes, matrix-vector
multiplication, matrix multiplication.

##  I Introduction

The inner product operation between two vectors captures the similarity
between vectors and allows us to describe the lengths, angles, projections,
vector norms, matrix norms induced by vector norms, orthogonality of vectors,
polynomials, and a variety of other functions as well [1]. Inner products are
widely used in geometry and trigonometry using linear algebra and in
applications spanning physics, engineering, and mathematics, e.g., to
determine the convolution of functions [1], and the Fourier transform
approximations, machine learning [2] and pattern recognition [3], e.g., the
linear regression and the least squares models [1], and quantum computing,
e.g., to describe the overlap between the two quantum states [4].

With the advent of massive parallelization techniques and coded computation
frameworks, modern distributed computing systems, e.g., MapReduce, Hadoop, and
Spark, have been devised to implement the computationally intensive task of
distributed matrix multiplication with low communication and computation cost
[5]. To that end, novel coded matrix-multiplication constructions, e.g.,
polynomial codes [6], [7], gradient coding [5], and Lagrange coded computing
[8, 9], have been designed to mitigate the costs, faulty nodes, and
stragglers. The inner product computation serves as the primary building block
of such settings.

In this paper, we devise structured encoding schemes for distributed computing
of inner products and symmetric and, more generally, square matrices via
distributed matrix products. Our main contributions are summarized as follows:

  * •

We devise a distributed encoding scheme that performs structured coding on
nonlinear mappings of two distributed sources ${\bf A}$ and ${\bf B}$ to
compute their inner product.

  * •

We showcase the conditions for which the sum rate achieved by this structured
coding is strictly less than the sum rate of distributed unstructured encoding
of the sources ${\bf A}$ and ${\bf B}$ [10]. Here, the performance criterion
is that both ${\bf A}$ and ${\bf B}$ cannot be decoded by the receiver.

  * •

We derive achievable rate regions for distributed computation of symmetric and
square matrices, a class of nonlinear transformations, with a vanishing
probability of error, and determine example scenarios — detailed in
Corollaries 1-2 — with special correlation structures to guarantee, via the
structured coding scheme of Körner and Marton [11], significant savings over
[10].

  * •

We contrast the achievable rates with the existing approaches (e.g., [10, 12,
13]) via numerical examples.

Connections to the state of the art. Slepian and Wolf have provided an
unstructured coding technique for the asymptotic lossless compression of
distributed source variables $X_{1}$ and $X_{2}$ at the minimum rate needed,
i.e., $H(X_{1},X_{2})$ [10]. Han-Kobayashi [14] have provided a
characterization to determine whether computing a general bivariate function
$f(X_{1},X_{2})$ of two random sequences $\\{X_{1i}\\}$ and $\\{X_{2i}\\}$
from two correlated memoryless sources requires a smaller rate than
$H(X_{1},X_{2})$. For distributed coding of a finite alphabet source $X$ with
side information $Y$, Orlitsky and Roche have devised an unstructured coding
scheme to achieve the minimum rate at which source $X$ has to compress for
distributed computing of $f(X,Y)$ with vanishing error [15], exploiting
Körner’s characteristic graph $G_{X}$ and its entropy [16]. This scheme is
equivalent to performing Slepian-Wolf encoding on the colors of the
sufficiently large OR powers of $G_{X}$ given $Y$ [17, 18, 19].

Körner and Marton have devised a structured encoding strategy that minimizes
the sum rate for distributed computing the modulo-two sum of doubly symmetric
binary source (DSBS) sequences with a low probability of error [11]. Ahlswede
and Han have tightened the rate region for general binary sources [20] that
embed the regions of [11] and [10]. The rate region for this problem has been
extended to a larger class of source distributions [21], and for
reconstructing the modulo-$q$ sum of the two sources in a $q$-ary prime finite
field $\mathbb{F}_{q}$ at a sum rate of $2H(X_{1}\oplus_{q}X_{2})$ [14]. For
computing a nonlinear function, the embedding of the function in a
sufficiently large prime $\mathbb{F}_{q}$ [12], [13], and finding an injective
mapping between the function and $X_{1}\oplus_{q}X_{2}$, known as _structured
binning_ , may provide savings over [10].

In this paper, leveraging these fundamental principles, we demonstrate further
savings in compression for computing inner products as well as matrix products
through devising nonlinear mappings of the sources followed by the linear
encoding scheme in [11], while requiring a smaller prime field size versus
[12], [13].

Notation. We denote by $H(X)$ the Shannon entropy of a discrete random
variable $X$, which is drawn from probability mass function (PMF) $P_{X}$.
Similarly, $H(\cdot,\cdot)$ and $H(\cdot\,|\,\cdot)$ denote the joint and
conditional entropies, given a joint PMF $P_{X_{1},X_{2}}$. Let $h(\epsilon)$
denote the binary entropy function for Bernoulli distributed $X$ with
parameter $\epsilon\in[0,1]$, i.e., $X\sim{\rm Bern}(\epsilon)$.

We denote by $X^{n}$ the length $n$ realization of $X$. The boldface notation
${\bf X}$ denotes a random matrix with elements in $\mathbb{F}_{q}$ and ${\bf
X}^{\intercal}$ is its transpose. ${\bm{1}}_{m}$ and ${\bm{0}}_{m}$ are length
$m$ all-ones and all-zeros vectors, respectively. $\mathbb{P}(A)$ is the
probability of an event $A$, and $1_{x\in A}$ is the indicator function which
takes the value $1$ if $x\in A$, and $0$ otherwise. The inner product of two
vectors in the vector space $V$ over a field $F$ is a scalar, which is a map
$\langle\cdot,\cdot\rangle:V\times V\to F$.

##  II System Model and Main Results

We consider a distributed coding scenario with two sources and a receiver. The
two sources are assigned matrix variables ${\bf A}\in\mathbb{F}_{q}^{m\times
l}$ and ${\bf B}\in\mathbb{F}_{q}^{m\times l}$, respectively, that model two
statistically dependent independent and identically distributed (i.i.d.)
finite alphabet source matrix sequences with entries from a field whose
characteristic is $q\geq 2$. The objective of the receiver is to compute a
function $f({\bf A},{\bf B})$, which is the matrix product of ${\bf A}$ and
${\bf B}$, i.e., ${\bf A}^{\intercal}{\bf B}:\mathbb{F}_{q}^{m\times
l}\times\mathbb{F}_{q}^{m\times l}\to\mathbb{F}_{q}^{l\times l}$. To that end,
we devise a distributed encoding scheme that performs structured coding on
nonlinear source mappings.

We demonstrate different regimes where the sum rate achieved by _structured
coding_ is strictly less than $H({\bf A},{\bf B})$ [10] so that the receiver
can recover ${\bf A}^{\intercal}{\bf B}$, while it is desired that sources
${\bf A}$ and ${\bf B}$ cannot be decoded by the receiver, implying a security
constraint. Exploiting [11], we first study the special scenario of
distributed computing of inner products, where the key idea is to implement
linear structured coding of [11] on vector-wise embeddings of the sources
(Propositions 1 and 2), and propose a hybrid scheme that combines [11] with
the _unstructured coding_ technique of [15]. We contrast the achievable rate
savings over unstructured coding [10], and over linear embedding of ${\bf
A}^{\intercal}{\bf B}$ in a prime finite field [12] and [13]. We generalize
our scheme for computing symmetric matrices via distributed multiplication of
the sources without being able to recover the sources in their entirety
(Proposition 3) and provide a necessary condition for not being able to
recover ${\bf A}$ and ${\bf B}$ (Proposition 4). Finally, we consider the
distributed computation of square matrices (Proposition 5 and Corollary 2). We
showcase the compression performance of our schemes versus the existing
results, demonstrating the achievable savings.

We next detail our main achievability results.

###  II-A Distributed Computation of Inner Products of Sources

The distributed sources hold the even-length vector variables ${\bf
A}=\begin{bmatrix}a_{1}&a_{2}&\ldots&a_{m}\end{bmatrix}^{\intercal}$ and ${\bf
B}=\begin{bmatrix}b_{1}&b_{2}&\ldots&b_{m}\end{bmatrix}^{\intercal}$,
respectively, with entries chosen from a prime field $\mathbb{F}_{q}$. The
receiver aims to compute the inner product $f({\bf A},{\bf B})=\langle\,{\bf
A},{\bf B}\,\rangle$.

We next present an achievable coding scheme for distributed computation of
$\langle\,{\bf A},{\bf B}\,\rangle=\sum\nolimits_{i=1}^{m}a_{i}b_{i}$. The
coding scheme relies on devising nonlinear mappings from each source and using
the linear encoding scheme of Körner-Marton in [11].

######  Proposition 1.

(Distributed inner product computation.) Given two sequences of random vectors
${\bf A}=\begin{bmatrix}{\bf A}_{1}^{\intercal}&{\bf
A}_{2}^{\intercal}\end{bmatrix}^{\intercal}\in\mathbb{F}_{q}^{m\times 1}$ and
${\bf B}=\begin{bmatrix}{\bf B}_{1}^{\intercal}&{\bf
B}_{2}^{\intercal}\end{bmatrix}^{\intercal}\in\mathbb{F}_{q}^{m\times 1}$ of
even length $m$, generated by two correlated memoryless $q$-ary sources, where
${\bf A}_{1},{\bf A}_{2},{\bf B}_{1},{\bf B}_{2}\in\mathbb{F}_{q}^{m/2\times
1}$, the following sum rate is achievable by the separate encoding of the
sources for the receiver to recover $\langle\,{\bf A},{\bf B}\,\rangle$ with a
small probability of error:

$\displaystyle R_{\rm KM}^{\Sigma}=2H({\bf U},\,{\bf V},\,W)\ ,$ (1)

where ${\bf U},{\bf V}\in\mathbb{F}_{q}^{m/2\times 1}$ are vector variables,
and $W\in\mathbb{F}_{q}$ is a random variable, and they satisfy the following
relations:

$\displaystyle{\bf U}$ $\displaystyle={\bf A}_{2}\oplus_{q}{\bf B}_{1}\ ,$
$\displaystyle{\bf V}$ $\displaystyle={\bf A}_{1}\oplus_{q}{\bf B}_{2}\ ,$
$\displaystyle W$ $\displaystyle={\bf A}_{2}^{\intercal}{\bf
A}_{1}\oplus_{q}{\bf B}_{1}^{\intercal}{\bf B}_{2}\ .$ (2)

###### Proof.

See Appendix -D. ∎

Figure 1: Gain, $\eta$ (cf. (3) in Corollary 1. The flat (yellow) line marks
$\eta=1$.

In the scheme of Proposition 1, the receiver may not recover ${\bf A}$ and
${\bf B}$ in their entirety, from the modulo-$q$ additions ${\bf U}$, ${\bf
V}$, and $W$, enabling secure distributed inner product computation.

We next contrast the sum rate for computing $\langle\,{\bf A},{\bf
B}\,\rangle$ given by Proposition 1 with the sum rate of Slepian-Wolf in [10]
to demonstrate the gap between $R_{\rm KM}^{\Sigma}$ and $R_{\rm
SW}^{\Sigma}$. The following result indicates that the rate needed to compute
$\langle\,{\bf A},{\bf B}\,\rangle$ may be substantially less than $H({\bf
A},{\bf B})$ for binary-valued sources.

######  Corollary 1.

Consider two sequences of correlated source vectors ${\bf
A}\in\mathbb{F}_{2}^{m\times 1}$ and ${\bf B}\in\mathbb{F}_{2}^{m\times 1}$,
with entries $a_{i}$ and $b_{i}$ that are i.i.d. across $i=1,\dots,m$ each,
where $a_{i}\sim{\rm Bern}(1/2)$ and $b_{i}\sim{\rm Bern}(1/2)$ are correlated
for a given $i=1,\dots,m$. Assume that ${\bf U}$ and ${\bf V}$, as defined in
(1), have entries $u_{i},v_{i}\sim{\rm Bern}(p)$, that are i.i.d. across
$i=1,\dots,m/2$.

For this setting, the gain of the encoding technique in Proposition 1 with a
sum rate $R_{\rm KM}^{\Sigma}$ given in (1) over the sum rate $R_{\rm
SW}^{\Sigma}$ for lossless compression of the sources is

$\displaystyle\eta=\frac{R_{\rm SW}^{\Sigma}}{R_{\rm
KM}^{\Sigma}}=\frac{m(1+h(p))}{2mh(p)+2(1-(1-p)^{m})}\ .$ (3)

###### Proof.

The proof follows from evaluating (1) and contrasting it with $R_{\rm
SW}^{\Sigma}=H({\bf A},{\bf B})=m(1+h(p))$. For details, we refer the reader
to Appendix -E. ∎

Note from Corollary 1 that when $\eta>1$, the receiver can compute
$\langle\,{\bf A},{\bf B}\,\rangle$ without recovering $({\bf A},{\bf B})$. It
holds that

$\displaystyle\lim\limits_{p\to 0}\eta=\infty\ ,\quad\lim\limits_{p\to
1}\eta=\frac{m}{2}\ ,\quad\lim\limits_{m\to\infty}\eta=\frac{1+h(p)}{2h(p)}\
,$ (4)

where the limit $\lim\limits_{m\to\infty}\eta$ is the same as the gain for the
DSBS model studied in [11], which tends to infinity as $p\to\\{0,1\\}$.

Figure 2: Rate comparisons for various source PMFs. (Left) Corollary 1 for
$m=2$. (Middle) $m=1$, and $a,b$ are DSBSs. (Right) $m=2$, and
$\\{a_{i},b_{i}\\}_{i=1}^{2}$ are DSBSs.

We illustrate the gain $\eta$ as a function of $(m,p)$ in Figure 1, indicating
that $R_{\rm KM}^{\Sigma}$ may be substantially less than the joint entropy of
the sources for this special class of source PMFs.

It has been shown in [12], [13] that via embedding the nonlinear function
$D_{k}=a_{k}b_{k}$, where $a_{k},b_{k}\in\mathbb{F}_{q}$, in a sufficiently
large prime $\mathbb{F}_{q}$, the decoder can reconstruct
$\tilde{D}^{n}=\\{a_{k}\oplus_{q}b_{k}\\}_{k=1}^{n}$, and hence, compute
$D^{n}=\\{a_{k}b_{k}\\}_{k=1}^{n}$ with high probability. For instance, if
$a_{k},b_{k}\in\mathbb{F}_{2}$, we can reconstruct $D^{n}$ from
$\tilde{D}^{n}=\\{a_{k}\oplus_{3}b_{k}\\}_{k=1}^{n}$ using a sum rate of
$R_{\rm S}^{\Sigma}=2H({\bf A}\oplus_{3}{\bf B})$.

Motivated by the notion of embedding in [12] and [13], we next devise an
achievability scheme for computing $\langle\,{\bf A},{\bf B}\,\rangle$, where
the key idea is to compress the vector-wise embeddings of the sources vectors
${\bf A}$ and ${\bf B}$ via employing the linear structured encoding scheme of
[11], in contrast to entry-wise embeddings that require a sum rate of $R_{\rm
S}^{\Sigma}$ (cf. [12] and [13]).

######  Proposition 2.

(Vector-wise embeddings followed by linear encoding for distributed
computation of $\langle\,{\bf A},{\bf B}\,\rangle$.) Given two sequences of
vectors ${\bf A}\in\mathbb{F}_{q}^{m\times 1}$ and ${\bf
B}\in\mathbb{F}_{q}^{m\times 1}$, generated by two correlated memoryless
sources, with entries from a field $\mathbb{F}_{q}$ with $q>2$, restricting
the source mappings to be linear, the following sum rate is achievable via the
Körner-Marton’s scheme to recover $\langle\,{\bf A},{\bf B}\,\rangle$ at the
receiver with a small probability of error:

$\displaystyle R_{\rm
SV}^{\Sigma}=2H\Big{(}\Big{\\{}a_{i}\oplus_{r}b_{i}\Big{\\}}_{i=1}^{m},\,{\bigoplus_{i=1}^{m}}{{}_{{}_{{\operatorname*{{\small
q}}}}}\,a_{i}^{2}\oplus_{q}b_{i}^{2}}\Big{)}\ ,$ (5)

where $r=2(q-1)m$ for $m$ even, and $r=2(q-1)m+1$ for $m$ odd, respectively,
and ${\bigoplus_{{}_{{\operatorname*{{\small q}}}}}}$ denotes a modulo-$q$
addition.

###### Proof.

The proof follows from noting that the receiver, upon receiving
$\big{\\{}a_{i}\oplus_{r}b_{i}\big{\\}}_{i=1}^{m}$ and
${\bigoplus\limits_{i=1}^{m}}{{}_{{}_{{\operatorname*{{\small
q}}}}}\,a_{i}^{2}\oplus_{q}b_{i}^{2}}$, can reconstruct

$\displaystyle
2c=qk+\Big{(}\sum\limits_{i=1}^{m}(a_{i}\oplus_{r}b_{i})^{2}-\,{\bigoplus\limits_{i=1}^{m}}{{}_{{}_{{\operatorname*{{\small
q}}}}}\,(a_{i}^{2}\oplus_{q}b_{i}^{2})\Big{)}}\mod q\ ,$

where there is a unique $k\in\mathbb{F}_{q}$ for which $c=\langle\,{\bf
A},{\bf B}\,\rangle\in\mathbb{F}_{q}$.

For details, we refer the reader to Appendix -F. ∎

We next describe a hybrid encoding scheme. Note that ${\bf A}^{\intercal}{\bf
B}$ can be rewritten as $g({\bf A},{\bf Y})={\bf A}^{\intercal}({\bf Y}-{\bf
A})\,\mod q$ if ${\bf Y}={\bf A}\oplus_{q}{\bf B}$ is known. Exploiting
_Körner’s characteristic graphs_ [16] to enable nonlinear encoding of ${\bf
A}$111For a detailed description of characteristic graphs and their entropies,
we refer the reader to [17, 18, 19]., the minimum compression rate of ${\bf
A}$ for computing $g({\bf A},{\bf Y})$ given side information ${\bf Y}$ equals
the _conditional graph entropy_ of ${\bf A}$ given ${\bf Y}$, denoted by
$H_{G_{\bf A}}({\bf A}\,|\,{\bf Y})$, as introduced by Orlitsky and Roche
[15]. Via concatenating the structured coding scheme of Körner and Marton [11]
to first compute ${\bf Y}$ and then the unstructured coding model of Orlitsky
and Roche [15] to next determine $g({\bf A},{\bf Y})$, it is possible to
achieve a sum rate $R_{\rm KM-OR}^{\Sigma}=2H({\bf Y})+H_{G_{\bf A}}({\bf
A}\,|\,{\bf Y})$. When ${\bf Y}={\bf B}$, the required rate is $H({\bf
B})+H_{G_{\bf A}}({\bf A}\,|\,{\bf B})$, which is smaller than $R_{\rm
SW}^{\Sigma}$ because $H_{G_{\bf A}}({\bf A}\,|\,{\bf B})\leq H({\bf
A}\,|\,{\bf B})$ [15].

In Figure 2, we contrast the sum rate performance of Propositions 1-2 for
distributed computing of $\langle\,{\bf A},{\bf B}\,\rangle$ of vectors ${\bf
A},{\bf B}\in\mathbb{F}_{2}^{m\times 1}$ with a small probability of error,
with existing unstructured and structured coding schemes. In Figure 2-(Left),
we use the PMF in Corollary 1 where $m=2$, i.e., the pairs $(a_{1},b_{2})$ and
$(a_{2},b_{1})$ represent DSBSs, each with a crossover probability $p$. Note
that the sum rates $R_{\rm S}^{\Sigma}$ and $R_{\rm KM-OR}^{\Sigma}$ perform
poorly versus $R_{\rm SV}^{\Sigma}$, and are not indicated. The sum rate
$R_{\rm KM}^{\Sigma}$ converges to $H(\langle\,{\bf A},{\bf B}\,\rangle)$ at
low $p$.

In Figure 2-(Middle), we use $m=1$, and the pair $(a,b)$ is a DSBS with a
crossover probability $p$. At low $p$, $R_{\rm KM-OR}^{\Sigma}$ and $R_{\rm
SW}^{\Sigma}$ converge to $H(\langle\,a,b\,\rangle)$ whereas $R_{\rm
S}^{\Sigma}$ performs poorly. For large $p$, structured coding yields low
rates ($R_{\rm S}^{\Sigma}$ and $R_{\rm KM-OR}^{\Sigma}$).

In Figure 2-(Right), we use $m=2$, and the pairs $(a_{1},b_{1})$ and
$(a_{2},b_{2})$ are DSBSs, each with a crossover probability $p$. The
performance of $R_{\rm KM}^{\Sigma}$ is worse than $R_{\rm SW}^{\Sigma}$ and
not shown. Similarly, $R_{\rm SV}^{\Sigma}$ is higher than $R_{\rm KM-
OR}^{\Sigma}$, and is not indicated. For any given $p$ value, $R_{\rm
S}^{\Sigma}$ is always smaller than $R_{\rm SW}^{\Sigma}$, and $R_{\rm KM-
OR}^{\Sigma}$ approaches $H(\langle\,{\bf A},{\bf B}\,\rangle)$ for small and
large $p$.

###  II-B Distributed Computation of Symmetric Matrices

We next consider a generalization of Proposition 1 for distributed computing
of inner products to distributed computing of a square symmetric matrix ${\bf
D}=(D_{ij})\in\mathbb{F}_{q}^{l\times l}$, given by the product ${\bf D}={\bf
A}^{\intercal}{\bf B}$, where ${\bf A},{\bf B}\in\mathbb{F}_{q}^{m\times l}$,
for $q>2$ and $l\geq 1$, and symmetry implies that $D_{ji}=D_{ij}$ for every
$i,j$.

######  Proposition 3.

(Computing symmetric matrices via distributed multiplication.) Given two
sequences of random matrices ${\bf A}=\begin{bmatrix}{\bf
A}_{1}^{\intercal}&{\bf
A}_{2}^{\intercal}\end{bmatrix}^{\intercal}\in\mathbb{F}_{q}^{m\times l}$ and
${\bf B}=\begin{bmatrix}{\bf B}_{1}^{\intercal}&{\bf
B}_{2}^{\intercal}\end{bmatrix}^{\intercal}\in\mathbb{F}_{q}^{m\times l}$
generated by two correlated memoryless $q$-ary sources, respectively, where
${\bf A}_{1},{\bf A}_{2},{\bf B}_{1},{\bf B}_{2}\in\mathbb{F}_{q}^{m/2\times
l}$ and $q>2$, the achievable sum rate by the separate encoding of the sources
for the receiver to recover the symmetric matrix ${\bf D}={\bf
A}^{\intercal}{\bf B}\in\mathbb{F}_{q}^{l\times l}$ with vanishing error is
given as

$\displaystyle R_{\rm KM}^{\Sigma}=2H({\bf U},\,{\bf V},\,{\bf W})\ ,$ (6)

where ${\bf U}={\bf A}_{2}\oplus_{q}{\bf B}_{1}\in\mathbb{F}_{q}^{m/2\times
l}$, ${\bf V}={\bf A}_{1}\oplus_{q}{\bf B}_{2}\in\mathbb{F}_{q}^{m/2\times
l}$, and ${\bf W}={\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{q}{\bf
B}_{1}^{\intercal}{\bf B}_{2}\in\mathbb{F}_{q}^{l\times l}$ are matrix
variables.

###### Proof.

The proof is on similar lines as Appendix -D. For details, we refer the reader
to Appendix -G. ∎

In the scheme of Proposition 3, similar to Proposition 1, the receiver may not
recover ${\bf A}$ and ${\bf B}$ in their entirety, using ${\bf U}$, ${\bf V}$,
and ${\bf W}$, hence rendering secure distributed compression of the sources
feasible. To achieve information-theoretically secure distributed matrix
multiplication, our future work includes exploiting the polynomial encoding
scheme in [22] that incorporates random matrices, and the techniques in [23],
[24].

In Figure 3-(Left), we showcase the sum rates $R_{\rm KM}^{\Sigma}$ and
$R_{\rm SW}^{\Sigma}$ versus $p$ (in log scale) for distributed computing of
symmetric matrices ${\bf D}={\bf A}^{\intercal}{\bf B}$ for $q=2$ under
assumptions222For $q>2$, it is clear (from Appendix -G of Proposition 3) that
the choices of ${\bf U}$, ${\bf V}$, and ${\bf W}$ guarantee the recovery of
${\bf A}^{\intercal}{\bf B}$ without assumptions i)-ii).: i) ${\bf
A}_{1}^{\intercal}{\bf B}_{1}={\bf B}_{1}^{\intercal}{\bf A}_{1}$, i.e., ${\bf
W}={\bm{0}}_{l\times l}$, and ii) ${\bf A}_{2}^{\intercal}{\bf A}_{1}={\bf
B}_{1}^{\intercal}{\bf B}_{2}$. Because ${\bf D}$ is symmetric, i) and ii)
lead to

$\displaystyle{\bf U}^{\intercal}{\bf V}$ $\displaystyle=({\bf
A}_{2}\oplus_{2}{\bf B}_{1})^{\intercal}({\bf A}_{1}\oplus_{2}{\bf B}_{2})$
$\displaystyle\overset{(a)}{=}{\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{2}{\bf
A}_{2}^{\intercal}{\bf B}_{2}\oplus_{2}{\bf A}_{1}^{\intercal}{\bf
B}_{1}\oplus_{2}{\bf B}_{1}^{\intercal}{\bf B}_{2}$
$\displaystyle\overset{(b)}{=}{\bf A}_{1}^{\intercal}{\bf B}_{1}\oplus_{2}{\bf
A}_{2}^{\intercal}{\bf B}_{2}={\bf D}\ ,$

where $(a)$ and $(b)$ follow from assumptions i) and ii), respectively. These
assumptions ensure a rate gain of $\eta={R_{\rm SW}^{\Sigma}}/{R_{\rm
KM}^{\Sigma}}$ that grows exponentially fast, as $p$ tends to $\\{0,1\\}$.

We next state a necessary condition for successful recovery of ${\bf D}={\bf
A}^{\intercal}{\bf B}$, where ${\bf A},{\bf B}\in\mathbb{F}_{q}^{m\times l}$,
without recovering ${\bf A}$ and ${\bf B}$. This result holds true for any
symmetric ${\bf D}\in\mathbb{F}_{q}^{l\times l}$.

######  Proposition 4.

(A necessary condition for the nonrecovery of the sources in distributed
computation of ${\bf A}^{\intercal}{\bf B}$.) For distributed encoding of
${\bf A},{\bf B}\in\mathbb{F}_{q}^{m\times l}$ to compute ${\bf D}={\bf
A}^{\intercal}{\bf B}$, which is symmetric, the below condition ensures that
the sum rate in (1) is less than the achievable sum rate of [10]:

$\displaystyle H({\bf A}^{\intercal}{\bf B})+H({\bf Q}\,|\,{\bf
A}^{\intercal}{\bf B})<H({\bf A}\,|\,{\bf Q},{\bf A}^{\intercal}{\bf B})\ ,$
(7)

where ${\bf Q}=\begin{bmatrix}{\bf U}\\\ {\bf
V}\end{bmatrix}\in\mathbb{F}_{q}^{m\times 1}$, and $U$ and $V$ are defined in
(1).

###### Proof.

The inequality in (7) can be obtained using the expansions of $R_{\rm
KM}^{\Sigma}$ and $R_{\rm SW}^{\Sigma}$. For the technical steps, we refer the
reader to Appendix -H. ∎

It is clear, using the expansion (-H) of $R_{\rm KM}^{\Sigma}$ in Appendix -H,
that the sum rate $R_{\rm KM}^{\Sigma}=2H({\bf A}^{\intercal}{\bf B})$ is
achievable when $H({\bf Q}\,|\,{\bf A}^{\intercal}{\bf B})=0$ in (7). In this
case, (7) implies that $H({\bf A}^{\intercal}{\bf B})<H({\bf A}\,|\,{\bf
Q},{\bf A}^{\intercal}{\bf B})$, meaning that it is possible to recover the
inner product ${\bf A}^{\intercal}{\bf B}$ while keeping ${\bf A}$ and ${\bf
B}$ unknown to the receiver.

###  II-C Distributed Computation of Square Matrices

Recall that Proposition 3 does not capture non-symmetric matrix products. We
here consider distributed computing of a square matrix ${\bf
D}=(D_{ij})\in\mathbb{F}_{q}^{l\times l}$, given by the product ${\bf D}={\bf
A}^{\intercal}{\bf B}$, where ${\bf D}$ is not symmetric, and ${\bf A},{\bf
B}\in\mathbb{F}_{q}^{m\times l}$, for $l,\,m>1$, and $q>2$. The following
proposition gives an achievable distributed encoding scheme of the sources
${\bf A}$ and ${\bf B}$ towards computing ${\bf D}$.

######  Proposition 5.

(Computing square matrices via distributed matrix multiplication.) Given two
sequences of random matrices ${\bf A}\in\mathbb{F}_{q}^{m\times l}$ and ${\bf
B}\in\mathbb{F}_{q}^{m\times l}$ generated by two correlated memoryless
$q$-ary sources, where $q>2$, the following sum rate is achievable by the
separate encoding of the sources for the receiver to recover a general square
matrix ${\bf D}={\bf A}^{\intercal}{\bf B}\in\mathbb{F}_{q}^{l\times l}$ with
vanishing error:

$\displaystyle R_{\rm KM}^{\Sigma}=2H(\\{{\bf
A}\oplus_{q}{\bf\tilde{B}}_{j}\\}_{j=1}^{l}\,,\,\\{{\bf\tilde{A}}_{j}^{\intercal}{\bf\tilde{A}}_{j}\\}_{j=1}^{l})\
,$ (8)

where we use the shorthand notation
${\bf\tilde{A}}_{j}^{\intercal}{\bf\tilde{A}}_{j}={\bf A}^{\intercal}{\bf
A}\oplus_{q}{\bf\tilde{B}}_{j}^{\intercal}{\bf\tilde{B}}_{j}$ for
$j\in\\{1,2,\dots,l\\}$, where ${\bf\tilde{B}}_{j}={\bf
B}_{j}{\bm{1}}_{1\times l}\in\mathbb{F}_{q}^{m\times l}$ are matrix variables,
${\bm{1}}_{1\times l}$ is a length $l$ row vector of all ones, where ${\bf
B}=\begin{bmatrix}{\bf B}_{1}&{\bf B}_{2}&\ldots&{\bf B}_{l}\end{bmatrix}$
with ${\bf B}_{j}\in\mathbb{F}_{q}^{m\times 1}$ for $j\in\\{1,2,\dots,l\\}$.

###### Proof.

We refer the reader to Appendix -I. ∎

Figure 3: Rate (in log scale) versus $p$ for distributed computing of (Left)
symmetric matrices ${\bf A}^{\intercal}{\bf B}={\bf B}^{\intercal}{\bf A}$ via
distributed multiplication of matrices ${\bf A},{\bf
B}\in\mathbb{F}_{2}^{m\times m}$ for different $m$, and (Right) square
matrices via distributed matrix multiplication for different $m$ and $l=2$,
where the joint source PMF is given in Corollary 2.

To demonstrate the performance of Proposition 5, we next consider a corollary,
where $l=2$, $m>1$, and with $q=3$.

######  Corollary 2.

Given ${\bf A}=\begin{bmatrix}{\bf A}_{1}&{\bf
A}_{2}\end{bmatrix}\in\mathbb{F}_{3}^{m\times 2}$ and ${\bf
B}=\begin{bmatrix}{\bf B}_{1}&{\bf
B}_{2}\end{bmatrix}\in\mathbb{F}_{3}^{m\times 2}$ with entries
$a_{ij}\sim\big{(}\frac{1}{2}-\epsilon,2\epsilon,\frac{1}{2}-\epsilon\big{)}$
for some $\epsilon\in\big{[}0,\frac{1}{2}\big{]}$, and $b_{i1}=b_{i2}=-a_{i2}$
that are i.i.d. across $i=1,\dots,m$ and the joint PMF of $a_{i1},b_{i1}$
satisfies

$\displaystyle
P_{a_{i1},b_{i1}}=\begin{bmatrix}(\frac{1}{2}-\epsilon)(1-p)&(\frac{1}{2}-\epsilon)p&0\\\
2\epsilon p&0&2\epsilon(1-p)\\\
0&(\frac{1}{2}-\epsilon)(1-p)&(\frac{1}{2}-\epsilon)p\end{bmatrix}\ .$

The sum rate for distributed encoding of $({\bf A},{\bf B})$ is given as

$\displaystyle R_{\rm SW}^{\Sigma}=m(h(2\epsilon)+(1-2\epsilon)+h(p))\ .$ (9)

Exploiting Proposition 5 to compute ${\bf A}^{\intercal}{\bf B}$, we can
achieve

$\displaystyle R_{\rm KM}^{\Sigma}\leq
2mh\Big{(}2\Big{(}\frac{1}{2}-\epsilon\Big{)}(1-p)+2\epsilon(1-p),$
$\displaystyle\hskip 59.75095pt2\Big{(}\frac{1}{2}-\epsilon\Big{)}p+2\epsilon
p\Big{)}+2\log_{2}(3)\ .$ (10)

###### Proof.

We refer the reader to Appendix -J. ∎

In Figure 3-(Right), we demonstrate the sum rate performance of Proposition 5
(in log scale) versus $p$ via contrasting the sum rates $H({\bf A},{\bf B})$
and (8) for the joint PMF model in Corollary 2, where we assume that
$\epsilon=0.2$. We observe that the rate gain $\eta={R_{\rm
SW}^{\Sigma}}/{R_{\rm KM}^{\Sigma}}$ grows exponentially fast, as $p$ tends to
$\\{0,1\\}$.

Discussion. We proposed structured coding techniques for nonlinear mappings of
distributed sources in a $q$-ary prime finite field to perform inner product-
based matrix computation toward realizing distributed multiplication for
special matrix classes, e.g., symmetric and square matrices, through imposing
structural constraints on sources. Our future work includes the study of
general matrix products, providing insights into the problems of distributed
rank computation, trace computation, and low-rank matrix factorization, as
well as the derivation of tighter achievability bounds for the distributed
multiplication of general matrices and higher dimensional matrices or tensors.

## Acknowledgment

The author gratefully acknowledges the constructive discussions with Prof.
Arun Padakandla at EURECOM.

###  -D Proof of Proposition 1

The receiver aims to compute the inner product $\langle\,{\bf A},{\bf
B}\,\rangle=\sum\nolimits_{i=1}^{m}a_{i}b_{i}$, with entries from a field of
characteristic $q\geq 2$. Here, we focus on $q=2$, and the generalization to
$q>2$ is straightforward [14].

Encoding: Sources devise mappings $g_{1}:{\bf A}\to{\bf X}_{1}$ and
$g_{2}:{\bf B}\to{\bf X}_{2}$, respectively, defined below, to determine the
binary-valued column vectors

$\displaystyle{\bf X}_{1}$ $\displaystyle=g_{1}({\bf A})=\begin{bmatrix}{\bf
A}_{2}\\\ {\bf A}_{1}\\\ {\bf A}_{2}^{\intercal}{\bf
A}_{1}\end{bmatrix}\in\mathbb{F}_{2}^{(m+1)\times 1}\ ,$ (11)
$\displaystyle{\bf X}_{2}$ $\displaystyle=g_{2}({\bf B})=\begin{bmatrix}{\bf
B}_{1}\\\ {\bf B}_{2}\\\ {\bf B}_{1}^{\intercal}{\bf
B}_{2}\end{bmatrix}\in\mathbb{F}_{2}^{(m+1)\times 1}\ .$ (12)

We denote by ${\bf X}_{1}^{n},{\bf X}_{2}^{n}\in\mathbb{F}_{2}^{(m+1)\times
n}$ the length $n$ source vectors. The encoders of the sources $\\{{\bf
X}_{1i}\\}$ and $\\{{\bf X}_{2i}\\}$ are defined by functions $f_{1}:{\bf
X}_{1}^{n}\to\mathcal{R}_{f_{1}}$ and $f_{2}:{\bf
X}_{2}^{n}\to\mathcal{R}_{f_{2}}$, where $\mathcal{R}_{f_{1}}$ and
$\mathcal{R}_{f_{2}}$ denote the ranges of $f_{1}$ and $f_{2}$, respectively.
The pair of functions $(f_{1},f_{2})$ is called an $(n,\epsilon)$-coding
scheme if there exists a function
$\phi:\mathcal{R}_{f_{1}}\times\mathcal{R}_{f_{2}}\to\mathcal{Z}^{n}$ such
that by letting

$\displaystyle{\bf\hat{Z}}^{n}\triangleq\phi(f_{1}({\bf X}_{1}^{n}),f_{2}({\bf
X}_{2}^{n}))\ ,$ (13)

we have $\mathbb{P}({\bf\hat{Z}}^{n}\neq{\bf Z}^{n})<\epsilon$. Here, ${\bf
Z}$ is the modulo-two sum of ${\bf X}_{1}$ and ${\bf X}_{2}$, i.e., ${\bf
Z}={\bf X}_{1}\oplus_{2}{\bf X}_{2}\in\mathbb{F}_{2}^{(m+1)\times 1}$.

Our encoding scheme requires a well-known lemma of Elias [25], which showed
that linear codes achieve the capacity of binary symmetric channels, and its
adaption to the problem of computing the modulo-two sum of DSBSs in [11].
Using a simple generalization of this result to vector variables, for fixed
$\epsilon>0$ and for sufficiently large $n$, there exists a binary matrix
${\bf C}\in\mathbb{F}_{2}^{k\times n}$, where $f_{1}({\bf
X}_{1}^{n})\triangleq{\bf C}({\bf X}_{1}^{n})={\bf C}\cdot({\bf
X}_{1}^{n})^{\intercal}\in\mathbb{F}_{2}^{k\times(m+1)}$ and $f_{2}({\bf
X}_{2}^{n})\triangleq{\bf C}\cdot({\bf
X}_{2}^{n})^{\intercal}\in\mathbb{F}_{2}^{k\times(m+1)}$ denote the modulo-two
product of the matrix ${\bf C}$ with the transpose of the binary vector
sequences ${\bf X}_{1}^{n}$ and ${\bf X}_{2}^{n}$, respectively, and a
decoding function $\psi:\\{0,1\\}^{k\times(m+1)}\to\\{0,1\\}^{n\times(m+1)}$
that satisfy

$\displaystyle\phi(f_{1}({\bf X}_{1}^{n}),f_{2}({\bf
X}_{2}^{n}))\triangleq\psi(f_{1}({\bf X}_{1}^{n})+f_{2}({\bf X}_{2}^{n}))$

such that i) $k<n(H({\bf Z})+\epsilon)$, and ii) $\mathbb{P}(\psi({\bf C}({\bf
Z}^{n}))\neq{\bf Z}^{n})<\epsilon$. Hence, application of Elias’s lemma [25]
and [11] yields that $({\bf C},{\bf C})$ is an $(n,\epsilon)$-coding scheme.

Decoding: Exploiting the achievability result of Körner-Marton [11], the sum
rate needed for the receiver to recover the vector sequence ${\bf Z}^{n}={\bf
X}_{1}^{n}\oplus_{2}{\bf X}_{2}^{n}\in\mathbb{F}_{2}^{(m+1)\times n}$ with a
vanishing error probability can be determined as [11]:

$\displaystyle R_{\rm KM}^{\Sigma}=2H({\bf U},\,{\bf V},\,W)\ .$ (14)

Using $\psi$ given in ii), the receiver can recover ${\bf Z}^{n}$. However,
lossless decoding of ${\bf X}_{1}^{n}$ and ${\bf X}_{2}^{n}$ is not
guaranteed.

We next show that the sum rate in (14) is sufficient to recover $\langle\,{\bf
A},{\bf B}\,\rangle$. Using ${\bf\hat{Z}}^{n}$, the receiver computes

$\displaystyle{\bf U}^{\intercal}{\bf V}-W$ $\displaystyle=({\bf
A}_{2}^{\intercal}\oplus_{q}{\bf B}_{1}^{\intercal})({\bf A}_{1}\oplus_{q}{\bf
B}_{2})$ $\displaystyle-({\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{q}{\bf
B}_{1}^{\intercal}{\bf B}_{2})$ $\displaystyle={\bf B}_{1}^{\intercal}{\bf
A}_{1}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf B}_{2}$
$\displaystyle\overset{(a)}{=}{\bf A}_{1}^{\intercal}{\bf B}_{1}\oplus_{q}{\bf
A}_{2}^{\intercal}{\bf B}_{2}=\langle\,{\bf A},{\bf B}\,\rangle$
$\displaystyle\ ,$ (15)

where $(a)$ follows from ${\bf B}_{1}^{\intercal}{\bf A}_{1}={\bf
A}_{1}^{\intercal}{\bf B}_{1}\in\mathbb{F}_{q}$.

Hence, (14) with (-D) gives the achievability result we seek.

###  -E Proof of Corollary 1

Employing the definitions of ${\bf X}_{1}$ and ${\bf X}_{2}$ in (11), (1) and
[11], we can determine the sum rate needed for the receiver to recover $({\bf
U},\,{\bf V},\,W)$ in an asymptotic manner:

$\displaystyle R_{\rm KM}^{\Sigma}$ $\displaystyle=2H({\bf U},\,{\bf V},\,W)$
$\displaystyle=2H({\bf U})+2H({\bf V})+2H({\bf A}_{2}^{\intercal}{\bf
A}_{1}\oplus_{2}{\bf B}_{1}^{\intercal}{\bf B}_{2}\,|\,{\bf U},{\bf V})$
$\displaystyle=2\frac{m}{2}h(p)+2\frac{m}{2}h(p)+2H({\bf
A}_{2}^{\intercal}{\bf A}_{1}\oplus_{2}{\bf B}_{1}^{\intercal}{\bf
B}_{2}\,|\,{\bf U},{\bf V})$ $\displaystyle\overset{(a)}{=}2mh(p)+2H({\bf
U}^{\intercal}{\bf A}_{1}\oplus_{2}{\bf A}_{2}^{\intercal}{\bf V}\,|\,{\bf
U},{\bf V})$ $\displaystyle\overset{(b)}{=}2mh(p)+2H({\bf Q}^{\intercal}{\bf
A}\,|\,{\bf U},{\bf V})$
$\displaystyle\overset{(c)}{=}2mh(p)+2H\Big{(}\sum\limits_{i=1}^{m}a_{i}z_{i}\,|\,{\bf
Q}\Big{)}$
$\displaystyle\overset{(d)}{=}2mh(p)+2\sum\limits_{j=1}^{m}\binom{m}{j}p^{j}(1-p)^{m-j}H\Big{(}\sum\limits_{i=1}^{j}a_{i}\Big{)}$
$\displaystyle\overset{(e)}{=}2mh(p)+2(1-(1-p)^{m})\ ,$ (16)

where $(a)$ follows from employing the relations ${\bf U}={\bf
A}_{2}\oplus_{2}{\bf
B}_{1}=\begin{bmatrix}u_{1}&u_{2}&\ldots&u_{m/2}\end{bmatrix}^{\intercal}$ and
${\bf V}={\bf A}_{1}\oplus_{2}{\bf
B}_{2}=\begin{bmatrix}v_{1}&v_{2}&\ldots&v_{m/2}\end{bmatrix}^{\intercal}$,
and simplification using conditioning, $(b)$ follows from employing ${\bf
Q}=\begin{bmatrix}{\bf U}\\\ {\bf V}\end{bmatrix}\in\mathbb{F}_{2}^{m\times
1}$ and ${\bf A}_{2}^{\intercal}{\bf V}={\bf V}^{\intercal}{\bf A}_{2}$, $(c)$
from the definition of ${\bf Q}$, and $(d)$ from $H({\bf Q}^{\intercal}{\bf
A}\,|\,{\bf Q})\leq H({\bf Q}^{\intercal}{\bf A}\,|\,{\bf
Q}^{\intercal}{\bm{1}}_{m})$, where ${\bf
Q}^{\intercal}{\bm{1}}_{m}=\sum\limits_{i=1}^{m}y_{i}\sim B(m,p)$, such that
$H({\bf Q}^{\intercal}{\bf A}\,|\,{\bf
Q}^{\intercal}{\bm{1}}_{m}=j)=H\Big{(}\sum\limits_{i=1}^{j}a_{i}\Big{)}$,
$j\geq 1$ exploiting that $a_{i}\overset{i.i.d.}{\sim}{\rm Bern}(1/2)$, and
$H({\bf Q}^{\intercal}{\bf A}\,|\,{\bf Q}^{\intercal}{\bm{1}}_{m}=0)=0$.
Finally, incorporating $H\Big{(}\sum\limits_{i=1}^{j}a_{i}\Big{)}=1$, $j\geq
1$ we obtain $(e)$.

The encoding rate for asymptotic lossless compression of ${\bf A}$ and ${\bf
B}$ is given by the Slepian-Wolf theorem [10]:

$\displaystyle R_{\rm SW}^{\Sigma}$ $\displaystyle=H({\bf A},{\bf B})$
$\displaystyle\overset{(a)}{=}2H({\bf A}_{1},{\bf B}_{2})$
$\displaystyle=2\frac{m}{2}(1+h(p))$ $\displaystyle=m(1+h(p))\ ,$ (17)

where $(a)$ follows from using ${\bf A}=\begin{bmatrix}{\bf A}_{1}\\\ {\bf
A}_{2}\end{bmatrix}$ and ${\bf B}=\begin{bmatrix}{\bf B}_{1}\\\ {\bf
B}_{2}\end{bmatrix}$, with ${\bf A}_{1}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}{\bf A}_{2}$, ${\bf
B}_{1}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}{\bf B}_{2}$, and ${\bf A}$ and ${\bf B}$
having i.i.d. entries.

From (-E) and (-E), $\eta={R_{\rm SW}^{\Sigma}}/{R_{\rm KM}^{\Sigma}}$ is
given by (3).

###  -F Proof of Proposition 2

Here, we provide a sketch of the proof. For binary sources, i.e., $q=2$, it is
easy to verify that

$\displaystyle\langle\,{\bf A},{\bf
B}\,\rangle=\Big{\lfloor}\frac{1}{2}\Big{(}\sum\limits_{i=1}^{m}a_{i}\oplus_{r}b_{i}-\sum\limits_{i=1}^{m}(a_{i}\oplus_{2}b_{i})\Big{)}\Big{\rfloor}\mod
2\ ,$

where $r=2m$ for $m$ even, and $r=2m+1$ for $m$ odd, respectively. Hence, the
following sum rate is achievable:

$\displaystyle R_{\rm
S}^{\Sigma}=2H\Big{(}\sum\limits_{i=1}^{m}a_{i}\oplus_{r}b_{i},\,\\{a_{i}\oplus_{2}b_{i}\\}_{i=1}^{m}\Big{)}\
.$ (18)

When the data is generated by two correlated memoryless $q$-ary sources for
$q\geq 2$, it is possible to achieve a sum rate

$\displaystyle R_{\rm
S}^{\Sigma}=2H\Big{(}\Big{\\{}a_{i}\oplus_{r}b_{i}\Big{\\}}_{i=1}^{m},\,{\bigoplus_{i=1}^{m}}{{}_{{}_{{\operatorname*{{\small
q}}}}}\,a_{i}^{2}\oplus_{q}b_{i}^{2}}\Big{)}\ ,$

where $r=2(q-1)m$ and $r=2(q-1)m+1$ for even and odd $m$, respectively. Upon
receiving $\big{\\{}a_{i}\oplus_{r}b_{i}\big{\\}}_{i=1}^{m}$ and
${\bigoplus\limits_{i=1}^{m}}{{}_{{}_{{\operatorname*{{\small
q}}}}}\,a_{i}^{2}\oplus_{q}b_{i}^{2}}$, the receiver can reconstruct

$\displaystyle
2c=qk+\Big{(}\sum\limits_{i=1}^{m}(a_{i}\oplus_{r}b_{i})^{2}-\,{\bigoplus\limits_{i=1}^{m}}{{}_{{}_{{\operatorname*{{\small
q}}}}}\,(a_{i}^{2}\oplus_{q}b_{i}^{2})\Big{)}}\mod q\ ,$

where there is a unique $k\in\mathbb{F}_{q}$ for which $c=\langle\,{\bf
A},{\bf B}\,\rangle\in\mathbb{F}_{q}$.

###  -G Proof of Proposition 3

Given two sequences of random matrices ${\bf A},\,{\bf
B}\in\mathbb{F}_{q}^{m\times l}$, the receiver aims to compute ${\bf
A}^{\intercal}{\bf B}={\bf B}^{\intercal}{\bf A}\in\mathbb{F}_{q}^{l\times
l}$.

Encoding: Each source uses mappings $g_{1}:{\bf A}\to{\bf X}_{1}$ and
$g_{2}:{\bf B}\to{\bf X}_{2}$, to determine the respective matrices:

$\displaystyle{\bf X}_{1}$ $\displaystyle=g_{1}({\bf A})=\begin{bmatrix}{\bf
A}_{2}\\\ {\bf A}_{1}\\\ {\bf A}_{2}^{\intercal}{\bf
A}_{1}\end{bmatrix}\in\mathbb{F}_{q}^{(m+l)\times l}\ ,$ (19)
$\displaystyle{\bf X}_{2}$ $\displaystyle=g_{2}({\bf B})=\begin{bmatrix}{\bf
B}_{1}\\\ {\bf B}_{2}\\\ {\bf B}_{1}^{\intercal}{\bf
B}_{2}\end{bmatrix}\in\mathbb{F}_{q}^{(m+l)\times l}\ .$ (20)

Following the steps of Appendix -D, there exists an $(n,\epsilon)$-coding
scheme $({\bf C},{\bf C})$ for a matrix ${\bf C}\in\mathbb{F}_{q}^{k\times
n}$, where $k=(m+l)\times l$, for decoding ${\bf Z}={\bf X}_{1}\oplus_{q}{\bf
X}_{2}\in\mathbb{F}_{q}^{(m+l)\times l}$ with a small probability of error.

Decoding: Exploiting the achievability result of Körner-Marton [11], the sum
rate needed for the receiver to recover the matrix sequence ${\bf Z}^{n}={\bf
X}_{1}^{n}\oplus_{q}{\bf X}_{2}^{n}$ with a vanishing error probability can be
determined as [11]:

$\displaystyle R_{\rm KM}^{\Sigma}=2H({\bf U},\,{\bf V},\,{\bf W})\ .$ (21)

For $q>2$, given a symmetric matrix ${\bf D}={\bf A}^{\intercal}{\bf
B}\in\mathbb{F}_{q}^{l\times l}$, the following relation holds:

$\displaystyle{\bf D}=\frac{1}{2}({\bf D}\oplus_{q}{\bf D}^{\intercal})={\bf
A}^{\intercal}{\bf B}={\bf B}^{\intercal}{\bf A}\ .$ (22)

Using ${\bf\hat{Z}}^{n}$, the receiver computes

$\displaystyle\frac{1}{2}(({\bf U}^{\intercal}{\bf V}-{\bf W})\oplus_{q}({\bf
U}^{\intercal}{\bf V}-{\bf W})^{\intercal})$ $\displaystyle=\frac{1}{2}(({\bf
A}_{2}^{\intercal}\oplus_{q}{\bf B}_{1}^{\intercal})({\bf A}_{1}\oplus_{q}{\bf
B}_{2})-({\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{q}{\bf
B}_{1}^{\intercal}{\bf B}_{2})$ $\displaystyle\hskip 12.80365pt\oplus_{q}({\bf
A}_{1}^{\intercal}\oplus_{q}{\bf B}_{2}^{\intercal})({\bf A}_{2}\oplus_{q}{\bf
B}_{1})-({\bf A}_{1}^{\intercal}{\bf A}_{2}\oplus_{q}{\bf
B}_{2}^{\intercal}{\bf B}_{1}))$ $\displaystyle=\frac{1}{2}(({\bf
B}_{1}^{\intercal}{\bf A}_{1}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf
B}_{2})\oplus_{q}({\bf A}_{1}^{\intercal}{\bf B}_{1}\oplus_{q}{\bf
B}_{2}^{\intercal}{\bf A}_{2}))$ $\displaystyle=\frac{1}{2}(({\bf
A}_{1}^{\intercal}{\bf B}_{1}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf
B}_{2})\oplus_{q}({\bf A}_{1}^{\intercal}{\bf B}_{1}\oplus_{q}{\bf
A}_{2}^{\intercal}{\bf B}_{2})^{\intercal})$ $\displaystyle=\frac{1}{2}({\bf
A}^{\intercal}{\bf B}\oplus_{q}({\bf A}^{\intercal}{\bf B})^{\intercal})={\bf
A}^{\intercal}{\bf B}\ ,$ (23)

where the last equality follows from that ${\bf D}={\bf A}^{\intercal}{\bf B}$
is a symmetric matrix and it satisfies (22).

Combining (21) with (-G) gives the achievability result.

###  -H Proof of Proposition 4

We first show that the encoding scheme of Proposition 1 does not allow the
recovery of ${\bf A},{\bf B}\in\mathbb{F}_{q}^{m\times 1}$ by the receiver,
i.e., $H({\bf A},{\bf B}\,|\,{\bf A}^{\intercal}{\bf B},{\bf Q})>0$, for
$m>1$.

The receiver can recover ${\bf Q}=\begin{bmatrix}{\bf U}\\\ {\bf
V}\end{bmatrix}\in\mathbb{F}_{q}^{m\times 1}$ and $W\in\mathbb{F}_{q}$ with a
small probability of error. The extra rate needed from the encoders for the
receiver to determine ${\bf A},\,{\bf B}\in\mathbb{F}_{q}^{m\times 1}$ is

$\displaystyle H({\bf A},{\bf B}\,|\,{\bf A}^{\intercal}{\bf B},{\bf Q})$
$\displaystyle\overset{(a)}{=}H({\bf A},{\bf B},{\bf A}^{\intercal}{\bf
B},{\bf Q})-H({\bf A}^{\intercal}{\bf B},{\bf Q})$
$\displaystyle\overset{(b)}{=}H({\bf A},{\bf B},{\bf Q})-H({\bf
A}^{\intercal}{\bf B},{\bf Q})$ $\displaystyle\overset{(c)}{=}H({\bf A},{\bf
Q})-H({\bf A}^{\intercal}{\bf Q},{\bf Q})$
$\displaystyle\overset{(d)}{=}H({\bf A},{\bf A}^{\intercal}{\bf Q},{\bf
Q})-H({\bf A}^{\intercal}{\bf Q},{\bf Q})$
$\displaystyle\overset{(e)}{=}H({\bf A}\,|\,{\bf A}^{\intercal}{\bf Q},{\bf
Q})\overset{(f)}{\geq}0\ ,$

where $(a)$ follows from using the definition of conditional entropy and that
the receiver can compute ${\bf A}^{\intercal}{\bf B}={\bf U}^{\intercal}{\bf
V}-W$ from ${\bf Q}$ and $W$, $(b)$ from $H({\bf A}^{\intercal}{\bf
B}\,|\,{\bf Q},{\bf A},{\bf B})=0$, $(c)$ from $H({\bf B}\,|\,{\bf A},{\bf
Q})=0$ and ${\bf A}^{\intercal}{\bf B}={\bf A}^{\intercal}{\bf Q}$ given ${\bf
Q}$, $(d)$ from $H({\bf A}^{\intercal}{\bf Q}\,|\,{\bf A},{\bf Q})=0$, $(e)$
from employing the definition of conditional entropy, and $(f)$ holds with
equality if the function $g({\bf A},{\bf Q})={\bf A}^{\intercal}{\bf Q}$ is
partially invertible, i.e., $H({\bf A}\,|\,g({\bf A},{\bf Q}),{\bf Q})=0$,
e.g., when $g$ is the arithmetic sum or the modulo sum of the two vectors.
Hence, the inequality in $(f)$ is strict for inferring ${\bf A},\,{\bf
B}\in\mathbb{F}_{q}^{m\times 1}$ from ${\bf Q}$ and $W$.

We next prove the main result of the proposition. Given matrix variables ${\bf
A},{\bf B}\in\mathbb{F}_{q}^{m\times l}$ such that ${\bf D}={\bf
A}^{\intercal}{\bf B}$ is symmetric, letting ${\bf Q}=\begin{bmatrix}{\bf
U}\\\ {\bf V}\end{bmatrix}$, we first expand $R_{\rm KM}^{\Sigma}$ as

$\displaystyle R_{\rm KM}^{\Sigma}$ $\displaystyle=2H({\bf U},\,{\bf V},\,W)$
$\displaystyle=2H({\bf U},\,{\bf V},\,{\bf A}^{\intercal}{\bf B})$
$\displaystyle=2H({\bf A}^{\intercal}{\bf B})+2H({\bf Q}\,|\,{\bf
A}^{\intercal}{\bf B})\ .$ (24)

We next expand $R_{\rm SW}^{\Sigma}$ as

$\displaystyle R_{\rm SW}^{\Sigma}=H({\bf A},{\bf B})=H({\bf A},{\bf B},{\bf
U},\,{\bf V},\,W,{\bf A}^{\intercal}{\bf B})$ $\displaystyle=H({\bf U},\,{\bf
V},\,{\bf A}^{\intercal}{\bf B})+H({\bf A}_{1},{\bf B}_{1},{\bf A}_{2},{\bf
B}_{2}\,|\,{\bf A}_{2}\oplus_{q}{\bf B}_{1},\,$ $\displaystyle\hskip
108.12054pt{\bf A}_{1}\oplus_{q}{\bf B}_{2},{\bf A}_{1}^{\intercal}{\bf
B}_{1}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf B}_{2})$ $\displaystyle=H({\bf
U},\,{\bf V},\,{\bf A}^{\intercal}{\bf B})+H({\bf A}_{1},{\bf A}_{2}\,|\,{\bf
U},\,{\bf V},{\bf A}_{1}^{\intercal}({\bf U}\oplus_{q}{\bf A}_{2})$
$\displaystyle\hskip 142.26378pt\oplus_{q}{\bf A}_{2}^{\intercal}({\bf
V}\oplus_{q}{\bf A}_{1}))$ $\displaystyle=H({\bf Q},\,{\bf A}^{\intercal}{\bf
B})+H({\bf A}\,|\,{\bf Q},{\bf A}^{\intercal}{\bf Q}\oplus_{q}{\bf
A}_{1}^{\intercal}{\bf A}_{2}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf A}_{1})$
$\displaystyle=H({\bf A}^{\intercal}{\bf B})+H({\bf Q}\,|\,{\bf
A}^{\intercal}{\bf B})$ $\displaystyle+H({\bf A}\,|\,{\bf Q},{\bf
A}^{\intercal}{\bf Q}\oplus_{q}{\bf A}_{1}^{\intercal}{\bf
A}_{2}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf A}_{1})\ ,$ (25)

where it is easy to observe that ${\bf A}^{\intercal}{\bf B}={\bf
A}_{1}^{\intercal}({\bf U}+{\bf A}_{2})+{\bf A}_{2}^{\intercal}({\bf V}+{\bf
A}_{1})={\bf A}^{\intercal}{\bf Q}\oplus_{q}{\bf A}_{1}^{\intercal}{\bf
A}_{2}\oplus_{q}{\bf A}_{2}^{\intercal}{\bf A}_{1}$.

Similarly, via exploiting $\tilde{\bf Q}=\begin{bmatrix}{\bf V}\\\ {\bf
U}\end{bmatrix}$, we can show that

$\displaystyle R_{\rm SW}^{\Sigma}$ $\displaystyle=H({\bf A}^{\intercal}{\bf
B})+H({\bf Q}\,|\,{\bf A}^{\intercal}{\bf B})$ $\displaystyle+H({\bf
B}\,|\,\tilde{\bf Q},\tilde{\bf Q}^{\intercal}{\bf B}\oplus_{q}{\bf
B}_{1}^{\intercal}{\bf B}_{2}+{\bf B}_{2}^{\intercal}{\bf B}_{1})\ ,$ (26)

where ${\bf A}^{\intercal}{\bf B}=\tilde{\bf Q}^{\intercal}{\bf
B}\oplus_{q}{\bf B}_{1}^{\intercal}{\bf B}_{2}\oplus_{q}{\bf
B}_{2}^{\intercal}{\bf B}_{1}$.

From (-H) and (-H), we note that $H({\bf A}\,|\,{\bf A}^{\intercal}{\bf
Q})=H({\bf B}\,|\,\tilde{\bf Q}^{\intercal}{\bf B})$. Contrasting (-H) with
(-H), the condition

$\displaystyle H({\bf A}^{\intercal}{\bf B})+H({\bf Q}\,|\,{\bf
A}^{\intercal}{\bf B})$ $\displaystyle<H({\bf A}\,|\,{\bf Q},{\bf
A}^{\intercal}{\bf B})$ $\displaystyle=H({\bf B}\,|\,\tilde{\bf Q},{\bf
A}^{\intercal}{\bf B})$

ensures that $R_{\rm KM}^{\Sigma}<R_{\rm SW}^{\Sigma}$.

When $l=1$, we have ${\bf A}_{1}^{\intercal}{\bf A}_{2}={\bf
A}_{2}^{\intercal}{\bf A}_{1}$ and ${\bf B}_{1}^{\intercal}{\bf B}_{2}={\bf
B}_{2}^{\intercal}{\bf B}_{1}$, hence, the above condition is equivalent to

$\displaystyle H({\bf A}^{\intercal}{\bf B})+H({\bf Q}\,|\,{\bf
A}^{\intercal}{\bf B})$ $\displaystyle<H({\bf A}\,|\,{\bf Q},{\bf
A}^{\intercal}{\bf Q})$ $\displaystyle=H({\bf B}\,|\,\tilde{\bf Q},\tilde{\bf
Q}^{\intercal}{\bf B})\ .$

###  -I Proof of Proposition 5

Note that ${\bf D}_{j}={\bf A}^{\intercal}{\bf B}_{j}$ for
$j\in\\{1,\dots,l\\}$, where ${\bf
D}_{j}=\begin{bmatrix}d_{1j}&d_{2j}&\ldots&d_{lj}\end{bmatrix}^{\intercal}\in\mathbb{F}_{q}^{l\times
1}$. Following the steps in Appendix -D, the receiver can recover $\\{{\bf
A}\oplus_{q}{\bf\tilde{B}}_{j}\\}_{j=1}^{l}\,,\,\\{{\bf A}^{\intercal}{\bf
A}\oplus_{q}{\bf\tilde{B}}_{j}^{\intercal}{\bf\tilde{B}}_{j}\\}_{j=1}^{l}$,
and then compute the following $l\times l$ matrix:

$\displaystyle({\bf A}\oplus_{q}{\bf\tilde{B}}_{j})^{\intercal}({\bf
A}\oplus_{q}{\bf\tilde{B}}_{j})-({\bf A}^{\intercal}{\bf
A}\oplus_{q}{\bf\tilde{B}}_{j}^{\intercal}{\bf\tilde{B}}_{j})$
$\displaystyle={\bf
A}^{\intercal}{\bf\tilde{B}}_{j}\oplus_{q}{\bf\tilde{B}}_{j}^{\intercal}{\bf
A}$
$\displaystyle=\begin{bmatrix}d_{1j}\oplus_{q}d_{1j}&d_{1j}\oplus_{q}d_{2j}&\ldots&d_{1j}\oplus_{q}d_{lj}\\\
d_{2j}\oplus_{q}d_{1j}&d_{2j}\oplus_{q}d_{2j}&\ldots&d_{2j}\oplus_{q}d_{lj}\\\
\vdots&\vdots&\ddots&\vdots\\\
d_{lj}\oplus_{q}d_{1j}&d_{lj}\oplus_{q}d_{2j}&\ldots&d_{lj}\oplus_{q}d_{lj}\end{bmatrix}\
,$

which is a symmetric matrix with $l$ unknowns and $\frac{l(l-1)}{2}\geq l$
linearly independent equations for $l\geq 2$ and $q>2$. Hence, ${\bf D}_{j}$,
for each $j=1,\dots,l$, as well as ${\bf D}$ can be recovered.

###  -J Proof of Corollary 2

The sum rate for distributed encoding of $({\bf A},{\bf B})$ is given as

$\displaystyle R_{\rm SW}^{\Sigma}=H({\bf A},{\bf B})$ $\displaystyle=H({\bf
A}_{1},{\bf B}_{1},{\bf A}_{2},{\bf B}_{2})$ $\displaystyle=H({\bf A}_{1},{\bf
B}_{1})$ $\displaystyle=m(h(2\epsilon)+(1-2\epsilon)+h(p))\ .$

Exploiting Proposition 5 to compute ${\bf A}^{\intercal}{\bf B}$, we can
achieve

$\displaystyle R_{\rm KM}^{\Sigma}$ $\displaystyle=2H({\bf
A}\oplus_{3}{\bf\tilde{B}}_{1}\,,\,{\bf A}^{\intercal}{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1}^{\intercal}{\bf\tilde{B}}_{1})$
$\displaystyle=2mh\Big{(}2\Big{(}\frac{1}{2}-\epsilon\Big{)}(1-p)+2\epsilon(1-p),$
$\displaystyle\hskip 85.35826pt2\Big{(}\frac{1}{2}-\epsilon\Big{)}p+2\epsilon
p\Big{)}$ $\displaystyle+2h({\bf A}^{\intercal}{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1}^{\intercal}{\bf\tilde{B}}_{1}\,|\,{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1})$ $\displaystyle\leq
2mh\Big{(}2\Big{(}\frac{1}{2}-\epsilon\Big{)}(1-p)+2\epsilon(1-p),$
$\displaystyle\hskip 39.83368pt2\Big{(}\frac{1}{2}-\epsilon\Big{)}p+2\epsilon
p\Big{)}+2\log_{2}(3)\ ,$

where the last step follows from using that

$\displaystyle{\bf A}^{\intercal}{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1}^{\intercal}{\bf\tilde{B}}_{1}=\begin{bmatrix}{\bf
A}_{1}^{\intercal}{\bf A}_{1}\oplus_{3}{\bf B}_{1}^{\intercal}{\bf B}_{1}&{\bf
A}_{1}^{\intercal}{\bf A}_{2}\oplus_{3}{\bf B}_{1}^{\intercal}{\bf B}_{1}\\\
{\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{3}{\bf B}_{1}^{\intercal}{\bf
B}_{1}&{\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{3}{\bf B}_{1}^{\intercal}{\bf
B}_{1}\end{bmatrix}\ ,$

and evaluating the conditional entropy as

$\displaystyle h({\bf A}^{\intercal}$ $\displaystyle{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1}^{\intercal}{\bf\tilde{B}}_{1}\,|\,{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1})$ $\displaystyle=h({\bf A}_{1}^{\intercal}{\bf
A}_{1}\oplus_{3}{\bf B}_{1}^{\intercal}{\bf B}_{1},{\bf A}_{1}^{\intercal}{\bf
A}_{2}\oplus_{3}{\bf B}_{1}^{\intercal}{\bf B}_{1},$ $\displaystyle\hskip
112.38829pt{\bf A}_{2}^{\intercal}{\bf A}_{1}\oplus_{3}{\bf
B}_{1}^{\intercal}{\bf B}_{1}\,|\,{\bf A}\oplus_{3}{\bf\tilde{B}}_{1})$
$\displaystyle=h(\big{\\{}\sum\limits_{i=1}^{m}a_{ij}^{2}\oplus_{3}b_{i1}^{2}\big{\\}}_{j=1}^{2},\sum\limits_{i=1}^{m}a_{i1}a_{i2}\oplus_{3}b_{i1}^{2}\,|\,{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1})$
$\displaystyle\overset{(a)}{=}h(\sum\limits_{i=1}^{m}a_{i1}^{2}\oplus_{3}b_{i1}^{2},\sum\limits_{i=1}^{m}a_{i2}^{2}\oplus_{3}b_{i1}^{2}\,|\,{\bf
A}\oplus_{3}{\bf\tilde{B}}_{1})$
$\displaystyle\overset{(b)}{\leq}2\log_{2}(3)\ ,$

where $(a)$ follows from that $a_{i1}a_{i2}\oplus_{3}b_{i1}^{2}$ can be
recovered given the side information ${\bf
A}\oplus_{3}{\bf\tilde{B}}_{1}=\big{\\{}a_{ij}\oplus_{3}b_{i1}\big{\\}}_{i,j}$
for $i=1,2,\dots,m$ and $j=1,2$, and given
$\big{\\{}a_{ij}^{2}\oplus_{3}b_{i1}^{2}\big{\\}}_{j}$, for $j=1,2$. More
specifically, the receiver can recover

$\displaystyle
2a_{ij}b_{i1}=(a_{ij}\oplus_{3}b_{i1})^{2}-(a_{ij}^{2}\oplus_{3}b_{i1}^{2})\
,\quad j=1,2\ .$ (27)

Hence, using the side information ${\bf A}\oplus_{3}{\bf\tilde{B}}_{1}$ and
(27), the receiver can recover $a_{i1}a_{i2}\oplus_{3}b_{i1}^{2}$ as follows:

$\displaystyle a_{i1}a_{i2}\oplus_{3}b_{i1}^{2}$
$\displaystyle=(a_{i1}\oplus_{3}b_{i1})(a_{i2}\oplus_{3}b_{i1})$
$\displaystyle-(a_{i1}b_{i1}\oplus_{3}a_{i2}b_{i1})\ .$

Finally, step $(b)$ follows from exploiting that both variable
$\sum\limits_{i=1}^{m}a_{i1}^{2}\oplus_{3}b_{i1}^{2}$ and variable
$\sum\limits_{i=1}^{m}a_{i2}^{2}\oplus_{3}b_{i1}^{2}$ reside in
$\mathbb{F}_{3}$.

## References

  * [1] G. Strang, _Introduction to linear algebra_. SIAM, 2022.
  * [2] M. H. Mousavi, M. A. Maddah-Ali, and M. Mirmohseni, “Private inner product retrieval for distributed machine learning,” in _Proc., IEEE Int. Symp. Inf. Theory (ISIT)_ , Paris, France, Jul. 2019, pp. 355–359.
  * [3] É. Bouscatié, G. Castagnos, and O. Sanders, “Pattern matching in encrypted stream from inner product encryption,” in _Proc., IACR Int. Conf. Public-Key Cryptogr._ Springer, May 2023, pp. 774–801.
  * [4] H. Buhrman, R. Cleve, J. Watrous, and R. De Wolf, “Quantum fingerprinting,” _Physical Review Letters_ , vol. 87, no. 16, p. 167902, Sep. 2001.
  * [5] R. Tandon, Q. Lei, A. G. Dimakis, and N. Karampatziakis, “Gradient coding: Avoiding stragglers in distributed learning,” in _International Conference on Machine Learning_. PMLR, Jul. 2017, pp. 3368–3376.
  * [6] Q. Yu, M. Maddah-Ali, and S. Avestimehr, “Polynomial codes: an optimal design for high-dimensional coded matrix multiplication,” in _Proc., Advances in Neural Information Processing Systems_ , vol. 30, Long Beach, CA, Dec. 2017, pp. 4403––4413.
  * [7] S. Dutta, M. Fahim, F. Haddadpour, H. Jeong, V. Cadambe, and P. Grover, “On the optimal recovery threshold of coded matrix multiplication,” _IEEE Trans. Inf. Theory_ , vol. 66, no. 1, pp. 278–301, Jul. 2019.
  * [8] M. Soleymani, H. Mahdavifar, and A. S. Avestimehr, “Analog lagrange coded computing,” _IEEE J. Sel. Areas Inf. Theory_ , vol. 2, no. 1, pp. 283–295, Feb. 2021.
  * [9] Q. Yu, S. Li, N. Raviv, S. M. M. Kalan, M. Soltanolkotabi, and S. A. Avestimehr, “Lagrange coded computing: Optimal design for resiliency, security, and privacy,” in _Proc., Int. Conf. Artif. Intell. Stat._ Naha, Okinawa, Japan: PMLR, Apr. 2019, pp. 1215–1225.
  * [10] D. Slepian and J. K. Wolf, “Noiseless coding of correlated information sources,” _IEEE Trans. Inf. Theory_ , vol. 19, no. 4, pp. 471–480, Jul. 1973\.
  * [11] J. Körner and K. Marton, “How to encode the modulo-two sum of binary sources (corresp.),” _IEEE Trans. Inf. Theory_ , vol. 25, no. 2, pp. 219–221, Mar. 1979.
  * [12] D. Krithivasan and S. S. Pradhan, “Distributed source coding using abelian group codes: A new achievable rate-distortion region,” _IEEE Trans. Inf. Theory_ , vol. 57, no. 3, pp. 1495–1519, Feb. 2011.
  * [13] S. S. Pradhan, A. Padakandla, F. Shirani _et al._ , “An algebraic and probabilistic framework for network information theory,” _Found. Trends Commun. Inf. Theory_ , vol. 18, no. 2, pp. 173–379, Dec. 2020.
  * [14] T. Han and K. Kobayashi, “A dichotomy of functions F(X, Y) of correlated sources (X, Y),” _IEEE Trans. Inf. Theory_ , vol. 33, no. 1, pp. 69–76, Jan. 1987.
  * [15] A. Orlitsky and J. R. Roche, “Coding for computing,” _IEEE Trans. Inf. Theory_ , vol. 47, no. 3, p. 903–917, Mar. 2001.
  * [16] J. Körner, “Coding of an information source having ambiguous alphabet and the entropy of graphs,” in _Proc., 6th Prague Conf. Inf. Theory_ , Prague, Czech Republic, Sep. 1973, pp. 411–425.
  * [17] S. Feizi and M. Médard, “On network functional compression,” _IEEE Trans. Inf. Theory_ , vol. 60, no. 9, pp. 5387–5401, Sep. 2014.
  * [18] D. Malak, “Fractional graph coloring for functional compression with side information,” in _Proc., IEEE Inf. Theory Wksh._ , Mumbai, India, Nov. 2022, pp. 750–755.
  * [19] M. R. Deylam Salehi and D. Malak, “An achievable low complexity encoding scheme for coloring cyclic graphs,” in _Proc., Annu. Allerton Conf. Commun. Control Comput._ , Monticello, IL, USA, Sep. 2023, pp. 1–8.
  * [20] R. Ahlswede and T. Han, “On source coding with side information via a multiple-access channel and related problems in multi-user information theory,” _IEEE Trans. Inf. Theory_ , vol. 29, no. 3, pp. 396–412, May 1983\.
  * [21] C. Nair and Y. N. Wang, “On optimal weighted-sum rates for the modulo sum problem,” in _Proc., IEEE Int. Symp. Inf. Theory_ , Jun. 2020, pp. 2416–2420.
  * [22] W.-T. Chang and R. Tandon, “On the capacity of secure distributed matrix multiplication,” in _Proc., IEEE Global Commun. Conf._ , Abu Dhabi, UAE, Dec. 2018, pp. 1–6.
  * [23] Z. Jia and S. A. Jafar, “On the capacity of secure distributed batch matrix multiplication,” _IEEE Trans. Inf. Theory_ , vol. 67, no. 11, pp. 7420–7437, Sep. 2021.
  * [24] Y. Zhao and H. Sun, “Expand-and-randomize: An algebraic approach to secure computation,” _Entropy_ , vol. 23, no. 11, p. 1461, Nov. 2021.
  * [25] R. G. Gallager, _Information theory and reliable communication_. Wiley: New York, Jan. 1968, vol. 588.


# Automatic Detection of Myocontrol Failures
Based upon Situational Context Information

Karoline Heiwolt1, Claudio Zito2, Markus Nowak3, Claudio Castellini3 and
Rustam Stolkin2 *This work was supported by UK Engineering and Physical
Sciences Research Council (EPSRC No. EP/R02572X/1).1Karoline Heiwolt is with
the Computational Neuroscience and Cognitive Robotics (CNCR), University of
Birmingham, UK<EMAIL_ADDRESS>Zito and Rustam Stolking are
with the Extreme Robotics Lab (ERL), University of Birmingham, UK
<EMAIL_ADDRESS>Nowak and Claudio Castellini are with the
Institute of Robotics and Mechatronics, DLR — German Aerospace Center, Germany
<EMAIL_ADDRESS>
###### Abstract

Myoelectric control systems for assistive devices are still unreliable. The
user’s input signals can become unstable over time due to e.g. fatigue,
electrode displacement, or sweat. Hence, such controllers need to be
constantly updated and heavily rely on user feedback. In this paper, we
present an automatic failure detection method which learns when plausible
predictions become unreliable and model updates are necessary. Our key insight
is to enhance the control system with a set of generative models that learn
sensible behaviour for a desired task from human demonstration. We illustrate
our approach on a grasping scenario in Virtual Reality, in which the user is
asked to grasp a bottle on a table. From demonstration our model learns the
reach-to-grasp motion from a resting position to two grasps (power grasp and
tridigital grasp) and how to predict the most adequate grasp from local
context, e.g. tridigital grasp on the bottle cap or around the bottleneck. By
measuring the error between new grasp attempts and the model prediction, the
system can effectively detect which input commands do not reflect the user’s
intention. We evaluated our model in two cases: i) with both position and
rotation information of the wrist pose, and ii) with only rotational
information. Our results show that our approach detects statistically highly
significant differences in error distributions with $\mathbf{p<0.001}$ between
successful and failed grasp attempts in both cases.

##  I INTRODUCTION

When designing a prosthesis, we aim to restore the full functionality of the
missing body part, but as the functionality increases, it becomes more
challenging to reliably control the device [13, 10, 5, 8, 4]. The demand for
advanced and dexterous control systems among patients is high. And yet, the
available myocontrol options are not very popular, with the main issue being
reliability [13, 8, 7].

To overcome this obstacle, Gijsberts et al. [6] proposed a supervised
incremental learning method, called iRR-RFF, that allows continuous adjustment
of the control model with minimal training effort. iRR-RFF is effective to
retain accurate control across signal shifts [15, 1], and the computational
cost of each update is not dependent on the amount of previously trained data,
so that model updates remain consistently quick. This is an important
efficiency benefit, that makes the method suitable for daily use. Although
this mechanism stabilises reliability over time, it still relies heavily on
human-driven feedback. The user will need to permanently monitor performance
and intervene for updates. In order to facilitate day-to-day use of an
electric prosthesis, it would be desirable to develop a complete system, that
assesses its state automatically and triggers updates when needed. In an
attempt to automate the incremental model updates, the authors in [12] trained
a standard linear classifier to detect failures in the sEMG input signal. The
model classified examples of myocontrol use in different tasks as good versus
poor control performance, and was able to match a human observer’s assessments
with an overall accuracy of 76.71%. This classification relies mostly on
detecting features that are generally not desired, such as oscillatory
behaviour or high accelerations. It remains unclear whether this classifier
could detect a shift in sEMG patterns, that results in plausible predictions,
but produces the wrong hand configuration.

The underpinning idea of this work is that automatic failure detection should
instead be able to spot every instance, where the myocontrol output does not
match the user’s intention. Our hypothesis is that we can detect more subtle
shifts in the control mapping by incorporating situational context
information. In the reach-to-grasp scenario, the hand configuration we choose
greatly depends on the shape of the object and how we position our hand
relative to the object (i.e. the local context). From human demonstration we
learn the relation between context and reach-to-grasp trajectories as a set of
generative models. These models are then used as a prior estimate of the user
intention in similar local contexts. By comparing new motions to our models’
predictions, we can assess which of the movement commands do not reflect the
intention and are more likely caused by a myocontrol failure. We developed and
demonstrated our approach in a virtual reality (VR) simulation, where both
patients and able-bodied users can control a model of a prosthetic hand. The
simulation provides a highly controlled test environment, and allows us to
accurately measure, record, evaluate and visualise the user’s movements in the
virtual scene.

This paper is structured as follows. We first describe the implementation of
the VR simulation. We then show how to construct the proposed approach for
learning a situational context model. Finally, we demonstrate that our
approach appropriately distinguishes between sets of successful and failed
control performance, and can be used to detect myocontrol failures.

##  II Method

###  II-A VR Simulation

####  II-A1 Materials

The virtual environment was developed using the Unity 3D game engine and was
displayed on the HTC Vive Virtual Reality Headset. One HTC Vive Tracker was
used for positional tracking of the user’s right forearm. The sEMG data was
gathered using the Myoband by Thalmic Labs, a stretchable bracelet fitted with
8 EMG sensors.

####  II-A2 Myocontrol Implementation

We specifically simulate control of transradial prostheses (i.e. below the
elbow joint). The virtual prosthesis is modelled as a biological hand with 20
degrees of freedom (DoF), that are grouped into two degrees of control (DoC)
for the purpose of this work, namely coupled flexion of the thumb, index
finger and middle finger, and coupled flexion of the ring finger and little
finger. Using the iRR-RFF as in [6], a model is trained to predict two
normalised values, representing the proportional activation of each DoC, from
8-dimensional sEMG patterns. This activation vector is predicted at a
frequency of 200 Hz. At every rendered frame of the simulation, we translate
the latest available activation vector into 20 joint angles using a mapping
matrix. The baseline hand configuration for a zero activation vector
corresponds to the resting position (configuration a) in Fig. 4).

####  II-A3 Object Interaction

Physics engines are unstable when multiple contacts act on an object, such as
in grasping (e.g. [11]). We thus rely on the physics engine only for collision
detection, but allow penetration and specify custom rules for grasp stability.
We classify a grasp as stable if there are at least two contacts with a
minimum angle of 90 degrees between their normals (see Fig. 1). The object is
then simply attached to the wrist’s coordinate frame and moves with it, as
long as the condition is satisfied. This criterion was experimentally tuned to
appear as natural as possible and only fulfills the purpose of enabling basic
interaction with a rigid object.

Figure 1: Example of a stable grasp according to our custom grasp stability
rules in the VR simulation. Two contact points are marked in yellow on the
object surface. The angle between their contact normals is 90 degrees.

###  II-B Failure Detection

We aim to detect potential failures in the myocontrol by learning a set of
generative models that encodes a correct behaviour for attempting a grasp. Our
training data is represented as a set of $N$ states recorded along all the
demonstrated grasping trajectories:

$\displaystyle
D=\left\\{s_{j}:s_{j}\in\mathbb{R}^{3}\times\mathbb{H}\times\left[0,1\right]^{k}\right\\}_{j\in\left[1,N\right]}\,,$
(1)

where each state $s_{j}$ consists of the position $p_{j}$ $\in\mathbb{R}^{3}$
of the wrist in space, the quaternion $q_{j}$ $\in\mathbb{H}$ representing the
orientation of the wrist, and a vector $\mathbf{M}_{j}$
$\in\left[0,1\right]^{k}$ of activation values for $k$ degrees of control that
encode the hand configuration.

The quality of our predictor is evaluated as a normalised weighted mean
squared error (MSE) between a new state $x$ and our model’s prediction based
on the training data $D$. The normalised weighted MSE can be seen as the
expected value of the squared error function of our model. We thus define an
error function $C$ as follows:

$\displaystyle
C\left(x,D\right)=\frac{1}{\sum_{j=1}^{N}w{{}_{j}}}\,\,\sum_{j=1}^{N}\text{MSE}\left(\mathbf{M}_{x},\mathbf{M}_{j}\right)w_{j}\,.$
(2)

The difference in hand configurations is measured by the MSE between two
activation vectors:

$\displaystyle\text{MSE}\left(\mathbf{M}_{x},\mathbf{M}_{j}\right)=\frac{\langle\left(\mathbf{M}_{x}-\mathbf{M}_{j}\right),\left(\mathbf{M}_{x}-\mathbf{M}_{j}\right)\rangle}{k}\,.$
(3)

The angled brackets $\langle\,,\rangle$ are used to denote the inner product
between two vectors. That way, the MSE compares activations for each DoC and
appropriately summarises the total degree of divergence in a single value
between $0$ and $1$, where $0$ expresses perfect congruence.

This error is calculated between $x$ and each state $s_{j}$ in $D$. But for
the overall error of $x$, we want to compare $x$ only to previously
demonstrated states in a similar context, that is states in which the
prosthesis was positioned relative to the object in a similar way. Therefore
we attach a weight $w_{j}$ to each training state, that represents its
relative importance with respect to $x$. This weight is a function of both the
spatial and angular distance between that training state and the current state
$x$:

$\displaystyle w_{j}=w_{\text{position}}\left(p_{x},p_{j}\right)\times
w_{\text{rotation}}\left(q_{x},q_{j}\right)$ (4a)
with $\displaystyle
w_{\text{position}}\left(p_{x},p_{j}\right)=\begin{cases}e^{-\alpha\,\,d\left(p_{x},p_{j}\right)^{2}}&\text{if
}d\left(p_{x},p_{j}\right)\leq\delta\\\ 0&\text{otherwise}\end{cases}$ (4b)
$\displaystyle
w_{\text{rotation}}\left(q_{x},q_{j}\right)=\begin{cases}e^{-\beta\,\,\theta\left(q_{x},q_{j}\right)^{2}}&\text{if
}\theta\left(q_{x},q_{j}\right)\leq\phi\\\ 0&\text{otherwise}\end{cases}$ (4c)
where $d\left(p_{x},p_{j}\right)$ is calculated as the euclidean distance be-
tween the two points, and
$\displaystyle\theta\left(q_{x},q_{j}\right)=\cos^{-1}\left(\frac{\langle
q_{x},q_{j}\rangle}{\sqrt{\left\langle
q_{x},q_{x}\right\rangle}\times\sqrt{\left\langle
q_{j},q_{j}\right\rangle}}\right).$ (4d)

These two weighting terms decrease exponentially with the square distance in
space and in angle respectively. Additionally, training states further than
the cut-off distances $\delta$ and $\phi$ from $x$ are excluded entirely.
Thus, these two parameters define a focal area (or area of similar context)
around $x$, in which training states are considered to be relevant. Parameters
$\alpha$ and $\beta$ can be tuned to specify how much points on the edge of
the focal area contribute to the evaluation relative to points in the centre:

$\displaystyle\alpha=-\frac{\text{ln}(\sqrt{r})}{(\delta)^{2}}$ (5a)
$\displaystyle\beta=-\frac{\text{ln}(\sqrt{r})}{(\phi)^{2}}$ (5b)

where $r$ is the relative contribution ratio or relative importance of
training states on the edge of the focal area (for example see Fig. 2).

$w_{\text{position}}(d)$ / $w_{\text{rotation}}(\theta)$

Figure 2: Weight functions $w_{\text{position}}$ with respect to the lower
x-axis and $w_{\text{rotation}}$ with respect to the upper x-axis, plotted for
the parameter values $r=0.25$, $\delta=0.02$m, $\phi=20^{\circ}$,
$\alpha=1733$ and $\beta=0.001733$. Both functions implement the same decrease
of relative importance for points within their respective cut-off distances
$\delta$ and $\phi$. In this example the values range between $0.5$ and $1$,
so that the combined weight $w_{j}$ can assume values between $0.25$ and $1$.

Then, we integrate all those single state comparisons into one overall error
function $C$ for the state $x$. In Eq. (2) this is computed as the normalised
weighted mean of the differences in hand configurations, weighted by their
similarity in context, rather than a sum or an arithmetic mean over the
comparisons. The weighted mean appropriately expresses the overall divergence
of the hand configuration from previous demonstrations in a similar context.
It is not distorted by the amount of available data and takes into account the
relative importance of training states.

In summary, the error function given in Eq. (2) can be used to calculate a
mismatch between the user’s motion and the model’s prediction in a similar
context. The learned models can be interpreted as an estimate of user
intention, thus the error values also encode how much the performed action
differs from the command we assume the user might have attempted. We can also
infer whether the activation vector in the new state is likely an accurate
reflection of the user’s intention (low values) or rather a failure in the
myocontrol (high values). To do so, a threshold can be defined, above which
states are classified as myocontrol errors.

###  II-C Experimental Method

The ability of the proposed function in Eq. (2) to detect failures was
analysed on a reach-to-grasp task. One set of training data and two sets of
test data (a success condition and a failure condition) were collected in a
single session by the same able-bodied subject, using the VR simulation
described in section II-A.

####  II-C1 Procedure

The sEMG bracelet and the Vive tracker were placed on the user’s forearm as
shown in Fig. 3.

a)b)

Figure 3: Placement of the sensors on the user’s arm: a) The Vive Tracker is
placed just below the wrist joint and secured using adhesive tape, b) The
Myoband is placed approximately 5cm below the elbow joint, where it picks up
electrical signals at 8 positions around the forearm. Exact placement of the
electrodes is not necessary for regression model approaches, as long as they
cover different muscles around the forearm (e.g. [6, 2, 16]).

A myocontrol model was trained on the three different hand configurations
shown in Fig. 4. Just as in [6] and [16], the subject performed the desired
grasps by copying the visual cues in Fig. 4 at the maximal comfortable level
of force. The resulting sEMG activation pattern was labeled with the two-
dimensional activation array corresponding to the displayed grasp. In order to
achieve an accurate reflection of the user’s intention, the model was
carefully adapted incrementally until the subject was satisfied with the
control performance and subjectively perceived no control failures. In total,
each grasp was recorded eight times, four times in a relaxed position with the
right arm held close to the body and bent by 90 degrees at the elbow, and four
times with the arm extended forward and the palm facing inward. Once the
system was trained, the subject put on the head-mounted display and started
the data collection.

For the training data, the subject performed a total of 22 trials of the
grasping task. In each trial a capsule-shaped object appeared on a table and
within reaching distance of the user in the virtual scene. The subject then
reached for the object and picked it up, while the application recorded the
trajectory and configuration of the virtual hand continuously until the first
stable grasp on the object was achieved. The capsule shape was chosen, as it
resembles many real life demands (such as picking up a bottle), which can be
solved in a number of different ways, including at least two different grasp
types (power grasp and tridigital grasp) and several different contact regions
and arm configurations. So instead of repeatedly performing the same motion,
the subject was encouraged to demonstrate a variety of different grasp
solutions on the object. For a visualisation of the training data set, see the
left panel of Fig. 5.

a)b)c)

Figure 4: Hand configurations used for training in the experiment:
a) Resting position b) Power grasp c) Tridigital grasp.

Because this training data serves as a sample of natural grasps a human user
would attempt, it is essential to avoid accidentally recording any myocontrol
failures within this sample. Thus, all indistinct trials, that is trials in
which the myocontrol did not exactly reproduce the subject’s intention, were
discarded immediately, per the judgment of the subject. Twice during this data
collection a single update was incrementally added to the myocontrol model to
counteract performance degradation.

Next, we collected test data for the validation of the proposed error
function. The success condition is a sample of successfully executed grasps,
representing the desired state of a reliable myocontrol system, and containing
no myocontrol failures. For this sample 8 trials were collected in exactly the
same way as described above, where the execution accurately reflects the
subject’s intention. For the failure condition, we aimed to collect data, that
represents the myocontrol failures in the unreliable state of the myocontrol
system following a shift in the sEMG inputs. A second myocontrol model was
trained exactly as before, but this time the labels (i.e. the activation
vectors) for the grasps b) and c) in Fig. 4 were swapped. By intentionally
mis-labelling the grasps, we effectively trained an action to be associated
with a different intention. The resulting activations and model predictions
were plausible, but produced a different hand configuration from what the user
intended. Another 20 trials were recorded with this second myocontrol model.

####  II-C2 Data Analysis

The three data samples were represented as arrays of states (as defined in Eq.
(1)), sampled at approximately every 12 milliseconds along each trajectory. We
then conducted two evaluations of both experimental conditions using the
function specified in Eq. (2). For the first evaluation, we calculated the
value of the function for every single test state in both experimental
conditions based on the full training data; a total of 2789 training states,
shown in the left panel of Fig 5.

Figure 5: Visualisation of the training data. The left panel shows 2789 states
sampled along the whole length of all demonstrated trajectories. The right
panel shows the 440 target area training states, representing only the last
200ms of each trajectory. Trajectories pictured in yellow ended in a
tridigital grasp, red trajectories ended in a power grasp.

The second time, we calculated the errors based on a subset of the training
data, containing only the last 20 recorded data points before the stable grasp
was reached in each demonstrated trial. This smaller target area training set
contains a total of 440 training states representing roughly the last 200ms of
each trajectory, as shown in the right panel of Fig 5. This second evaluation
was conducted, because in qualitative observations, most trajectories in both
conditions started in a very similar way and only diverged into different hand
configurations close to the object. We expected that the difference in error
distributions between the conditions also becomes more pronounced closer to
the object surface. Thus, evaluating the test data only on the smaller but
more determining and relevant target area training sample should suffice to
distinguish between the two conditions. Finally, all test states for which
there were not enough training states $n$ available in a similar context for
an evaluation ($n<5$) were excluded from further analyses. All function
parameters were tuned experimentally and kept fixed throughout these
experiments (Table I). For an overview of the final sample sizes, see Table
II.

TABLE I: Parameter Choices for the Error Function $k$ | $n$ | $r$ | $\delta$ | $\phi$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---
$2$ | $5$ | $0.25$ |  $0.02$m | $20^{\circ}$ | $1733$ | $0.001733$
TABLE II: Overview of test sample sizes based on both positional and orientational information of the wrist. |  Success condition |  Failure condition
---|---|---
Number of grasp trials | 8 | 20
Total number of sampled test states | 1099 | 2066
Number of evaluable test states based on the full training data | 605 | 765
Number of evaluable test states based on the target area training data | 306 | 338

We then compared the sets of calculated error values across the two
conditions. Overall differences in the distribution of errors were analysed
using statistical tests. Furthermore, we closely examined the absolute range
of the values, as this feature in particular determines how appropriate a
threshold would be as a final failure detection criterion. One possible
concern of this approach is, that single trials could be over-represented
within their samples. Because of different movement speeds during the data
collection, the number of states within single trials varies, and we need to
ensure that the overall results in the direct comparison are not distorted,
but that similar results are obtained when the data is grouped into trials for
analysis. Therefore, we also examined the overall distributions, the mean
error and the ranges along each trial. It is important to note however, that
the error is intended to be interpreted with respect to the single current
state, and this grouped analysis only serves as a validation of direct
comparison results. And finally, we tuned an exemplary threshold value to the
results and determined the sensitivity and specificity of the eventual failure
classification of our test data set.

Figure 6: Comparison of error distributions between the two experimental
conditions when both positional and rotational information were available. The
chart shows Tukey boxplots, where the whiskers represent data within 1.5
inter-quartile-ranges of the first and third quartile. The bars illustrate the
absolute range of error values. The two groups to the left display the results
obtained in evaluation 1 and 2, while the the group on the right represents
just one example of how error distributions can be influenced by tuning other
variables of the error function.

To translate this approach into a mobile prosthesis, the device needs to be
fitted with a depth camera and on-board positional and rotational tracking
that does not require external reference (for example inertial tracking, e.g.
[18]). Since this positional information is much harder to obtain than
orientational information, we also assessed the distribution of errors and the
potential failure detection performance based only on the orientation of the
wrist, in case there is no positional information of the wrist available at
all.

##  III Results

###  III-A Based on both positional and rotational information

The distributions of the errors for both conditions in both evaluations are
shown in the left two groups of Fig. 6. Most of the data accumulates around
the low end of the spectrum. Overall, however, the failure condition yielded
higher values than the success condition (Table III). This difference was
further analysed using a permutation test.

TABLE III: Comparison of error distributions based on both positional and
orientational information of the wrist.

Evaluation 1 | Success condition | Failure condition
---|---|---
$M$ | 0.0129 | 0.0437
$SD$ | 0.0269 | 0.1066
Median | 0.0054 | 0.0120
$p$ |  $<$0.001

Evaluation 2 | Success condition | Failure condition
---|---|---
$M$ | 0.0247 | 0.0948
$SD$ | 0.0317 | 0.1489
Median | 0.0104 | 0.0391
$p$ |  $<$ 0.001

Note. Values reported in this table are the mean (M), standard deviation (SD),
and the median of all calculated errors in each experimental condition, as
well as the one-sided p-value of a permutation test between the conditions
after 5000 permutations.

Specifically, we used a MatLab implementation by Ehninger [3], based on the
corrections proposed by Phipson and Smyth [14]. The test detected
statistically highly significant differences in error distributions with p $<$
0.001 after 5000 permutations for both evaluations.

We then examined the same data grouped into trials, and the effect persists.
The mean error values, when corrected for the different amounts of evaluable
states in single trials, describe the same relation between the two
conditions. The mean range along the trials also differed between the
conditions just as expected from the overall comparison between the groups on
both evaluations (Fig. 7). After a closer look at the descriptive data for
each trial, we conclude that there is no reason to believe that the difference
found in the overall comparison misrepresents the underlying data.

Figure 7: Results of the error evaluation for all test states grouped into
trials. The horizontal lines indicate the mean value, corrected for different
numbers of evaluable test states across single trials. The vertical lines
indicate the average range along single trials in the condition.

The range in error values in the failure condition is substantially wider than
in the success condition by a factor of 2.65 in the first evaluation, and 5.15
in the second evaluation. Fig. 6 also shows that the difference between
conditions is more pronounced in the second evaluation (based on only target
area training data) than in the first evaluation (based on all training data).
Accordingly, we observed that the error development along single trials also
shows a greater divergence between the two conditions towards the end of each
reach-to-grasp trajectory (Fig. 8).

Figure 8: Error values along one single exemplary grasp trajectory of each
experimental condition.

We also found that the maximal error value along each trial is reached very
close to the object surface. In the first evaluation, across all grasp trials
the highest error value is reached on average only 9.75 states away from the
final grasp in the success condition, and 3.15 states in the failure
condition. Similarly, in the second evaluation, the highest errors are on
average 14.63 (success) and 4.85 (failure) states away from the grasp. All of
these averages lie within the last 176ms or less of each trial. Additionally,
we found that the reported differences and especially the difference in range
can be further tuned, by varying other parameters of the error function. As an
example, the third group in Fig. 6 displays error distributions for an
evaluation using the target area training data and $n$ = 20.

Using a threshold of $C_{T}=0.0249$ on our test data (evaluated on the target
area training data with $n=5$), the method correctly detected 18 out of 20
failures trials and correctly classified 6 out of 8 correct trials as non-
failures (Sensitivity = 0.9, Specificity = 0.75).

###  III-B Based on rotational information only

TABLE IV: Sample sizes based on orientational information
of the wrist only. |  Success condition |  Failure condition
---|---|---
Number of grasp trials | 8 | 20
Total number of sampled test states | 1099 | 2066
Number of evaluable test states based on the full demonstration data | 1099 | 2066
Number of evaluable test states based on the target area demonstration data | 1099 | 2051

In this second analysis the focal area of similar context was determined only
by the angular distance between a test state and the training data. As a
result, more training states were considered to be relevant to the current
situational context. This greater generalisation increased the number of
evaluable states (Table IV), but at the same time slightly reduced the
divergence between the two experimental conditions (Table V, Fig. 9). However,
there is still a clear distinction between the success and the failure
condition. Using a threshold value of $C_{T}=0.0511$ on our test data
(evaluated on all training data) the method achieved a sensitivity of 0.6, and
a specificity of 0.875.

TABLE V: Comparison of error distributions based on orientational information
of the wrist only.

Evaluation 1 | Success condition | Failure condition
---|---|---
$M$ | 0.0192 | 0.0436
$SD$ | 0.0231 | 0.1022
Median | 0.0189 | 0.0217
$p$ |  $<$0.001

Evaluation 2 | Success condition | Failure condition
---|---|---
$M$ | 0.0533 | 0.0719
$SD$ | 0.0307 | 0.0802
Median | 0.0454 | 0.0546
$p$ |  $<$ 0.001

Note. Values reported in this table are the mean (M), standard deviation (SD),
and the median of all calculated errors in each experimental condition, as
well as the one-sided p-value of a permutation test between the conditions
after 5000 permutations.

##  IV Discussion

The experimental results indicate that the proposed approach can potentially
be used to monitor the status of a prosthetic control system, detect and even
anticipate failures, and therefore increase its accuracy and reliability over
time. The model can clearly distinguish between the two conditions; while
lower values occur in both conditions, for high values it becomes more and
more likely that the observed state is a failure, as there is a large range of
error values that only ever occurred in the failure condition (see Fig. 6).
The way to achieve an optimal failure detection must therefore be to minimise
the overlap of the two conditions, or the area of uncertainty, and tune the
error calculation in such a way, that the data in the two groups diverge as
much as possible.

One simple way of determining whether we are facing a failure is then to use a
threshold for triggering a model update; such a threshold will always
constitute a trade-off between sensitivity and specificity, but the more
distinctly we can isolate the two conditions, the more accurately and
confidently we can place the threshold and identify failures. The exact value
could even be tuned individually, depending on the available data and user
preference.

In this work we have made several assumptions. Firstly, the experiment was
limited to two grasp types only, and secondly, the task only included a single
object shape. Additionally, the approach in this paper was tested only on one
non-disabled expert user. Future work will further explore the proposed method
in a broader experiment with more participants, including patients from the
target population. We also assumed that the user attempts a sensible grasping
motion, as opposed to e.g. punching the object. In such cases the system will
detect a failure, but the user should be able to dismiss the alert.

Figure 9: Comparison of error distributions between the two experimental
conditions, based on only orientational information of the wrist. The chart
shows Tukey boxplots, where the whiskers represent data within 1.5 inter-
quartile-ranges of the first and third quartile. The bars illustrate the
absolute range of error values.

As future extension we will also build on previous work of Kopicki and Zito
[9, 17, 19, 20] who have demonstrated that a set of generative models can be
efficiently learned, for a robot manipulator, in one shot such that
manipulative contacts and trajectories are computed for previously unseen
objects. By integrating these methods in our system, we will replace the human
demonstration samples as a prior estimate of which grasps the user might
attempt, and more importantly it would allow generalisation to new object
shapes.

##  V Conclusion

In this work we introduced an automatic failure detection method that
continuously monitors the system performance. Although we are still in an
early stage of development, our results indicate that the approach performs
well and has potential for future extensions. We believe that this is a
promising first step towards a truly reliable self-correcting control system,
highly interactive and customisable, without the need for excessive training.

## References

  * [1] C. Castellini, G. Passig, and E. Zarka.  Using ultrasound images of the forearm to predict finger positions.  IEEE Trans. on Neural Systems and Rehabilitation Engineering, 20(6):788–797, 2012.
  * [2] Claudio Castellini, Emanuele Gruppioni, Angelo Davalli, and Giulio Sandini.  Fine detection of grasp force and posture by amputees via surface electromyography.  Journal of Physiology, 103(3-5):255–262, 2009.
  * [3] B. Ehinger.  Permtest, 2016.
  * [4] Dario Farina, Ning Jiang, Hubertus Rehbaum, Ales Holobar, Bernhard Graimann, Hans Dietl, and Oskar C. Aszmann.  The extraction of neural information from the surface EMG for the control of upper-limb prostheses: Emerging avenues and challenges.  IEEE Trans. on Neural Systems and Rehabilitation Engineering, 22(4):797–809, 2014.
  * [5] A. Fougner, O. Stavdahl, P. J. Kyberd, Y. G. Losier, and P. A. Parker.  Control of upper limb prostheses: Terminology and proportional myoelectric control—a review.  IEEE Trans. on Neural Systems and Rehabilitation Engineering, 20(5):663–677, 2012.
  * [6] Arjan Gijsberts, Rashida Bohra, David Sierra González, Alexander Werner, Markus Nowak, Barbara Caputo, Maximo A. Roa, and Claudio Castellini.  Stable myoelectric control of a hand prosthesis using non-linear incremental learning.  Frontiers in Neurorobotics, 8, 2014.
  * [7] Janne M. Hahne, Meike A. Schweisfurth, Mario Koppe, and Dario Farina.  Simultaneous control of multiple functions of bionic hand prostheses: Performance and robustness in end users.  Science Robotics, 3(19), 2018.
  * [8] Ning Jiang, Strahinja Dosen, Klaus-Robert Muller, and Dario Farina.  Myoelectric control of artificial limbs—is there a need to change focus? [in the spotlight].  IEEE Signal Processing Magazine, 29(5):152–150, 2012.
  * [9] Marek Kopicki, Renaud Detry, Maxime Adjigble, Rustam Stolkin, Ales Leonardis, and Jeremy L. Wyatt.  One-shot learning and generation of dexterous grasps for novel objects.  The Int’l Journal of Robotics Research, 35(8):959–976, sep 2015\.
  * [10] Silvestro Micera, Jacopo Carpaneto, and Stanisa Raspopovic.  Control of hand prostheses using peripheral information.  IEEE Reviews in Biomedical Engineering, 3:48–68, 2010.
  * [11] Kiran Nasim and Young J. Kim.  Physics-based interactive virtual grasping.  In HCI Korea 2016. The HCI Society of Korea, 2016.
  * [12] Markus Nowak, Sarah Engel, and Claudio Castellini.  A preliminary study towards automatic detection of failures in myocontrol.  In MEC17 - A Sense of What’s to Come, 2017.
  * [13] Bart Peerdeman, Daphne Boere, Heidi Witteveen, Rianne Huis in ’t Veld, Hermie Hermens, Stefano Stramigioli, Hans Rietman, Peter Veltink, and Sarthak Misra.  Myoelectric forearm prostheses: State of the art from a user-centered perspective.  The Journal of Rehabilitation Research and Development, 48(6):719, 2011.
  * [14] Belinda Phipson and Gordon K. Smyth.  Permutation $p$-values should never be zero: Calculating exact $p$-values when permutations are randomly drawn.  Statistical Applications in Genetics and Molecular Biology, 9(1), 2010.
  * [15] David Sierra González and Claudio Castellini.  A realistic implementation of ultrasound imaging as a human-machine interface for upper-limb amputees.  Frontiers in Neurorobotics, 7, 2013.
  * [16] Ilaria Strazzulla, Markus Nowak, Marco Controzzi, Christian Cipriani, and Claudio Castellini.  Online bimanual manipulation using surface electromyography and incremental learning.  IEEE Trans. on Neural Systems and Rehabilitation Engineering, 25(3):227–234, mar 2017.
  * [17] J. Stuber, M. Kopicki, and C. Zito.  Feature-based transfer learning for robotic push manipulation.  In Proceeding of IEEE International Conference on Robotics and Automation (ICRA), 2018.
  * [18] Xiaoping Yun, Eric R. Bachmann, Hyatt Moore, and James Calusdian.  Self-contained position tracking of human movement using small inertial/magnetic sensor modules.  In Proceedings IEEE Int’l Conf. on Robotics and Automation, 2007\.
  * [19] C. Zito, M. Kopicki, R. Stolkin, C. Borst, F. Schmidt, M. A. Roa, and J. L. Wyatt.  Sequential trajectory re-planning with tactile information gain for dextrous grasping under object-pose uncertainty.  In Proceeding of IEEE International Conference on Intelligent Robots and Systems (IROS), pages 2013–2040, 2013.
  * [20] C. Zito, V. Ortenzi, M. Adjigble, M. S. Kopicki, R. Stolkin, and J. L. Wyatt.  Hypothesis-based belief planning for dexterous grasping.  arXiv preprint, arXiv:1903.05517 [cs.RO] (cs.AI), 2019.