footnote

# Linear Computation Coding

Ralf R. Müller1, Bernhard Gäde1, Ali Bereyhi1
This paper was/will be presented in part at the Information Theory &
Applications Workshop, San Diego, CA, Feb. 2020 and the IEEE Conference on
Acoustics, Speech, and Signal Processing, Toronto, Canada, Jun. 2021.This work
was partially supported by Deutsche Forschungsgemeinschaft (DFG) under grant
MU-3735/6-1.  1Institute for Digital Communications, Friedrich-Alexander
Universität Erlangen-Nürnberg, Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We introduce the new concept of computation coding. Similar to how rate-
distortion theory is concerned with the lossy compression of data, computation
coding deals with the lossy computation of functions.

Particularizing to linear functions, we present an algorithm to reduce the
computational cost of multiplying an arbitrary given matrix with an unknown
column vector. The algorithm decomposes the given matrix into the product of
codebook and wiring matrices whose entries are either zero or signed integer
powers of two.

For a typical implementation of deep neural networks, the proposed algorithm
reduces the number of required addition units several times. To achieve the
accuracy of 16-bit signed integer arithmetic for 4k-vectors, no multipliers
and only 1.5 adders per matrix entry are needed.

###### Index Terms:

approximate computing, computational complexity, estimation error, fixed-point
arithmetic, linear systems, rate-distortion theory, quantization.

OAS

    oversampled adaptive sensing
CSI

    channel state information
AWGN

    additive white Gaussian noise
i.i.d.

    independent and identically distributed
r.h.s.

    right hand side
l.h.s.

    left hand side
w.r.t.

    with respect to
RS

    replica symmetry
RSB

    replica symmetry breaking
MSE

    mean squared error
MMSE

    minimum MSE
SINR

    signal to interference and noise ratio
MF

    matched filtering

##  I Motivation

Neural networks are becoming an integral part of modern day’s reality. This
technology consists of two stages: A training phase and an inference phase.
The training phase is computationally expensive and typically outsourced to
cluster or cloud computing. It takes place only now and then, eventually only
once forever. The inference phase is implemented on the device running the
application. It is repeated whenever the neural network is used. This work
solely targets the inference phase after the neural network has been
successfully trained.

The inference phase consists of scalar nonlinearities and matrix-vector
multiplications. The computational burden is overwhelmingly at the latter. The
target of this work is to reduce the computational cost of the following task:
Multiply an arbitrary unknown vector with a known, but arbitrary matrix. At
the first layer of the neural network, the unknown vector is the input to the
neural network. At a subsequent layer, it is the output of the respective
previous layer. The known matrices are the outcome of the training phase in
the data fusion center and fixed for all inference cycles of the neural
network.

The computing unit running the inference phase need not be a general-purpose
processor. With neural networks being more and more frequently deployed in
low-energy devices, it is attractive to employ dedicated hardware. For some of
them, e.g, field programmable gate arrays or application-specific integrated
circuits with memory, the data center has the option to update the known
matrices, whenever it wants to reconfigure the neural network. Still, the
matrices stay fixed for most of the time. In this work, we will not address
those updates, but focus on the most computationally costly effort: The
frequent matrix-vector multiplications within the dedicated hardware.

Besides the matrix-vector multiplications, memory access is currently also
considered a major bottleneck in the inference phase of neural networks.
However, technological solutions to the memory access problem, e.g., stacked
dynamical random access memory utilizing through-silicon vias [1] or emerging
non-volatile memories [2], are being developed and expected to be available
soon. Thus, we will not address memory-access issues in this work.

Various works have addressed the problem of simplifying matrix-matrix
multiplications utilizing certain recursions that result in sub-cubic time-
complexity of matrix-matrix multiplication (and matrix inversion) [3, 4].
However, these algorithms and their more recent improvements, to the best of
our knowledge, do not help for matrix-vector products. This work is not
related to that group of ideas.

Various other works have addressed the problem of simplifying matrix-vector
multiplications in neural networks utilizing structures of the matrices, e.g.,
sparsity [5, 6]. However, this approach comes with severe drawbacks: 1) It
does not allow to design training phase and inference phase independently from
each other. This restricts interoperability, hinders efficient training, and
compromises performance [7]. 2) Sparsity alone does not necessarily reduce
computational cost, as it may require higher accuracy, i.e. larger word-length
for the nonzero matrix elements. In this work, we will neither utilize
structures of the trained matrices nor structures of the input data. The
vector and matrix to be multiplied may be totally arbitrary. They may, but
need not, contain independent identically distributed (IID) random variables,
for instance.

It is not obvious that, without any specific structure in the matrix,
significant computational savings are possible over state-of-the-art methods
implementing matrix-vector multiplications. In the sequel, we will explain why
such savings are possible and how they can be achieved. We also show that
these savings are very significant for typical matrix-sizes in present day’s
deep networks: By means of linear computation coding, the computational cost,
if measured in number of additions and bit shifts, is reduced several times.

The paper is organized as follows: In Section II, the general concept of
computation coding is introduced. In Section III, we review the state-of-the-
art and define a benchmark for comparison. In Section IV, we propose our new
algorithm. Sections V and VI study its performance by analytic and simulative
means, respectively. Section VII summarizes conclusions and gives an outlook
to open problems and applications of linear computation coding beyond the area
of neural networks.

Matrices are denoted by boldface upper letters, vectors are not explicitly
distinguished from scalar variables. The sets $\mathbb{Z}$ and $\mathbb{R}$
denote the integers and reals, respectively. The identity matrix, the all zero
matrix, the all one matrix, the expectation operator, the sign function,
Landau’s big O-operator, and the Frobenius norm are denoted by $\mathbf{I}$,
$\mathbf{0}$, $\mathbf{1}$, $\mathop{\mbox{$\mathsf{E}$}}[\cdot]$,
$\text{sign}(\cdot)$, ${\mathsf{O}}(\cdot)$, and $||\cdot||_{\text{F}}$,
respectively. Indices to constant matrices express their dimensions. The
notation $||\cdot||_{0}$ counts the number of non-zero entries of the vector-
or matrix-valued argument.

##  II Computation Coding for General Functions

The approximation by a deep network is the current state-of-the-art to compute
a multi-dimensional function efficiently. There may be other ones, yet
undiscovered, as well. Thus, we define computation coding for general multi-
dimensional functions. Subsequently, we discuss the practically important case
of linear functions, i.e., matrix-vector products, in greater detail.

The best starting point to understand computation coding is rate-distortion
theory in data compression. In fact, computation coding can be interpreted as
a lossy encoding of functions with a side constraint on the computational cost
of the decoding algorithm. It shares a common principle with lossy source
coding: Random codebooks, if suitably constructed, usually perform well. In
contrast to rate-distortion theory, however, random codebooks turn out useful
even for practical applications of computation coding.

Computation coding consists of computation encoding and computation decoding.
Roughly speaking, computation encoding is to find some approximate
representation $m(x)$ for a given and known function $f(x)$ such that $m(x)$
can be calculated for any arbitrary $x\in{\cal X}$ with low computational cost
and high accuracy. Computation decoding is the calculation of $m(x)$. Formal
definitions are as follows:

######  Definition 1:

Given a probability space $({\cal X},P_{\cal X})$ and a metric $d:{\cal
F}\times{\cal F}\mapsto\mathbb{R}$, a computation encoding with distortion $D$
for given function $f:{\cal X}\mapsto{\cal F}$ is a mapping $m:{\cal
X}\mapsto{\cal M\subseteq F}$ such that
$\mathop{\mbox{$\mathsf{E}$}}_{x\in{\cal X}}[d(f(x),m(x))]\leq D$.

######  Definition 2:

A computation decoding with computational cost $C$ for given operator
$\mathsf{C}$ is an implementation of the mapping $m:{\cal X}\mapsto{\cal M}$
such that $\mathsf{C}[m(x)]\leq C$ for all $x\in{\cal X}$.

The computational cost operator $\mathsf{C}[m(\cdot)]$ measures the cost to
implement the function $m(\cdot)$. It reflects the properties of the hardware
that executes the computation.

Computation coding can be regarded as a generalization of lossy source coding.
If we consider the identity function $f(x)=x$ and the limit $C\to\infty$,
computation coding specializes to lossy source coding with $m(x)$ being the
lossy codeword for $x$. Rate-distortion theory analyzes the trade-off between
distortion $D$ and the code rate. In computation coding, we are interested in
the trade-off between distortion $D$ and computational cost $C$. The code rate
is of no or at most subordinate concern.

The expectation operator in the distortion constraint of Definition 1 is
natural to readers familiar with rate-distortion theory. From a computer
science perspective, it follows the philosophy of approximate computing [8].
Nevertheless, hard constraints on the accuracy of computation can be addressed
via distortion metrics based on the infinity norm, which enforces a maximum
tolerable distortion.

The computational cost operator may also include an expectation. Whether this
is appropriate or not, depends on the goal of the hardware design. If the
purpose is minimum chip area, one usually must be able to deal with the worst
case and an expectation can be inappropriate. Power consumption, on the other
hand, overwhelmingly correlates with average computational cost.

The above definitions shall not be confused with related, but different
definitions in the literature of approximation theory [9]. There, the purpose
is rather to allow for proving theoretical achievability bounds than
evaluating algorithms. The approach to distortion is similar. Complexity,
however, is measured as the growth rate of the number of bits required to
achieve a given upper bound on distortion. This is quite different from the
computational cost in Definition 2.

##  III State of the Art

###  III-A Scalar Linear Functions

The principles and benefits of computation coding are most easily explained
for linear functions. Let us start with scalar linear functions, before we get
to the multi-dimensional case.

######  Example 1:

Let $x\in\mathbb{R}$ have unit variance and $f(x)=tx,t\in\mathbb{R}$. Let the
computation encoding take the form

$\displaystyle m(x)={\rm sign}(t)\sum\limits_{b\in{\cal B}}2^{b}x$ (1)

where the set ${\cal B}\subset\mathbb{Z}$ simply contains the positions of all
the 1-bits in the appropriately rounded binary representation of $t$. Define
the computational cost operator as ${\mathsf{C}}[m(x)]=\max{\cal B}-\min{\cal
B}+1$, independent of $x$. This is the standard way a linear function is
implemented on a modern computer by means of additions and bit shifts.
Considering $t$ as a random variable with uniform distribution on $[-1,+1]$,
the trade-off between average mean-squared distortion and computational cost
is well-known to be [10]

$\displaystyle\mathop{\mbox{$\mathsf{E}$}}\limits_{x\in\mathbb{R}}[f(x)-m(x)]^{2}=4^{-C}/3,$
(2)

i.e., every additional bit of resolution reduces the quantization error by 6
dB.

On average, half of the additions are actually additions of zero. For 16-bit
signed integer arithmetic, we only need to add $(16-1)/2=7.5$ terms, on
average, thus, compute 6.5 additions. The multiplication of a $512\times 4096$
matrix of 16-bit signed fixed-point numbers by an unknown column vector, thus,
requires $512\cdot(4096\cdot 6.5+4095)\approx 15.7$ million additions, i.e.,
7.5 additions per matrix entry, on average. It achieves an average distortion
of $4^{-15}/3$ which is equivalent to $-95$ dB. This is used as benchmark in
the sequel.

For very high precision, i.e. very small distortion, more efficient algorithms
are known from literature that are all based on some version of the Chinese
Remainder Theorem [11, 12, 13, 14, 15, 16]. For precisions relevant in neural
networks, however, these algorithms are not useful.

######  Example 2:

Consider the setting of Example 1, except for the computation encoding to take
the canonical signed digit form [17]

$\displaystyle m(x)=\sum\limits_{(s,b)\in{\cal P}}s2^{b}x$ (3)

with ${\cal P}\subset\\{\pm 1\\}\times\mathbb{Z}$. Like in Example 1, the
number $t$ is approximated by a weighted sum of powers of 2. However, the
weighting coefficients are chosen from the set $\\{-1,0,+1\\}$ instead of
$\\{0,1\\}$. The computational cost operator is defined as
${\mathsf{C}}[m(x)]=|{\cal P}|$. In contrast to Example 1, it does not impose
any cost for additions of zeros, but only counts the number of additions and
subtractions. For the optimum choice of the set ${\cal P}$, the trade-off
between average mean-squared distortion and computational cost is shown in the
appendix to be

$\displaystyle\mathop{\mbox{$\mathsf{E}$}}\limits_{x\in\mathbb{R}}[f(x)-m(x)]^{2}=28^{-C}/3.$
(4)

Thus, every unit increment of the computational cost reduces the quantization
error by 14.5 dB. In this case, the average distortion of $-95$ dB is achieved
for an average of $C=\log_{28}4^{15}\approx 6.24$ additions per matrix entry.
That is a 17% reduction over the benchmark in Example 1.

There are various further improvements to computation coding of scalar linear
functions in the literature, e.g., redundant logarithmic arithmetic [18, 19].
One can also search for repeating bit patterns in the representation of $t$
and reuse the initial computation for later occurrences [20, 21]. One can
represent $t$ by multiple bases at the same time [22]. However, the conversion
into the multiple base number system does not come for free. Area and power
savings were studied in [23] for 32-bit and 64-bit integer arithmetic. Only
for 64-bit arithmetic improvements were reported.

###  III-B Multidimensional Linear Functions

Linear computation coding for products of structured matrices with unknown
real or complex vectors is widely spread in the literature. The most known
example is probably the fast Fourier transform. By contrast, computation
coding for linear functions with unstructured matrices has rarely been studied
at all with the notable exception of [24, 25]. Reference [24] finds up to 40%
improvements for structured matrices, but the performance for random matrices
is reported as only “slightly better”. In [25], the processing on the matrix
entries is optimized utilizing various algorithms including linear programming
and pattern search. Computational cost for random matrices is reported to be
reduced by at most 28% over the scalar canonically signed digit form.

Accelerations to products of unstructured matrices with unknown vectors is
well-studied for Boolean semi-rings [26, 27]. For Boolean semi-rings, [27]
shows that lossless linear computation coding in $K$ dimensions requires at
most ${\mathsf{O}}(K^{2}(\log K)^{-2})$ operations. The mailman algorithm [28]
is inspired by these ideas. It allows to compute the product of a matrix
composed of entries from a finite field with an arbitrary (even real or
complex) vector by at most ${\mathsf{O}}(K^{2}(\log K)^{-1})$ operations. For
this purpose, the target matrix

$\mathbf{T}=\mathbf{B}\mathbf{W}$ (5)

is decomposed into a codebook matrix $\mathbf{B}$ and a wiring matrix
$\mathbf{W}$ in order to simplify the computation of the linear function
$f(x)=\mathbf{T}x$ by $f(x)=\mathbf{B}(\mathbf{W}x)$ with appropriate choices
of the codebook and the wiring matrix.

Let $\mathbf{T}\in{\cal T}^{N\times K}$ for the finite set ${\cal
T}=\\{0,\dots,T-1\\}$ with cardinality $T$. Let $K=T^{N}$ and fix the codebook
matrix $\mathbf{B}\in{\cal T}^{N\times K}$ in such a way that entry
$\mathbf{B}_{nk}$ is the $n^{\text{th}}$ digit of the $T$-ary representation
of $k-1$. Thus, the $k^{\text{th}}$ column of the codebook matrix is the
$T$-ary representation of $k-1$. Since $\mathbf{B}$ contains all the $T^{N}$
possible columns, any column of $\mathbf{T}$ is also a column of $\mathbf{B}$.
The purpose of the wiring matrix is to pick the columns of the codebook matrix
in the right order similar to a mailman who orders the letters suitably prior
to delivery. Since the wiring matrix only reorders the columns of the codebook
matrix, it contains a single one per column while all other elements are zero.
Thus, the product $h=\mathbf{W}x$ does not require any arithmetic. On a
circuit board, it just defines a wiring. For the product $\mathbf{B}h$, [28]
gives a recursive algorithm that only requires ${\mathsf{O}}(K)$ operations.
Decomposing a $K\times K$ matrix into $K/\log_{T}K$ submatrices of size
$\log_{T}K\times K$, the overall complexity scales as ${\mathsf{O}}(K^{2}/\log
K)$.

The drawback of the mailman algorithm is that it can be used only for very low
accuracy of computation. To reach the accuracy of $q$-bit unsigned integer
arithmetic, we need matrices whose size is at least $2^{q}\times 2^{(2^{q})}$.
Arbitrary accuracy can only be achieved for infinite matrix size. Furthermore,
the matrix size must grow exponentially with the desired signal-to-
quantization-noise-ratio.

##  IV Proposed Algorithm

We start with the decomposition of the target matrix into codebook matrix and
wiring matrix as in (5). However, we design at least the wiring matrix in a
manner different from [28]. We also drop the purity of wiring and allow some
computational effort to calculate the product $\mathbf{W}x$.

###  IV-A Single Wiring Matrix

For given target matrix $\mathbf{T}=[t_{1},\dots,t_{K}]\in\mathbb{R}^{N\times
K}$ and codebook matrix $\mathbf{B}\in\\{0,\pm 2^{\mathbb{Z}}\\}^{N\times K}$,
we find the wiring matrix $\mathbf{W}=[w_{1},\dots,w_{K}]$ such that it is a
good solution to the sparse recovery problem

$\mathbf{W}=\operatorname*{\mathrm{argmin}\hskip
2.84526pt}\limits_{\mathbf{\Omega}\in\\{0,\pm 2^{\mathbb{Z}}\\}^{K\times
K}:||\mathbf{\Omega}||_{0}=K+S}||\mathbf{T}-\mathbf{B}\mathbf{\Omega}||_{\text{F}}$
(6)

for some parameter $S$ controlling the computational cost. This wiring matrix
requires $S$ additions and $K+S$ shift operations.

Since the optimum solution to the optimization problem (6) is hard to find in
polynomial time (NP-hard), we propose a greedy algorithm that is demonstrated
in Sections V and VI to perform well. First, we break the matrix optimization
problem into $K$ column-wise optimization problems:

$w_{k}=\operatorname*{\mathrm{argmin}\hskip
2.84526pt}\limits_{\omega\in\\{0,\pm
2^{\mathbb{Z}}\\}^{K}:||\omega||_{0}=1+s}||t_{k}-\mathbf{B}\omega||_{2}\qquad\forall
k$ (7)

with $s=S/K$. Since the column-wise optimization (7) is still NP-hard, we take
the following greedy approach for each column:

  1. 1.

Start with $s=0$ and $\omega=\mathbf{0}_{N\times 1}$.

  2. 2.

Update $\omega$ such that it changes in at most a single component.

  3. 3.

Increment $s$.

  4. 4.

If $s\leq S/K$, go to step 2.

The codebook matrix need not be a binary mailman matrix. Any matrix that can
be multiplied with little effort to the product $h=\mathbf{W}x$ is welcome.
Details are discussed in Section IV-C.

###  IV-B Multiple Wiring Matrices

The wiring matrix can be further decomposed into $L$ sub-wiring matrices

$\mathbf{W}=\mathbf{W}_{1}\mathbf{W}_{2}\cdots\mathbf{W}_{L}.$ (8)

This means that $\mathbf{B}$ serves as codebook for $\mathbf{W}_{1}$ and
$\mathbf{B}\mathbf{W}_{1}\cdots\mathbf{W}_{\ell}$ serves as codebook for
$\mathbf{W}_{\ell+1}$. Repeating the greedy algorithm of Section IV-A for
$\ell$ increasing from $1$ to $L$, the wiring matrices $\mathbf{W}_{\ell}$ are
recursively found.

Multiple wiring matrices are useful, if the codebook matrix $\mathbf{B}$ is
computationally cheap, but poor from a distortion point of view. The product
of a computationally cheap codebook matrix $\mathbf{B}$ with a computationally
cheap wiring matrix $\mathbf{W}_{1}$ can serve as a codebook
$\mathbf{B}\mathbf{W}_{1}$ for subsequent wiring matrices that performs well
with respect to both distortion and computational cost.

Multiple wiring matrices can also be useful if the hardware in the inference
phase favors some serial over fully parallel processing. In this case,
circuitry for multiplying with $\mathbf{W}_{\ell}$ can be reused for
subsequent multiplication with $\mathbf{W}_{\ell-1}$. Note that in the design
phase, wiring matrices are preferably calculated in increasing order of the
index $\ell$, while in the inference phase, they are used in decreasing order
of $\ell$.

###  IV-C Codebook Matrices

For codebook matrices, the computational cost depends on the way they are
designed. Besides being easy to multiply onto a given vector, a codebook
matrix should be designed such that pairs of columns are not collinear. A
column that is collinear to another one is obsolete: It does not help to
reduce the distortion while it requires to compute additions. In an early
conference version of this work [29], we proposed to find the codebook matrix
by sparse quantization of the target matrix. While this results in significant
savings of computational cost over the state of the art, there are even better
designs for codebook matrices. Three of them are detailed in the sequel.

####  IV-C1 Binary Mailman Codebook

In the binary mailman codebook, only the all zero column is obsolete. It is
shown in the appendix that the multiplication of the binary mailman matrix
with an arbitrary vector requires less than $2K$ additions. While performing
well, it lacks flexibility, as it requires the matrix dimensions to fulfill
$K=2^{N}$.

####  IV-C2 Two-Sparse Codebook

We choose the alphabet ${\cal S}\subset\\{0,\pm 2^{0},\pm 2^{1},\pm
2^{2},\dots\\}$ as a subset of the signed positive powers of two augmented by
zero and find $K$ vectors of lengths $N$ such that no pair of vectors is
collinear and each vector has zero norm equal to either 1 or 2. For
sufficiently large size of the subset, those vectors always exist. These
vectors are the columns of the codebook matrix. The ordering is irrelevant. It
turns out useful to restrict the magnitude of the elements of ${\cal S}$ to
the minimum that is required to avoid collinear pairs.

####  IV-C3 Self-Designing Codebook

We set $\mathbf{B}=\mathbf{B}_{0}\mathbf{B}_{1}\mathbf{B}_{2}$ with
$\mathbf{B}_{0}=[\mathbf{I}_{N}\ \mathbf{0}_{N\times(K-N)}]$ and find the
$K\times K$ matrices $\mathbf{B}_{1}$ and $\mathbf{B}_{2}$ recursively via (6)
for $S=K$ interpreting them as wiring matrices for some given auxiliary target
matrix $\tilde{\mathbf{T}}$. The auxiliary target matrix may, but need not be
identical to $\mathbf{T}$. The codebook designs itself taking the auxiliary
target matrix as model.

##  V Performance Analysis

In order to analyze the expected distortion, we resort to IID Gaussian random
codebooks and target vectors, as well as mean-square distortion averaged over
the codebook ensemble. The IID Gaussian codebook is solely chosen, as it
simplifies the performance analysis. In practice, the IID Gaussian random
matrix $\mathbf{B}$ must be replaced by a codebook matrix with low
computational cost, but similar performance. Simulation results in Section VI
will show that practical codebooks perform very similar to IID Gaussian ones.

###  V-A Logarithmic Aspect Ratio

A key point to good performance of the multiplicative matrix decomposition (5)
is the logrithmic aspect ratio. The number of rows of the codebook matrix $N$
scales logarithmically with the number of columns $K$. For a linear
computation code, we define the code rate as

$R=\frac{1}{N}\log_{2}K.$ (9)

The code rate is a design parameter that, as we will see later on, has some
minor impact on the trade-off between distortion and computational cost.

The logarithmic scaling of the aspect ratio is fundamental. This is a
consequence of extreme-value statistics of large-dimensional random vectors:
Consider the correlation coefficients (inner products normalized by their
norms) of $N$-dimensional real random vectors with IID entries in the limit
$N\to\infty$. For any set of those vectors whose size is polynomial in $N$,
the squared maximum of all correlation coefficients converges to zero, as
$N\to\infty$ [30]. Thus, the angle $\alpha$ in Fig. 1

Figure 1:  Decomposition of the approximation error.

becomes a right angle and the norm of the angle error is lower bounded by the
norm of the target vector. However, for an exponentially large set of size
$2^{RN}$ with rate $R>0$, the limit for $N\to\infty$ is strictly positive and
given by rate-distortion theory as $1-4^{-R}$ [31]. The asymptotic squared
relative error of approximating a target vector by an optimal scaling of the
best codeword, thus, is $4^{-R}$. The error itself can be approximated by
another vector of the exponentially large set to get the total squared error
down to $4^{-2R}$. Repeating that procedure for $s$ times, the (squared) error
decays exponentially in $s$. In practice, the scale factor cannot be a real
number, but must be quantized. This additional error is illustrated in Fig. 1
and labelled distance error as opposed to the previously discussed angle
error. As a result, the total squared error does not decay as $4^{-(s+1)R}$,
but with a slightly reduced rate.

The logarithmic aspect ratio has the following impact on the trade-off between
distortion and computational cost: For $K+S$ choices from the codebook, the
wiring matrix has $K+S$ nonzero entries. This implies $S$ additions in the
product $h=\mathbf{W}x$. In return, we get approximated an $N\times K$ target
matrix $\mathbf{T}$ with $N=\frac{1}{R}\log_{2}K=\mathsf{O}(\log K)$ rows. For
any desired distortion $D$, the computational cost of the product
$\mathbf{W}x$ is by a factor $\mathsf{O}(\log K)$ smaller than the number of
entries in the target matrix. This is the same scaling as in the mailman
algorithm. Given such a scaling, the mailman algorithm allows for a fixed
distortion $D$ which depends on the size of the target matrix. The proposed
algorithm, however, can achieve arbitrarily low distortion by setting $S$
appropriately large irrespective of the matrix size. The connection between
the computational cost $C=2K+S$ and distortion $D$ is addressed in the sequel.

###  V-B Angle Error

Consider a unit norm target vector $t\in\mathbb{R}^{N}$ that shall be
approximated by a scaled version of one out of $K$ codewords
$b_{k}\in\mathbb{R}^{N}$ which are random and jointly independent. Denoting
the angle between the target vector $t$ and the codeword $b_{k}$ as
$\alpha_{k}$, we get (norm of) the angle error as

$a_{k}=|\sin\alpha_{k}|.$ (10)

The correlation coefficient between target vector $t$ and codeword $b_{k}$ is
given as

$\rho_{k}=\frac{<t;b_{k}>}{||t||_{2}\,||b_{k}||_{2}}=\cos\alpha_{k}$ (11)

with $<\cdot;\cdot\cdot>$ denoting the inner product. The minimum angle error
$a$ and the correlation coefficients are related by

$a^{2}=1-\max\limits_{k}\rho_{k}^{2}.$ (12)

Next, we will study the behavior of the correlation coefficients in order to
learn about the minimum angle error.

Let ${\text{P}}_{\rho^{2}|t}(r,t)$ denote the cumulative distribution function
(CDF) of the squared correlation coefficient given target vector $t$. As the
columns of the codebook matrix are jointly independent, we conclude that for

$\varrho^{2}=\max_{k}\rho_{k}^{2},$ (13)

we have

${\text{P}}_{\varrho^{2}|t}(r,t)=\left[{\text{P}}_{\rho^{2}|t}(r,t)\right]^{K}.$
(14)

The target vector $t$ follows a unitarily invariant distribution. Thus, the
conditional CDF does not depend on it. In the sequel, we choose $t$ to be the
first unit vector of the coordinate system, without loss of generality.

The squared correlation coefficient $\rho_{k}^{2}$ is known to be distributed
according to the beta distribution with shape parameters $\frac{1}{2}$ and
$\frac{N-1}{2}$ [32, Sec. III.A] and given by

${\text{P}}_{\rho^{2}}(r)={\text{B}}\left(\frac{1}{2},\frac{N-1}{2},r\right)$
(15)

with

${\text{B}}(a,b,x)=\frac{\int\limits_{0}^{x}\xi^{a-1}(1-\xi)^{b-1}{\text{d}}\xi}{\int\limits_{0}^{1}\xi^{a-1}(1-\xi)^{b-1}{\text{d}}\xi}$
(16)

denoting the regularized incomplete Beta function [33]. The distribution of
the squared angle error is, thus, given by

$\displaystyle{\text{P}}_{a^{2}}(r)$
$\displaystyle=1-\left[{\text{B}}\left(\frac{1}{2},\frac{N-1}{2},1-r\right)\right]^{K}.$
(17)

It is shown in the appendix to converge to

$\lim\limits_{K\to\infty}\lim\limits_{N\to\frac{1}{R}\log_{2}K}{\text{P}}_{a^{2}}(r)=\left\\{\begin{array}[]{ll}0&r<4^{-R}\\\
1&r>4^{-R}\end{array}\right..$ (18)

for logarithmic aspect ratios confirming the considerations in Section V-A.
The CDF is depicted in Fig. 2.

Figure 2: CDF of the squared angle error for rate $R=\frac{1}{2}$ and various
numbers of columns $K$.

For finite dimensions, we find the mean squared angle error as

$\displaystyle\overline{a^{2}}$
$\displaystyle=\int\limits_{0}^{1}a\,\text{dP}_{a^{2}}(a)=1-\int\limits_{0}^{1}\text{P}_{a^{2}}(a){\text{d}}a$
(19)
$\displaystyle=\int\limits_{0}^{1}{\text{B}}\left(\frac{1}{2},\frac{N-1}{2},r\right)^{K}{\text{d}}r.$
(20)

using integration by parts.

###  V-C Distance Error

For a target vector of unit norm, the distance error $d$ is bounded from above
by $\frac{1}{3}|\cos\alpha|$. To see this, note that the norm of the optimally
scaled codeword is $|\cos\alpha|$. Furthermore, the maximum error occurs, if
the magnitude of the optimum scale factor $v$ is exactly in the middle of two
adjacent powers of two, say $p$ and $2p$. In this case, we have

$d=|v|-p=2p-|v|\qquad\Rightarrow\qquad|v|-p=\frac{|v|}{3}.$

Assuming the error to be uniformly distributed, it is easily concluded that
the average squared distance error is given by

$\overline{d^{2}}=\frac{1}{27}\,\overline{\varrho^{2}}=\frac{1-\overline{a^{2}}}{27}.$
(21)

Note that the factor $1/27$ slightly differs from the factor $1/28$ in Example
2. Like in Example 2, the number to be quantized is uniformly distributed
within some interval. Here however, the interval boundaries are not signed
powers of two. This leads to a minor increase in the power of the quantization
noise.

###  V-D Total Error

Since distance error and angle error are orthogonal to each other, the average
total squared error is simply given as

$\overline{\epsilon^{2}}=\overline{a^{2}}+\overline{d^{2}}=\frac{1+26\,\overline{a^{2}}}{27}.$
(22)

The average total squared error for a single approximation step is depicted in
Fig. 3 for various rates $R$.

Figure 3: $R^{\text{th}}$ root of the average total squared error vs. the
number of columns $K$ in logarithmic scale for various rates $R$.

Note that the computational cost per matrix entry linearly scales with $R$ for
fixed $K$. In order to have a fair comparison, the average total squared error
is exponentiated with $1/R$. Despite the finite size of the matrices and the
additional distance error due to quantization of the scale factor, it does not
deviate very much from the asymptotic value of $\frac{1}{4}$ found in (18).

Every reduction of error compounds on top of the previous one. If the errors
in $s+1$ subsequent steps of approximation were statistically independent, the
remaining average total squared error would be

$D_{\text{LB}}=\left(\overline{\epsilon^{2}}\right)^{s+1}.$ (23)

for $s+1$ steps. As discussed in the next paragraph, they are not independent
and (23) is just a lower bound.

The angle between the target vector and the best codeword (denoted by
$\alpha$) and the angle between target vector and angle error (denoted by
$\beta$) add to $90^{\circ}$, see Fig. 1. The larger the angle error, the
smaller the angle $\beta$. That means a large angle error implies that the
next target vector is close to the previous one. In other words, a target
vector that is difficult to approximate implies that the next target vector is
likely to be difficult to approximate, as well. On the other hand, a target
vector that is easy to approximate results in a new target vector that is
close to orthogonal to the previous target vector. In this case, it is pretty
much a fair coin flip to decide whether the next approximation is easy or
difficult. Thus, it is more likely to get a difficult target vector than an
easy one, on average. This effect is all the more pronounced the smaller the
size of the codebook $K$ is. This is because smaller values of $K$ increase
the probability of having a large angle error, see Fig. 2.

###  V-E Computational Cost

For sake of simplicity, we use the number of additions to measure
computational cost. Sign changes and shifts are cheaper than additions and
their numbers are also very close to the number of additions.

For wiring matrix $\mathbf{W}_{\ell}$ with $1+s_{\ell}$ nonzero entries per
column, we have to sum $K$ times $1+s_{\ell}$ terms. So, we need $s_{\ell}K$
additions in total for wiring matrix $\mathbf{W}_{\ell}$. All three codebook
matrices discussed in Section IV-C require at most $2K$ additions.

Adding the computational costs of wiring and codebook matrices, we get

$\displaystyle C=2K+\sum\limits_{\ell=1}^{L}s_{\ell}K=(s+2)K$ (24)

with $s=\sum_{\ell=1}^{L}s_{\ell}$. Normalizing to the number of elements of
the $N\times K$ target matrix $\mathbf{T}$, we have

$\displaystyle\tilde{C}=\frac{C}{KN}=\frac{s+2}{N}.$ (25)

The computational cost per matrix entry vanishes with increasing matrix size.
This behavior fundamentally differs from state-of-the-art methods discussed in
Examples 1 and 2, where the matrix size has no impact on the computational
cost per matrix entry.

##  VI Simulation Results

###  VI-A Accuracy of Lower Bound

Consider $L$ wiring matrices each with two nonzero entries per column, i.e.,
$s_{\ell}=1$ for all $\ell$. The total average distortion for IID Gaussian
codebook matrices is shown in Fig. 4

Figure 4: Average relative total squared error vs. total number of additions
per column required for the wiring matrices. The solid lines refer to the
lower bound (23). The markers to simulation results. All results are for IID
Gaussian codebook matrices with unit rate. Results are averaged over $10^{6}$
random matrix entries, except for the circular markers which refer to $10^{7}$
random matrix entries.

vs. $s$, i.e., the total number of additions per column required for the
wiring matrices. Results are averaged over $10^{6}$ random matrix entries. For
unit rate, the lower bound (23) is very tight for $N>6$ and otherwise poor
over a wide range of $s$. Remarkably, the jittery behavior for $N=6$ seems not
to result from insufficient averaging, as going for ten times more matrix
samples does not reduce it.

###  VI-B Codebook Comparison

For an IID Gaussian target matrix with $K=256$ columns, we compare various
codebook matrices in Fig. 5. The wiring matrices are designed as described in
Section VI-A.

Figure 5: Average relative total squared error vs. total number of additions
per column required for the wiring matrices and Gaussian IID target matrices
of size $N\times 256$. Results are averaged over $10^{6}$ random matrix
entries.

For all numbers of rows $N$ that are shown in the figure, self-designing
codebooks perform best while two-sparse codebooks perform worst. The mailman
codebook, which only exists for $N=8$, performs slightly worse than the IID
Gaussian codebook. The choice of the codebook does not affect the slope of the
curve, but is responsible for an offset shift in the average distortion.

###  VI-C Number of Additions per Matrix Entry

Due to their superior performance, we restrict the subsequent considerations
to self-designing codebooks. First, we address IID Gaussian target matrices.
Here, the target matrix can be used to self-design the codebook matrix, as IID
Gaussian codebooks work well. As in the previous subsections, we use multiple
wiring matrices with $s_{\ell}=1$ for all of them. Table I

TABLE I: Number of additions per matrix entry that are required to reach the accuracy of 2-, 4-, 8-, 16-, and 24-bit signed integer arithmetic or better. Best values for given matrix width shown in boldface.  $\tilde{C}$ | 2 bit | 4 bit | 8 bit | 16 bit | 24 bit
---|---|---|---|---|---
$4\times 256$ | 0.5 | 0.75 | 1.25 | 2.25 | 3.25
$5\times 256$ | 0.4 | 0.8 | 1.2 | 2.2 | 3.2
$6\times 256$ | 0.333 | 0.667 | 1.167 | 2.333 | 3.333
$7\times 256$ | 0.427 | 0.714 | 1.143 | 2.286 | 3.286
$8\times 256$ | 0.375 | 0.625 | 1.25 | 2.25 | 3.375
$9\times 256$ | 0.333 | 0.667 | 1.222 | 2.333 | 3.444
$10\times 256$ | 0.3 | 0.6 | 1.2 | 2.3 | 3.5
$11\times 256$ | 0.364 | 0.636 | 1.182 | 2.364 | 3.546
$12\times 256$ | 0.333 | 0.583 | 1.25 | 2.417 | 3.583
$13\times 256$ | 0.307 | 0.615 | 1.231 | 2.462 | 3.615
$14\times 256$ | 0.286 | 0.643 | 1.214 | 2.5 | 3.714
$15\times 256$ | 0.333 | 0.6 | 1.267 | 2.533 | 3.8
$5\times 1024$ | 0.4 | 0.6 | 1 | 2 | 2.8
$6\times 1024$ | 0.333 | 0.667 | 1 | 1.833 | 2.667
$7\times 1024$ | 0.429 | 0.571 | 1 | 1.857 | 2.714
$8\times 1024$ | 0.375 | 0.625 | 1 | 1.875 | 2.75
$9\times 1024$ | 0.333 | 0.556 | 1 | 1.889 | 2.778
$10\times 1024$ | 0.3 | 0.5 | 1 | 1.9 | 2.7
$11\times 1024$ | 0.273 | 0.546 | 1 | 1.818 | 2.727
$12\times 1024$ | 0.333 | 0.5 | 1 | 1.833 | 2.75
$13\times 1024$ | 0.308 | 0.539 | 1 | 1.846 | 2.769
$14\times 1024$ | 0.286 | 0.5 | 1 | 1.857 | 2.786
$15\times 1024$ | 0.267 | 0.533 | 1 | 1.867 | 2.8
$16\times 1024$ | 0.25 | 0.5 | 1 | 1.875 | 2.813
$17\times 1024$ | 0.294 | 0.529 | 1 | 1.882 | 2.824
$7\times 4096$ | 0.429 | 0.571 | 0.857 | 1.571 | 2.286
$8\times 4096$ | 0.375 | 0.5 | 0.875 | 1.625 | 2.25
$9\times 4096$ | 0.333 | 0.556 | 0.889 | 1.556 | 2.333
$10\times 4096$ | 0.3 | 0.5 | 0.9 | 1.6 | 2.3
$11\times 4096$ | 0.273 | 0.455 | 0.818 | 1.546 | 2.273
$12\times 4096$ | 0.25 | 0.5 | 0.833 | 1.583 | 2.333
$13\times 4096$ | 0.308 | 0.462 | 0.846 | 1.539 | 2.308
$14\times 4096$ | 0.286 | 0.427 | 0.857 | 1.571 | 2.286
$15\times 4096$ | 0.267 | 0.467 | 0.8 | 1.533 | 2.333
$16\times 4096$ | 0.25 | 0.438 | 0.813 | 1.563 | 2.313
$17\times 4096$ | 0.294 | 0.412 | 0.824 | 1.588 | 2.294
$18\times 4096$ | 0.222 | 0.444 | 0.833 | 1.556 | 2.333
$19\times 4096$ | 0.263 | 0.421 | 0.790 | 1.579 | 2.316
$20\times 4096$ | 0.25 | 0.45 | 0.8 | 1.55 | 2.35

shows the average number of additions per matrix entry to achieve at least the
accuracy of standard signed integer arithmetic evaluated by (2). For 8-bit
accuracy, the code rate has little, for $K=1024$ even no impact. For larger
and lower accuracies, higher and lower code rates are favored, respectively.
This trend is slightly distorted by quantization effects, but clear enough to
be spotted. In comparison to the benchmark defined in Section III, the
computational cost is reduced by 80 %.

This behavior is explained as follows: There is the fixed computational cost
of $2K$ additions for the codebook matrix. This cost is shared by the fewer
rows, the larger is the code rate. This increases the computational burden per
row (and also per matrix entry) for larger code rates. For high accuracy, the
wiring matrices clearly dominate the computational cost, so the computation of
the codebook is secondary. Thus, higher rates are favored, as they generally
result in lower distortions, see Fig. 3. For low accuracy, only few wiring
matrices are needed and the relative cost of the codebook is more important.
This shifts the optimum rate towards lower values.

Consider now a random matrix whose entries are uniform IID within $[0,1)$. If
the same algorithm is used as before, the performance degrades significantly,
as such a matrix is not a good codebook. Since all entries are positive, all
codewords are constrained to a single orthant of the full space. In order to
circumvent this problem, we choose an auxiliary target matrix
$\tilde{\mathbf{T}}$ with Gaussian IID entries to self-design the codebook. In
order to keep that codebook for all subsequent wiring, we use only a single
wiring matrix and adapt the parameter $s_{1}$ in such a way as to reach the
desired accuracy. The results are shown in Table II.

TABLE II: Number of additions per matrix entry that are required to reach the accuracy of 2-, 4-, 8-, 16-, and 24-bit signed integer arithmetic or better.  $\tilde{C}$ | 2 bit | 4 bit | 8 bit | 16 bit | 24 bit
---|---|---|---|---|---
$8\times 256$ | 0.5 | 0.75 | 1.25 | 2.375 | 3.5
$10\times 1024$ | 0.4 | 0.6 | 1.1 | 2 | 2.9
$12\times 4096$ | 0.333 | 0.5 | 0.917 | 1.75 | 2.5

A comparison to Table I shows that this procedure only causes a minor
degradation, if any at all.

##  VII Conclusions & Outlook

Linear computation coding by means of codebook and wiring matrices is a
powerful tool to reduce the computational cost of matrix-vector
multiplications in deep neural networks. The number of additions to achieve
the accuracy of $q$-bit signed integer arithmetic is roughly given by

$\displaystyle\frac{2+q}{\log_{2}K}.$ (26)

This means no multiplication and only a single addition for 8-bit integer
arithmetic and vectors with at least $K=1000$ dimensions.

The idea of computation coding is not restricted to linear functions. Its
direct application to multidimensional nonlinear functions promises even
greater reductions in the computational cost of neural networks. We conjecture
that neural networks need neither have precise weights nor be densely
connected, if they are sufficiently deep. As in computation coding of linear
functions, any quantization errors at intermediate layers and any level of
sparsity may be compensated for with more layers and larger dimensions.

Fast matrix-vector products are also important for applications in other areas
of signal processing, e.g., beamforming for wireless multiple-input multiple-
output systems [34], compressive sensing, numerical solvers for partial
differential equations, etc. This opens up for many future research directions
based on linear computation coding.

## Acknowledgement

The authors like to thank Veniamin Morgenshtern, Marc Reichenbach, and Hermann
Schulz-Baldes for helpful discussions.

## Appendix

### Average Distortion of Canonically Signed Digit Form

For the canonically signed digit form, the set of reconstruction values is
${\cal R}=\\{\pm 2^{k}:k\in\mathbb{Z}\\}$. Let the target variable $t$ be
uniformly distributed within $[0,1]$. Thus, only reconstruction values for
$k\leq 0$ are used.

Consider the interval $[2^{k-1},2^{k}]$ with $k\leq 0$. The probability
$p_{k}$ that the variable $t$ falls within that interval is equal to its width
$w_{k}=2^{k-1}$. The average mean-squared quantization error within that
interval is well known to be $w_{k}^{2}/12$. Averaging over all intervals, the
average mean-squared distortion for quantization with a single digit is given
by

$\displaystyle\frac{1}{12}\sum\limits_{k=-\infty}^{0}p_{k}w_{k}^{2}$
$\displaystyle=\frac{1}{12}\sum\limits_{k=-\infty}^{0}w_{k}^{3}=\frac{1}{12}\sum\limits_{k=-\infty}^{0}2^{3k-3}$
(27)
$\displaystyle=\frac{1}{96}\sum\limits_{k=0}^{\infty}8^{-k}=\frac{1}{96}\cdot\frac{8}{7}=\frac{1}{3\cdot
28}.$ (28)

By symmetry, the same considerations hold true, if $t$ is uniformly
distributed within $[-1,0]$. Due to the signed digit representation, these
considerations even extend to $t$ uniformly distributed within $[-1,+1]$.

For representation by zero digits, $t\in[-1,+1]$ is quantized to 0 with
average mean-squared distortion equal to $\frac{1}{3}$. This is 28 times
larger than for quantization with one digit.

The canonically signed digit representation is invariant to any scaling by a
power of two. Thus, any further digit of representation also reduces the
average mean-squared distortion by a factor of 28.

### Additions Required for Binary Mailman Codebook

Let $\mathbf{B}_{N\times K}$ denote the binary mailman matrix with $N$ rows
and $K$ columns. It can be decomposed recursively as

$\displaystyle\mathbf{B}_{N\times
K}=\left[\begin{array}[]{cc}\mathbf{B}_{(N-1)\times\frac{K}{2}}&\mathbf{B}_{(N-1)\times\frac{K}{2}}\\\
\mathbf{0}_{1\times\frac{K}{2}}&\mathbf{1}_{1\times\frac{K}{2}}\end{array}\right].$
(31)

Following [28], we decompose $h$ into its first half $h_{1}$ and its second
half $h_{2}$ such that we get

$\displaystyle\mathbf{B}_{N\times
K}h=\left[\begin{array}[]{c}\mathbf{B}_{(N-1)\times\frac{K}{2}}(h_{1}+h_{2})\\\
\mathbf{1}_{1\times\frac{K}{2}}h_{2}\end{array}\right].$ (34)

Let $c(N)$ denote the number of additions to compute the product
$\mathbf{B}_{N\times K}h$. The recursion (34), implies

$\displaystyle c(N)$ $\displaystyle=\frac{K}{2}+c(N-1)+\frac{K}{2}-1$ (35)
$\displaystyle<c(N-1)+2^{N}.$ (36)

These are $\frac{K}{2}$ additions for $h_{1}+h_{2}$, $c(N-1)$ additions for
the matrix-vector product, and $\frac{K}{2}-1$ additions to sum the components
of $h_{2}$. Note that $\mathbf{B}_{1\times 2}h=h_{2}$. Thus, $c(1)=0$ and
$c(N)<2^{N+1}=2K$.

### Asymptotic Cumulative Distribution Function

In order to show the convergence of the CDF of the squared angle error to the
unit step function, recall the following limit holding for any positive $x$
and $u$

$\lim\limits_{K\to\infty}\left(1-\frac{x}{K^{u}}\right)^{K}=\left\\{\begin{array}[]{ll}0&u<1\\\
\exp(x)&u=1\\\ 1&u>1\end{array}\right..$ (37)

The limiting behavior of ${\text{P}}_{a^{2}}(r)$ is, thus, decided by the
scaling of $1-{\text{P}}_{\rho^{2}}(1-r)$ with respect to $K$. The critical
scaling is $\frac{1}{K}$. Such a scaling implies a slope of $-1$ in doubly
logarithmic scale. Thus,

$\lim\limits_{N\to\infty}\frac{\partial}{\partial(NR)}\log_{2}\left[1-{\text{P}}_{\rho^{2}}(1-r)\right]=-1.$
(38)

Explicit calculation of the derivative yields

$\lim\limits_{N\to\infty}\frac{\int\limits^{1}_{1-r}(1-\xi)^{\frac{N-3}{2}}\xi^{-\frac{1}{2}}\log_{2}\left(1-\xi\right){\text{d}}\xi.}{2R\int\limits^{1}_{1-r}(1-\xi)^{\frac{N-3}{2}}\xi^{-\frac{1}{2}}{\text{d}}\xi.}=-1$
(39)

and saddle point integration gives [35, Chapter 4]

$\frac{1}{2R}\log_{2}[1-(1-r)]=-1.$ (40)

This immediately leads to $r=4^{-R}$.

## References

  * [1] W.-W. Shen, Y.-M. Lin, S.-C. Chen, H.-H. Chang, C.-C. Lin, Y.-F. Chou, D.-M. Kwai, and K.-N. Chen, “3-D stacked technology of DRAM-logic controller using through-silicon via (TSV),” IEEE Journal of the Electron Devices Society, vol. 6, pp. 396–402, Mar. 2018.
  * [2] S. Hong, O. Auciello, and D. W. (Eds.), Emerging Non-Volatile Memories. Springer-Verlag, 2014.
  * [3] V. Strassen, “Gaussian elimination is not optimal,” Numerische Mathematik, vol. 13, pp. 354–356, 1969.
  * [4] D. Copperfield and S. Winograd, “Matrix multiplication via arithmetic progressions,” Journal of Symbolic Computation, vol. 9, no. 3, pp. 251–280, Mar. 1990.
  * [5] S. Han, J. Pool, J. Tran, and W. Dally, “Learning both weights and connections for efficient neural network,” in Proceedings of the Advances in Neural Information Processing Systems, Montreal, Canada, Dec. 2015.
  * [6] C. Louizos, K. Ullrich, and M. Welling, “Bayesian compression for deep learning,” in Proceedings of the Advances in Neural Information Processing Systems, Long Beach, CA, Dec. 2017.
  * [7] U. Evci, F. Pedregosa, A. Gomez, and E. Elsen, “The difficulty of training sparse neural networks,” 2019, arXiv:1906.10732v2.
  * [8] K. Palem and A. Lingamneni, “Ten years of building broken chips: The physics and engineering of inexact computing,” ACM Transactions on Embedded Computing Systems, vol. 12, no. 2s, pp. 87:1–87:23, May 2013.
  * [9] R. A. DeVore, “Nonlinear approximation,” Acta Numerica, vol. 7, pp. 51–150, Jan. 1998.
  * [10] R. M. Gray and D. L. Neuhoff, “Quantization,” IEEE Transactions on Information Theory, vol. 44, no. 6, pp. 2325–2383, Oct. 1998.
  * [11] A. Karatsuba and Y. Ofman, “Multiplication of multidigit numbers on automata (in Russian),” Doklady Akademii Nauk SSSR, vol. 145, no. 2, pp. 293–294, 1962.
  * [12] A. L. Toom, “The complexity of a scheme of functional elements simulating the multiplication of integers (in Russian),” Doklady Akademii Nauk SSSR, vol. 150, no. 3, pp. 496–498, 1963.
  * [13] S. A. Cook and S. O. Aanderaa, “On the minimum computation time of functions,” Transactions of the American Mathematical Society, vol. 142, pp. 291–314, Aug. 1969.
  * [14] A. Schönhage and V. Strassen, “Schnelle Multiplikation großer Zahlen,” Computing, vol. 7, pp. 281–292, Sep. 1971.
  * [15] M. Fürer, “Faster integer multiplication,” SIAM Journal on Computing, vol. 39, no. 3, pp. 979–1005, 2009.
  * [16] D. Harvey, J. van der Hoeven, and G. Lecerf, “Even faster integer multiplication,” Journal of Complexity, vol. 36, pp. 1–30, Oct. 2016\.
  * [17] A. D. Booth, “A signed binary multiplication technique,” The Quarterly Journal of Mechanics and Applied Mathematics, vol. 4, no. 2, pp. 236–240, Jan. 1951.
  * [18] M. G. Arnold, T. A. Bailey, J. R. Cowles, and J. J. Cupal, “Redundant logarithmic arithmetic,” IEEE Transactions on Computers, vol. 39, no. 8, pp. 1077–1086, Aug. 1990.
  * [19] H. Huang, M. Itoh, and T. Yatagai, “Modified signed-digit arithmetic based on redundant bit representation,” Applied Optics, vol. 33, no. 26, pp. 6146–6156, Sep. 1994.
  * [20] R. I. Hartley, “Subexpression sharing in filters using canonic signed digit multipliers,” IEEE Transactions on Circuits and Systems-II: Analog and Digital Signal Processing, vol. 43, no. 10, pp. 677–688, Oct. 1996.
  * [21] V. Lefèvre, Moyens arithmétiques pour un calcul fiable. École Normale Supérieure de Lyon: Ph. D. thesis, Jan. 2000, see also INRIA Research Report no 4192, May 2001.
  * [22] V. S. Dimitrov, G. A. Jullien, and W. C. Miller, “Theory and applications of the double-base number system,” IEEE Transactions on Computers, vol. 48, no. 10, pp. 1098–1106, Oct. 1999.
  * [23] V. S. Dimitrov, K. U. Järvinen, and J. Adikari, “Area-efficient multipliers based on multiple-radix representations,” IEEE Transactions on Computers, vol. 60, no. 2, pp. 189–201, Feb. 2011.
  * [24] N. Boullis and A. Tisserand, “Some optimizations of hardware multiplication by constant matrices,” IEEE Transactions on Computers, vol. 54, no. 10, pp. 1271–1282, Oct. 2005.
  * [25] L. Aksoy, P. Flores, and J. Monteiro, “A novel method for the approximation of multiplierless constant matrix vector multiplication,” EURASIP Journal on Embedded Systems, May 2016.
  * [26] M. A. Kronrod, V. L. Arlazarov, E. A. Dinic, and I. A. Faradzev, “On economic construction of the transitive closure of a direct graph (in Russian),” Doklady Akademii Nauk SSSR, vol. 194, no. 3, pp. 487–488, 1970.
  * [27] R. Williams, “Matrix-vector multiplication in sub-quadratic time (some preprocessing required),” in Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 2007, pp. 995–1001.
  * [28] E. Liberty and S. W. Zucker, “The mailman algorithm: A note on matrix-vector multiplication,” Information Processing Letters, vol. 109, pp. 179–182, Jan. 2009.
  * [29] R. Müller, B. Gäde, and A. Bereyhi, “Efficient matrix multiplication: The sparse power-of-2 factorization,” in Proc. of Information Theory & Applications Workshop, San Diego, CA, Feb. 2020, https://arxiv.org/abs/2002.04002v2.
  * [30] T. Jiang, “The asymptotic distribution of the largest entries of sample correlation matrices,” The Annals of Applied Probability, vol. 14, no. 2, pp. 865–880, 2004.
  * [31] T. Berger, Rate Distortion Theory. Eaglewood Cliffs, NJ: Prentice-Hall, 1971.
  * [32] R. Müller, “On approximation, bounding & exact calculation of block error probability for random code ensembles,” IEEE Transactions on Communications, 2021, early access.
  * [33] J. Spanier and K. B. Oldham, An Atlas of Functions. Berlin, Germany: Springer-Verlag, 1987.
  * [34] O. Castañeda, S. Jacobsson, G. Durisi, T. Goldstein, and C. Studer, “Finite-alphabet MMSE equalization for all-digital massive MU-MIMO mmWave communication,” IEEE Journal on Selected Areas in Communications, vol. 38, no. 9, pp. 2128–2141, Sep. 2020.
  * [35] N. Merhav et al., “Statistical physics and information theory,” Foundations and Trends® in Communications and Information Theory, vol. 6, no. 1–2, pp. 1–212, 2010.