###### Abstract

In this paper a variation of the classic vector quantization problem is
considered.  In the standard formulation, a quantizer is designed to minimize
the distortion between input and output when the number of reconstruction
points is fixed. We consider, instead, the scenario in which the number of
comparators used in quantization is fixed. More precisely, we study the case
in which a vector quantizer of dimension $d$ is comprised of $k$ comparators,
each receiving a linear combination of the inputs and producing the output
value one/zero if this linear combination is above/below a certain threshold.
In reconstruction, the comparators’ output is mapped to a reconstruction
point, chosen so as to minimize a chosen distortion measure between the
quantizer input and its reconstruction. The Comparison-Limited Vector
Quantization (CLVQ) problem is then defined as the problem of optimally
designing the configuration of the compactors and the choice of reconstruction
points so as to minimize the given distortion. In this paper, we design a
numerical optimization algorithm for the CLVQ problem. This algorithm
leverages combinatorial geometrical notions to describe the hyperplane
arrangement induced by the configuration of the comparators. It also relies on
a genetic genetic meta heuristic to improve the selection of the quantizer
initialization and avoid local minima encountered during optimization. We
numerically evaluate the performance of our algorithm in the case of input
distributions following uniform and Gaussian i.i.d. sources to be compressed
under quadratic distortion and compare it to the classic Linde–Buzo–Gray (LBG)
algorithm.

###### keywords:

Constrained quantization; Source compression; Hyperplane arrangement; LBG
algorithm; Gaussian quantizer; Uniform quantizer.

xx 1 5 Received: date; Accepted: date; Published: date Comparison-limited
Vector Quantization Joseph Chataignon 1, Stefano Rini 2 Joseph Chataignon,
Stefano Rini

##  1 Introduction

Quantization represents one of the crucial operations in signal processing as
it allows for the mapping of analog signals into discrete values, i.e. A2D
conversion. Being such a fundamental processing operation, quantization has
been an important subject of research and numerous results have been derived
for problem. In this work we consider a variation of the classic quantization
setting which appears to be so far neglected in the literature.  Vector
quantizers are typically produced using operational amplifier (op-amp)
comparators that take a signal and a bias as inputs and produce a zero/one
voltage as output, depending on whether the signal is larger/smaller than the
bias. Generally speaking, op-amp comparators have high power consumption and
manufacturing cost. Additionally, generating a precise set of bias values for
each comparator gives rise to a rather complex circuit structure.  For this
reason, in scenarios in which complexity, fabrication cost, or energy
consumption are severely restricted, it is desirable to evaluate the quantizer
performance in terms of the number of comparators it requires to produce its
output. This is in contrast with the classic approach, in which the quantizer
performance is optimized under a constraint on the output cardinality, i.e.
the number of bits used to represent the output in A2D conversion.

As a concrete example of the difference between the classic setting and the
setting we propose, let us consider the case of a two dimensional quantizer
($d=2$) in which each dimension is quantized with a rate of $1.5$ bits-per-
sample ($R=1.5\ \rm bps$), so that In the classic setting, the quantizer is
designed under the constraint that each pair of inputs is represented as one
of three possible reconstruction points. Under this constraint, a quantizer
can be represented as in  Fig. 1: here each of the three reconstruction points
is separated by the others through a edge in the Voronoi region, so that
quantizer design requires three comparisons to distinguish the three
reconstruction points. Accordingly, $k=3$ comparisons are required to
implement this quantization rule. More generally, since every two
reconstruction points are separated by an edge of the Voronoi region, a vector
quantizer might require up to $2^{R}(2^{R}-1)/2$ comparators. This scaling of
the number of comparators with the output cardinality is generally valid only
for low-rate regime, since the number of neighbors of each reconstruction
point is large. Given that this quantization regime is naturally associated
with low-cost and low-complexity devices, this cost scaling is particularly
disadvantageous. A question that naturally arises is whether a better scaling
of the number of comparators with the output cardinality can be attained. To
address this question, note that number of comparators, $k$, equals the total
number of hyperplane segments in the Voronoi regions. In this view, the best
scaling is attained when the $k$ hyperplanes induce the largest number of
partitions of the space of dimension $d$. After some rather straightforward
geometrical consideration, one realizes that for the case of $d=2$ and
$R=1.5$, the largest number of partition is $7$. This corresponds to the
quantizer in Fig. 1. It is now apparent that there exists a large gap in the
optimal quantizer design whether one considers a constraint on the number of
points used in reconstruction or on the number of comparators employed.  It is
indeed the aim of this paper to develop a first understanding of this problem,
which we term the “Comparison-Limited Vector Quantization (CLVQ) problem”.

Figure 1: Schematic representation of the vector quantization in which two
input symbol are mapped into reconstruction points. In (a), three
reconstruction points are obtained using three comparators while, in (b) seven
reconstruction points are obtained from three comparators

###  1.1 Relevant Results

Vector Quantization (VQ) Gray (1984); Gersho and Gray (2012) is a topic of
vast interest, so that we shall briefly review only the literature which is
most closely related to the topic of this publication. Suffices to say that VQ
is widely adopted in speech coding Makhoul et al. (1985); Paliwal and Atal
(1993); Yahampath and Rondeau (2007), image coding Gersho and Ramamurthi
(1982); Nasrabadi and King (1988); Chiranjeevi and Jena (2018), and video
coding Lee and Woods (1995); Regunathan et al. (2014). VQ has been
successfully used for speaker identification Soong et al. (1987); Hasan et al.
(2004), digital watermarking Wu and Chang (2005); Lu et al. (2005), and
clustering Equitz (1989); Lughofer (2008).

In the following, we follow the information theoretical study of the quantizer
performance as described by rate-distortion theory Berger (2003), originally
developed by Shannon in his seminal work Shannon (1948). As originally
formulated, the _rate distortion function_ is the function describing the
minimum rate at which a source of information can be compressed so that the
reconstruction is to within a prescribed distortion from the original value,
in the limit of infinitely many source samples are compressed.  This approach
is concerned with the compression rate needed to attain a prescribed
distortion level, that is the number of bits required on average to represent
a source sample after quantization. The idea of accounting for the number of
comparators required for A2D conversion emerges from the literature of mmWave
MIMO communication Mezghani and Nossek (2007), in which high-resolution
quantization is no longer feasible, duo to the high energy consumption in
high-speed A2D conversion Sundstrom et al. (2008). For this reason, authors
have investigated the transmission performance achievable when quantization
only preserves the sign of the channel output Mezghani and Nossek (2007); Mo
and Heath (2014). Generalizing the one-bit quantization case, in Rini et al.
(2017) the authors investigate the capacity of a MIMO channel with output
quantization constraint, i.e. the MIMO channel in which the channel output is
processed at the receiver using a finite number of one-bit threshold
quantizers. The channel model of Rini et al. (2017) is relevant in mm-wave
communications which allows for a large number of receiver antennas, while the
number of A2D conversion modules remains small due to limitation in the energy
and costs of radio-frequency modules. Building on an idea of Mo and Heath
(2014), a connection between combinatorial geometry and the MIMO channel with
output quantization constraints of Rini et al. (2017) is drawn in Khalili et
al. (2018). In particular, in Khalili et al. (2018) it is shown that each
quantizer can be interpreted as a hyperplane bisecting the transmitter signal
space; for this reason the largest rates are attained by the configuration
allowing for the largest number of partitions induced by the set of
hyperplanes. To the best of our knowledge, the problem of hardware limited
quantization has so far only been considered in Shlezinger et al. (2018). The
difference between our approach and that of Shlezinger et al. (2018) is our
focus on vector quantization, and our consideration of comparator limitations
on the vector, rather than scalar, quantizer input.

Deep learning techniques have been used successfully in the domain of
quantization, for task-based quantization problems Shlezinger and Eldar (2019)
as well as particular use cases (MIMO channel estimation in Shohat et al.
(2019) or high-dimensional signal recovery from one-bit quantization in
Khobahi et al. (2019)). The idea of adapting the quantizer design to the
subsequent task, or task-based quantization, is explored in Shlezinger et al.
(2019) and studied in the context of hardware limitations in Shlezinger et al.
(2018) Shlezinger et al. (2019). Tree-Structured Vector Quantization uses a
set of scalar quantizers to gain speed, but comes at a disadvantage on the
shape of quantization regions.

Signal recovery from one-bit quantized signals is a topic previously
investigated in the literature, such as in Kamilov et al. (2012) and Plan and
Vershynin (2013). An approach to one-bit signal recovery using deep network is
considered in Khobahi et al. (2019). These approaches can be interpreted as
one-dimensional formulations of the CLVQ quantization problem we study here.
Algorithm for the design of quantizers under quadratic distortion is studied
by Max and Lloyd for scalar quantization Max (1960); Lloyd (1982), and
extended to vector quantization by Linde,Buzo, and Gray Linde (1980).

###  1.2 Contribution

Our main contributions, in the following paper, consist in the definition of
the Comparison-Limited Vector Quantization (CLVQ) problem in its full
generality and identify the connections between this problem and the classic
quantization problem. In particular, we highlight a connection between the
CLVQ problem and combinatorial geometry and provide novel insights on the
ultimate performance of low-resolution quantization.

We also provide a first class of algorithms for the CLVQ design: although not
optimal, this first approach investigates the combinatorial geometric aspects
of the CLVQ design.  Numerical evaluations are provided for the case of
uniform and Gaussian i.i.d. sources. The performance of the proposed quantizer
is compared with the classic optimal vector quantizer design obtained by
Linde-Buzo-Gray (LBG) algorithm Linde (1980).

###  1.3 Paper Organization

The paper is organized as follows: Sec. 2 presents related results in the
literature and introduces some useful combinatorial notions. Sec. 3 presents
the CLVQ model and the performance evaluation for this model. Sec. 4
introduces a class of algorithms for the CLVQ problem. Sec. 5 provides
relevant numerical evaluations for the case of an iid Gaussian source to be
reconstructed under quadratic distortion.  Finally, Sec. 6 concludes the
paper.

###  1.4 Notation

In the remainder of the paper, all logarithms are taken in base two. With
$\mathbf{x}=[x_{1},\ldots,x_{N}]\subseteq{\cal X}^{N}$ we indicate a sequence
of elements from ${\cal X}$ with length $N$. The notation $\mathbf{x}_{i}^{j}$
indicates the substring $[x_{i},\ldots,x_{j}]$ of $\mathbf{x}$. The function
${\rm sign}(\mathbf{x})$ returns a vector with values in $\\{-1,+1\\}$ which
equals the sign of each entry of the vector $\mathbf{x}$. The set
$\\{1,\ldots,N\\}$ is indicated as $[N]$.  Finally, the 2-norm operator is
denoted as $\|\cdot\|$, expectation as $\mathbb{E}[\cdot]$, and ${\bf
1}(\cdot)$ is used for the indicator function.

##  2 Related Results

In this section let us review three topics of relevance in the remainder of
the paper: (i) rate-distortion theory, that is the information theoretical
study of the quantization, (ii) quantizer design algorithms, such as the LBG
algorithms which we use as a point of comparison of our results, and (iii)
some combinatorial geometrical notions that are useful in describing quantizer
configuration in the CLVQ problem.

###  2.1 Rate-Distortion Theory

The theoretical framework we introduce for the study of the CLVQ problem–see
Sec. 3– expands on the classic formulation of rate-distortion theory (Cover,
1999, Ch. 10 ). Let us briefly review some relevant results in this section.

Broadly speaking, rate-distortion theory is the information theoretical study
of A2D conversion and addresses the quantization performance in the limit of
an infinite source sequence when the number of bit-per-sample is kept
constant.

Consider the system model in Fig. 2: here the source sample
$X^{n}=[X_{1},\ldots,X_{n}]$ with support ${\cal X}^{n}$ is represented
through the index $m_{n}$ using the source encoder mapping

$\displaystyle f_{\rm enc}:{\cal X}^{n}\rightarrow[2^{\lfloor nR\rfloor}],$
(1)

so that the cardinality of $m_{n}$ is $2^{nR}$ when the source sequence has
length $n$. The index $m_{n}$ is communicated to a source decoder through a
noiseless channel.  In turn, the source encoder produces the reconstruction
$\mathbf{\widehat{X}}_{n}$ through the mapping

$\displaystyle f_{\rm dec}:\ [2^{\lfloor nR\rfloor}]\rightarrow\widehat{{\cal
X}}^{n},$ (2)

so as to minimize the distortion measure $\rho_{n}(X^{n},\widehat{X}^{n})$
between the original source sequence $X^{n}$ and its reconstruction
$\widehat{X}^{n}$. In (2), $\rho_{n}$ for $n\in\mathbb{N}$ is a positive and
bounded function which captures the distortion, or degradation, between
$X^{n}$ and $\widehat{X}^{n}$. The best attainable performance when
reconstructing the original symbol for a given value of $R$ is evaluated
through the _rate-distortion function_ at rate $R$ and block-length $n$,
defined as

$\displaystyle D_{n}(R)=\inf\rho_{n}(X^{n},\widehat{X}^{n}),$ (3)

where the infimization is over all source encoder/decoder mappings as in
(1)/(2) . One is generally interested in the performance when $n$ goes to
infinity, since it provides an upper bound to the ultimate compression
performance. We refer to this limit as the _rate distortion function_ at rate
$R$:

$\displaystyle D(R)=\liminf_{n\rightarrow\infty}\ \ D_{n}(R).$ (4)

It is well known that, when (i) the source sample is i.i.d. from the
distribution $P_{X}$ and (ii) the distortion function is an additive
distortion function, i.e.

$\displaystyle\rho_{n}(X^{n};\widehat{X}^{n})=\sum_{i\in[n]}\rho(X_{i},\widehat{X}_{i}),$
(5)

for some positive and bounded $\rho(x_{i},\widehat{x}_{i})$, we have

$\displaystyle D(R)=\inf_{P_{\widehat{X}|X},\ R{\geq}I(X,\widehat{X})}\
\mathbb{E}\left[\rho(X,\widehat{X})\right],$ (6)

where $I(X,\widehat{X})$ indicates the _mutual information_ between the source
$X$ and the reconstruction $\widehat{X}$. As an example of the result in (6)
consider the case in which $X$ has a standard normal distribution to be
reconstructed under quadratic distortion, then

$\displaystyle D(R)=2^{-2R},$ (7)

so that $D\in[0,1]$.

source encoder $\mathbf{X}_{n}$ source decoder
$m_{n}\in[2^{dR}]$$\mathbf{\widehat{X}}_{n}$

Figure 2: The rate-distortion problem in Sec. 2.1.

###  2.2 Quantization Design Algorithms

The approach in Sec. 2.1 provides the theoretical limit to the quantization
performance. In order to approach such performance, various algorithms have
been proposed in the literature for the design of high-dimensional quantizers.

In the following, we shall compare the performance of the CLVQ design with the
standard LBG algorithm which we shall describe next. The LBG algorithm, also
known as the Max-Lloyd algorithm or $k$-means clustering algorithm, is an
iterative algorithm for the design of quantizers for the case of quadratic
distortion Max (1960); Lloyd (1982); Linde (1980). It is known to approach the
optimal solution although this problem is NP-hard in the Euclidean space. The
pseudo-code for the algorithm is provided in Algo. 1: the algorithm performs
an alternate optimization between the reconstruction regions and the position
of the reconstruction points. In the first step, the sample space ${\cal X}$
is partitioned by a Voronoi tessellation with the reconstruction points as
generators. This means that each point in the sample space is assigned to the
closest reconstruction points. reconstruction point has a region, and every
point in space belongs to the region of the nearest reconstruction point. This
procedure results in a set of regions $\\{{\cal V}_{m}\\}_{m\in[M]}$ so that
each point in the region ${\cal V}_{m}$ contains the points that have the
minimum distortion when reconstructed as the point $\widehat{c}_{m}$, instead
of any other reconstruction point. In the second step, the position of the
reconstruction points $c_{m}$ is optimized: each reconstruction point $c_{m}$
is chosen as the mean of the sample points in the region ${\cal V}_{m}$. The
two steps are repeated in succession until convergence is achieved. Various
procedures can be designed to initialize the initial positions of the set of
reconstruction points $\\{\widehat{c}_{m}\\}_{m\in[M]}$ as it is known that
the choice of initialization set its crucial to attain fast convergence.

In the literature, many variations of the LBG algorithm as been studied with
applications ranging from unsupervised learning to low-ranking approximations.
Let us mention two variations of interest as a generalization of the settings
in this paper.

Entropy-constrained vector quantizers (ECVQs) are quantizers in which the
quantization objective is the distortion minimization subject an additional a
bound on the entropy of the quantizer output. ECVQ are useful when they are
used in tandem with variable-rate noiseless coding systems to provide locally
optimal variable-rate block source coding. Another class of quantizers of
interest are lattice quantizers. In lattice quantizers are quantizers in which
the set of reconstructions points form a lattice. Note that the problem of
determining the closest reconstruction point to a given quantizer input is, in
general, a nearest neighbor problem. This class of problems has a complexity
linear in the number of neighbors: when a quantizer has high dimension, this
complexity is not amenable to real-time applications. For this reason, one is
often interested in introducing further structure in the set
$\\{\widehat{c}_{m}\\}_{m\in[M]}$ with the aim of speeding up the quantization
operation. One class of such quantizers are lattice vector quantizers, that
is, quantizers in which the set of reconstruction points form a lattice. It
has been shown that lattice quantizers have good performance while also
allowing for a fast quantization Conway and Sloane (1982); Zamir and Feder
(1992).

Algorithm 1 Max-Lloyd/$k$-means clustering/Linde-Buzo-Gray (LBG) Algorithm

1: Input: source distribution $P_{X^{n}}$, number of reconstruction points
$M=2^{\lfloor nR\rfloor}$, maximum number of iterations $T_{\max}$

2: Output: position of the reconstruction points
$\\{\widehat{c}_{m}\\}_{m\in[M]}$, Voronoi region of each reconstruction point
$\\{{\cal V}_{m}\\}_{m\in[M]}$

3: randomly pick the initial reconstruction points
$\\{\widehat{c}_{m}(0)\\}_{m\in[M]}$, i.e.
$\displaystyle\widehat{c}_{m}(0)\sim P_{X^{n}}$

4: for $t$ in $[T_{\max}]$ do

5: update the Voronoi regions $\\{{\cal V}_{m}(t)\\}_{m\in[M]}$ as
$\displaystyle{\cal V}_{m}(t)=\left\\{x\in{\cal X}{\ \rm{s.t.}\
}(c_{m}(t-1)-x)^{2}\leq(c_{j}(t-1)-x)^{2},\ \forall j\in[M]\setminus
m\right\\}$ for all $m\in[M]$

6: update the position of the reconstruction points $\\{c_{m}(t)\\}_{m\in[M]}$
as  $\displaystyle c_{m}(t)=\frac{\mathbb{E}[X|X\in{\cal
V}_{m}(t)]}{\mathbb{E}[1_{\\{X\in{\cal V}_{m}(t)}\\}]}$

7: end for

###  2.3 Some combinatorial geometric notions

This section briefly introduces a few combinatorial concepts useful in the
remainder of the paper to describing the geometrical properties of the
proposed quantization scenario.  In this section and through the paper we
mostly use the nomenclature and notation of Stanley (2006) A _hyperplane
arrangement_ ${\cal A}$ is a finite set of $n$ affine hyperplanes in
$\mathbb{R}^{m}$ for some $n,m\in\mathbb{N}$. where $\mathbf{A}$ is obtained
by letting each row $i$ equals $\mathbf{a}_{i}^{T}$ and letting
$\mathbf{b}=[b_{1}\ldots b_{n}]^{T}$. A _region_ , ${\cal R}_{i}$ of an
arrangement ${\cal A}$ is a connected component of the complement
$\overline{{\cal A}}$ of the hyperplanes, defined as

$\displaystyle\overline{{\cal A}}=\mathbb{R}^{m}-{\cal A}.$ (8)

Let $\mathbf{r}({\cal A})$ be the number of regions in which the hyperplane
arrangement ${\cal A}$ divides the space $\mathbb{R}^{m}$, so that

$\displaystyle\overline{{\cal A}}=\bigcup_{i=1}^{\mathbf{r}({\cal A})}{\cal
R}_{i}.$ (9)

A plane arrangement ${\cal A}$ is said to be in General Position (GP) if and
only if:

$\\{H_{1},...,H_{p}\\}\subseteq{\cal A},p\leq n\Rightarrow
dim(H_{1}\cap...\cap H_{n})=n-p$ $\\{H_{1},...,H_{p}\\}\subseteq{\cal
A},p>n\Rightarrow H_{1}\cap...\cap H_{n}=\emptyset$

(Dimca and PARUSINSKI, 2017, Th. 1.2). A hyperplane arrangement of size $n$ in
$\mathbb{R}^{m}$ divides $\mathbb{R}^{m}$ into at most $\mathrm{r}(m,n)$
region for

$\mathrm{r}(m,n)=\sum_{i=0}^{m}{n\choose i}\leq 2^{n}.$ (10)

(Stanley, 2006, Prop. 2.4). A hyperplane arrangement of size $n$ in
$\mathbb{R}^{m}$ where all the hyperplanes pass through the origin divides
$\mathbb{R}^{m}$ into at most $\mathrm{r_{{\bf 0}}}$ regions for

$\displaystyle\mathrm{r_{{\bf 0}}}(n,m)=2\sum_{i=0}^{m-1}{{n-1}\choose{i}}.$
(11)

(Stanley, 2006, Prop. 2.4). Let ${\cal A}$ be a hyperplane arrangement of size
$l$ in $\mathbb{R}^{m}$ and consider a hyperplane arrangement ${\cal B}$ of
size $dl$ with $d\in\mathbb{N}$ hyperplanes parallel to each of the
hyperplanes in ${\cal A}$, then ${\cal B}$ divides $\mathbb{R}^{m}$ into at
most $\mathrm{r_{p}}$ regions for

$\mathrm{r_{p}}(m,n,d)=\sum_{i=0}^{m}\binom{l}{i}d^{i}\leq(1+d)^{l}.$ (12)

A necessary condition to attain the equality in Lem. 2.3, 10, and Lem. 2.3 is
for the hyperplane arrangement ${\cal A}$ to be in GP. As we shall see in the
next section,

##  3 Comparison-Limited Vector Quantization (CLVQ) Problem

After the introductory notion in Sec. 2 let us introduce the Comparison-
Limited Vector Quantization (CLVQ) problem, also conceptually depicted in Fig.
3.

$d$-parser$X_{i}$$\mathbf{v}_{1}$$\mathbf{v}_{2}$$\mathbf{v}_{3}$$\mathbf{v}_{4}$$\mathbf{X}_{n}=X_{di+1}^{d(i+1)}$$+$$+$$+$$+$$t_{4}$$t_{3}$$t_{2}$$t_{1}$${\rm
sign}(\cdot)$${\rm sign}(\cdot)$${\rm sign}(\cdot)$${\rm
sign}(\cdot)$$n$-parser$Y_{1n}$$Y_{2n}$$Y_{3n}$$Y_{4n}$ source encoder
$\mathbf{Y}_{n}=[Y_{1n},\ldots,Y_{kn}]$ source decoder
$m_{n}\in[2^{dR}]$$\mathbf{\widehat{X}}_{n}$

Figure 3: The Comparison-Limited Vector Quantization (CLVQ) problem with $k=4$
quantizers.

A source sequence $\\{X_{i}\\}_{i\in\mathbb{N}}$, where each $X_{i}$ has
support ${\cal X}$, is parsed in super-symbols
$\\{\mathbf{X}_{n}\\}_{n\in\mathbb{N}}$ of dimension $d$ with
$\mathbf{X}_{n}=[X_{d(n-1)+1},\ldots,X_{dn}]$ for $n\in\mathbb{N}$ where $d$
is referred to as the _dimension_ of the vector quantizer. The $j^{\rm th}$
comparator computes a linear combination of each super-symbol $\mathbf{X}_{n}$
and outputs a signal $Y_{jn}$ as

$\displaystyle Y_{jn}={\rm
sign}\left(\mathbf{v}_{j}\mathbf{X}_{n}+t_{j}\right),$ (13)

for $j\in[k]$. The value $k$ is called _resolution_ of the vector quantizer;
$\mathbf{v}_{j}\in\mathbb{R}^{d}$, $t_{j}\in\mathbb{R}$ are fixed and known.
For more convenience, we express the $k$ outputs of the quantizers of (13) as

$\displaystyle\mathbf{Y}_{n}={\rm
sign}\left(\mathbf{V}\mathbf{X}_{n}+\mathbf{t}\right),$ (14)

where $\mathbf{V}=\in\mathbb{R}^{k\times d}$ is such that its $i^{\rm th}$ row
is equal to the vector $\mathbf{v}_{i}$ from (13). Similarly,
$\mathbf{t}\in\mathbb{R}^{k}$ has the $i^{\rm th}$ entry equal to $t_{i}$ and
$\mathbf{Y}_{n}=[Y_{1n},\ldots,Y_{kn}]$. The set $[\mathbf{V},\mathbf{t}]$ is
referred to as the configuration of the quantizer. The super-symbol
$\mathbf{Y}_{n}$ is provided to a source encoder that produces a bit-
restricted representation of the quantizers’ output as $m_{n}\in[2^{\lfloor
dR\rfloor}]$ where $R$ is the name given to the _rate_ of the quantizer
through the _source encoding mapping_

$\displaystyle f_{\rm enc}:\ \\{-1,+1\\}^{k}\rightarrow[2^{\lfloor
dR\rfloor}].$ (15)

The message $m_{n}$ is sent into the source decoder which outputs a
reconstruction of the source super-symbol $\mathbf{X}_{n}$ as
$\mathbf{\widehat{X}}_{n}=[\widehat{X}_{dn+1},\ldots,\widehat{X}_{d(n+1)}]$
with $\widehat{X}_{i}\in\widehat{{\cal X}}$ , through the _source decoding
mapping_

$\displaystyle f_{\rm dec}:\ [2^{\lfloor dR\rfloor}]\rightarrow\widehat{{\cal
X}}^{d}.$ (16)

We measure the effectiveness of a vector quantizer with a measure of
distortion

$\displaystyle\rho_{n}(X^{n};\widehat{X}^{n}):\ {\cal
X}^{n}\times\widehat{{\cal X}}^{n}\rightarrow\mathbb{R}^{+},$ (17)

for $n\in\mathbb{N}$ which is assumed positive and non-decreasing in $n$. For
a given configuration of the linear combiners $[\mathbf{V},\mathbf{t}]$,
source encoder/decoder mappings $f_{\rm enc}$/$f_{\rm dec}$, and given a
distortion measure between input and reconstruction sequence
$\rho_{n}(X^{n};\widehat{X}^{n})$, the performance of the quantizer is
evaluated as

$\displaystyle\overline{\rho}=\limsup_{n\rightarrow\infty}\
\rho_{n}(X^{n};\widehat{X}^{n}).$ (18)

The optimal quantizer performance regarding the distortion $\overline{\rho}$,
dimension $d$, resolution $k$ and rate $R$ can be written as

$\displaystyle D(d,k,R)=\inf\ \ \overline{\rho},$ (19)

where the infimization is over all configuration of the linear combiners [V,t]
and all the source encoder/decoder mappings.

Similarly to (4), one is often interested in determining the limiting
performance as the dimension of the quantizer grows to infinity. In this case,
it is natural to let the number of quantizer grow exponentially with the
dimension, so that the number of comparators per dimension is kept constant,
with ration approximatively equal to $\alpha$.  Let $k=\lceil\alpha d\rceil$,
then the comparison-limited distortion-rate function is defined as

$\displaystyle D_{\alpha}(R)=\lim_{d\rightarrow\infty,\ k=\lceil\alpha
d\rceil}D(d,k,R),$ (20)

for $D(d,k,R)$ in (19), that is $D_{\alpha}(R)$ is the minimum distortion
attainable as the quantizer dimension grows to infinity while
$k/d\approx\alpha$.

###  3.1 Underlying Assumptions

The underlying assumptions in our problem formulation are as follows: op-amp
voltage comparators are employed in nearly all analog-to-digital converters to
obtain multilevel quantization. A reconfigurable receiver front-end might be
able to re-configure the comparators’ inputs so as to perform more complex
operations. Given the receiver’s ability to partially reconfigure its
circuitry depending on the channel realization, we wish to determine which
configuration of the comparators yields the largest capacity. We restrict our
results to the case of linear analog pre-coding, as the capacity of the model
without output quantization constraints can be attained through linear pre-
coding strategies. The extension to non-linear processing is possible but not
pursued here.

###  3.2 Discussion

Let us connect the CLVQ problem as defined above with the results in Sec. 2.
As it can be gathered by comparing Fig. 2 versus Fig. 3, the CLVQ problem
explicitly embeds a model of the hardware architecture performing A2D.  In
particular, in this formulation, we assume that quantization is performed
through a bank of op-amp amplifier. Each op-amp amplifier receives a linear
combination of the quantizer input and a bias.

It is immediate to connect the problem formulation in Fig. 3 with the setting
in Sec. 2.3. The set of the comparators output $\mathbf{Y}_{n}$ in (14)
represent the membership function of each quantizer input $\mathbf{X}_{n}$
with respect to each of the hyperplanes $[\mathbf{v}_{i},t_{i}]$. This output
is then compressed as in the classic rate-distortion problem in Sec. 2.1. When
the number of sign quantizers is sufficiently large, a hyperplane can be used
to separate any pair of reconstruction points, as in the Voronoi region of
Algo. 1. In this regime, the rate in (6) is determined by the number of
regions induced by the hyperplane arrangement. Given the considerations above,
it is now useful to imagine the CLVQ problem in the context of Fig. 1: each
comparator induces a hyperplane partitioning the sample space. The set of
comparators then can be imagined as a hyperplane arrangement partitioning the
sample space. The set of comparator outputs is then compressed as in the
classic rate-distortion problem, thus further reducing the number of bits
required to represent $\mathbf{Y}^{n}$ by exploiting its distribution.

Let us next provide some more specific remarks.

Extension to the vector version In the following we consider the case in which
the source sequence $\\{X_{i}\\}_{i\in\mathbb{N}}$ is scalar. The case of a
vector input sequence is not considered here. For the vector case, the
dependency among different dimensions can be addressed in different manner:
two come readily to mind. One approach would be to choose the linear
combination as the whitening transformation that make the covariance of the
input symbol unitary. Alternatively, the different dimension can be processed
in parallel by the comparator network while the source encoder is tasked with
accounting for the correlation across dimensions.

Entropy coding: Note that the CLVQ problem formulation also captures the
digital processing following A2D conversion. After quantization, the quantized
samples are further compressed using a digital source compressor, such as an
entropy coder. This is captured by the source encoder/decoder in (15)/(16) and
as in Fig. 3. This is similar to the ECVQ discussed in Sec. 2.2.

Limitations in the comparator configurations: In (19) we assume that the CLVQ
performance can be optimized over all possible configuration of the linear
combiners $[\mathbf{V},\mathbf{t}]$. As a configuration of a linear combiner
corresponds to an analog circutry implementation, this assumption might be
unfeasible in some scenarios, as generating precise linear combinations and
voltage references might require a high circuit complexity. To accommodate for
such further limitations, on can further restrict the optimization in (19)
over a set of configurations which can be efficiently implemented. Such an
approach is similar to that in Rini et al. (2017) and further discussion on
this approach can be found there.

Lattice quantization for CLVQ: Note that one can extend the construction of
lattice codes for quantization to the CLVQ problem by considering regular
hyperplane arraignments. In this case, it is possible to similarly define
computationally-efficient quantization algorithms as discussed in Sec. 2.2.

##  4 A CLVQ Quantizer Design Algorithm

In the following, we shall introduce an algorithm for the design of the CLVQ
linear combiner configuration $[\mathbf{V},\mathbf{t}]$ for given value of
$d$, $k$, and $R$ in (19). In other words, the algorithm attempts to
numerically determine the optimal solution to (19) by assuming that the
optimal choice of source encoder/decoder function is performed in a successive
optimization step.

For the design of this numerical optimization algorithm, we take inspiration
from the quantizer design algorithms discussed in Sec. 2.2, although the CLVQ
problem is more general than considered is Sec. 2.2. Indeed, we shall leverage
the notion in Sec. 2.3 to more efficiently represent and manipulate the linear
combiner configuration during the optimization. The theoretical analysis of
the asymptotically optimal solution as in (20) is left for future research.

For simplicity, in the design of our algorithm, we simplify the general
setting of Sec. 3 as follows:

  1. (a)

i.i.d. sources – with distribution unimodal distribution $P_{X}$,

  2. (b)

quadratic distortion – also known as Mean Squared Error (MSE) distortion, i.e.

$\displaystyle\rho_{n}(X^{n};\widehat{X}^{n})=\frac{1}{n}\sum_{i=1}^{n}\|X_{i}-\widehat{X}_{i}\|^{2}.$
(21)

  3. (c)

infinite compression rate – that is $R=\infty$ in (19), i.e. $D(d,k,\infty)$.

Regarding assumption (b), this assumption is relaxed in Sec. 5.3.3, here we
discuss the use of entropy as a measure of distortion measure in the numerical
optimization. Assumption (c) is chosen so as to simply the design of the
comparator configuration. An effective, but sub-optimal, approach to address
this constraint would be choose the encoding and decoding function which
implement an efficient lossy-compression algorithm for discrete data – see
Nelson and Gailly (1995).

Note that a more efficient approach to address the case $R<\infty$ would be to
consider an entropy constrain, as discussed in Rem. 3.2.  The design of such
an algorithm would be similar to the problem of entropy-constrained
quantization as in Chou et al. (1989). This further design step is not
considered at the present.

Given the three assumptions above, (19) simplifies to

$\displaystyle
D(d,k,\infty)=\sum_{i=1}^{d}\mathbb{E}\left[\|X_{i}-\widehat{X}_{i}\|^{2}\right],$
(22)

where $\mathbf{X}=[X_{i},\ldots,X_{d}]$ and
$\mathbf{\widehat{X}}=[\widehat{X}_{1},\ldots,\widehat{X}_{d}]$ for
$X_{i},\widehat{X}_{i}$ in (22) are i.i.d. distributed, so that the subscript
$n$ can be dropped.111That is, every super-symbol $\mathbf{X}_{i}$ is
quantized in the same manner, regardless of $n$. For the problem in (22), the
linear combiner optimization problem is mapped to a combinatorial geometric
optimization problem as in Sec. 2.3 in which one is interested in choosing the
hyperplane arrangement which divides the space in $\\{{\cal R}_{i}\\}$ regions
for which MSE distortion  is minimized as in Algo. 1 for the case in which
each source vector $\mathbf{X}_{n}$ is mapped to the closest centroid in the
set $\\{c_{i}\\}_{i\in[\mathbf{r}({\cal A})]}$.  The set of centroids is
obtained from the hyperplane arrangement ${\cal A}$ as

$\displaystyle c_{i}=\frac{\mathbb{E}[X|X\in{\cal
R}_{i}]}{\mathbb{E}[1_{\\{X\in{\cal R}_{i}}\\}]},$ (23)

similarly to the LBG algorithm.

Let us further detail the proposed approach: our algorithm is divided in two
steps: (i) an initialization step –discussed in Sec. 4.1 – and (ii) an
optimization step – detailed in Sec. 4.2.

Generally speaking, the optimization step is computationally complex, so the
role of the initialization step is to choose of a set of initial
configurations of quantizer that can span the large set of choices for the
combiner configuration. By considering multiple restarts of the optimization
point with a different initialization, one hopes to converge to an optimal
solution and avoid local minima.

Note that there is a large number of symmetries in the linear combiner
configurations, since one can permute each quantizer configuration and input
values and obtain an equivalent configuration. This induces a number of local
minima which hinder the convergence of the optimization step to the globally
optimal solution. For this reason, the choice of the set of combiner
configuration used for initialization is crucial to span the space of possible
solutions.

###  4.1 Initialization methods

In the following, we explore two possibilities for initialization:

– random initialization: in which the set of initial linear combiner
configurations selected using random samples from the source distribution.

– genetic initialization: in which we utilize a genetic algorithm to generate
the set of initial configurations.

####  4.1.1 Random initialization

This first approach to the generation of initialization steps is the simplest,
but yet provides good results. To obtain random hyperplanes that cut through
dense zones of the source distribution, it is sufficient to randomly generate
$d$ points following that same source distribution and the hyperplane passing
through these points can be used as a linear combiner configuration. By
repeating this operation $k$ times, one obtains an initial configuration for
every comparator. This random initialization has constant time complexity and
allows us to quickly generate a large number of possible candidate
arrangements.

####  4.1.2 Genetic Initialization

Given the poor performance of the random initialization in Sec. 4.1.1 and the
computational complexity of the exhaustive listing of all possible
configuration in App. C, we next consider a genetic algorithm that retains
elements of the two previous approaches. A genetic algorithm is a meta-
heuristic taking inspiration from genetic evolution to apply some processes
such as selection, crossover and mutation to optimization problems. It
consists of a few steps that are repeated over several iterations. An initial
population of potential solutions has to be generated. At each iteration, a
part of the current generation is selected to breed the next one, based on a
fitness function measuring its efficiency regarding the problem. Then, another
generation is formed by performing crossover (which corresponds to mixing some
elements between pairs of solutions) and mutation (small alterations of
solution) on the current generation.

In our case, after generating a few starting configurations according to the
random approach in Sec. 4.1.1, we proceed to genetic selection in three
distinct stages:

  1. (a)

select a set of configurations based on their fitness (the MSE in our case)

  2. (b)

perform crossover on these configurations to obtain the next generation with
improved fitness

  3. (c)

apply random mutations to the resulting configurations

$\bullet\ $ In step (a), the configurations are reordered by increasing
distortion (we use the opposite of distortion as fitness function).The $k$
first configurations are kept for the next step, $k$ being a fixed parameter.
In our simulations, we set $k$ to $80\%$ of the number of configurations.

$\bullet\ $ In step (b), we cross pairs of arrangements following one of two
policies. The first one consists in randomly selecting hyperplanes to form
pairs. With the second policy, we use a dissimilarity value to form the pairs
of hyperplanes. This dissimilarity value is given by $\theta\times
d(p_{1},p_{2})$, where $\theta$ is the angle between the hyperplanes; $d$ is
the euclidean distance; and $p_{1},p_{2}$ are the nearest point of each
hyperplane to the centroid of the source distribution. Measures of the
dissimilarity of each pair of hyperplanes are placed in a matrix, and
hyperplanes are associated two by two in a way that produces the lowest
average dissimilarity. Once pairs are formed, a new configuration is formed by
selecting one hyperplane from each pair.

$\bullet\ $ Step (c) consists in mutating the resulting configurations.To
apply mutation, we first form vectors of random numbers following a normal
probability law with mean $\mu=1$ and with a standard deviation $\sigma=0.2$
that was empirically chosen. Mutation is only applied to the half of
configurations that have the lowest fitness score. The intuition for applying
mutation this way is that solutions with a good fitness value might be made
worse by a mutation, while the ones which fitness is already not good have
better chances of being improved by the mutation.

These steps are repeated for a pre-determined number of iterations. Each step
has a constant complexity, so the complexity of each iteration is linear with
the number of initial hyperplane arrangements. This process provides us a set
of arrangements with low distortion

By repeating the process multiple times, we can accordingly generate a set of
good initialization points.

Pseudo-code of the proposed genetic algorithm is also described in Algo. 2.

Algorithm 2 Genetic algorithm

1: Input: ${T_{\max}}$, $M$ and $k<M\in\mathbb{N}$, a set $\\{{\cal
A}_{0,m}\\}_{m\in[M]}$ of random hyperplane arrangements, a fixed parameter

2: Output: one hyperplane arrangement with low distortion

3: for $t$ in $[T_{\max}]$ do

4: order $\\{{\cal A}_{t,m}\\}_{m\in[M]}$ by decreasing fitness

5: $m\leftarrow k$

6: form pairs $({\cal A}_{t,m_{1}},{\cal A}_{t,m_{2}})_{i}$ (read below about
the 2 existing policies to chose $(m_{1},m_{2})$ pairs)

7: for each pair $(m_{1},m_{2})$ do

8: perform crossover between ${\cal A}_{t,m_{1}}$ and ${\cal A}_{t,m_{2}}$

9: end for

10: order $\\{{\cal A}_{t,m}\\}_{m\in[M]}$ by decreasing fitness

11: for $m\in[M]$ do

12: if $m<\frac{M}{2}$ then

13: ${\cal A}_{t+1,m}\leftarrow{\cal A}_{t,m}$

14: else

15: ${\cal A}_{t+1,m}\leftarrow mutate({\cal A}_{t,m})$

16: end if

17: end for

18: end for

19: select ${\max(\\{{\cal A}_{T_{\max},m}\\}}_{m\in[M]})$

An idea that naturally arises would be to consider all possible hyperplane
configurations and simply optimize each of them. In App. A, we explain some
notions of matroid theory that provide a framework for the exhaustive
generation of hyperplane configurations. In App. B, we define a way to
describe the structure of hyperplane arrangements with graphs. Then in App. C,
we expose how one can progressively add nodes to this graph representation in
order to generate graph configurations, and the reasons for this approach not
being applicable in practice.

###  4.2 Optimization step

Being inspired from LBG algorithm (described in 1), our quantizer design
algorithm also features two steps that are being iterated until convergence:
(i) the optimization of the set of reconstruction points
$\mathbf{\widehat{X}}$ for a given combiner configuration $[\mathbf{V}\
\mathbf{t}]$ and (ii) the optimization of the combiner configuration for a
given set of reconstruction points. Step (i) consists in choosing the
reconstruction points as the centroids (with regards to the source
distribution) of the regions defined by the hyperplane arrangement
$[\mathbf{V}\ \mathbf{t}]$ in the space $\mathbb{R}^{d}$ Orlik and Terao
(2013). The optimization step (ii) presents more complexity. In that step we
update the hyperplane arrangement $[\mathbf{V}\ \mathbf{t}]$ using two
approaches: a _global_ configuration update and a _local_ one as described in
Algorithm 3. The global configuration update consists in adding a random
perturbation to every hyperplane, perturbation which variance decreases over
iterations. With the local configuration update, one hyperplane is randomly
selected among the hyperplanes of the configuration. Its position and
orientation are then set so as to minimize the distortion of the resulting
output. In order to do that, we calculate the resulting MSE for several values
of the hyperplane coefficients and apply interpolation between these values to
find the best ones. At each iteration, one of the two configuration update
approaches is chosen at random. The probability of choosing the local
configuration update augments with the number of iterations, at a speed
depending of an empirically fixed parameter $s$.

Figure 4:  Hyperplane arrangement at different steps of the optimization
algorithm. The source distribution is unitary uniform and its support is
represented by the grey square. The crosses show the centroids of each region.
The regions are becoming more and more regular over iterations.

The reasoning behind differentiating local and global update is the following.
Consider a lower triangular matrix $M$ of size
$\mathrm{r}(d,k)\times\mathrm{r}(d,k)$ and let the element in position
$i\times j$ in $M$ be equal to one of the hyperplanes separating
$\widehat{X}_{i}$ and $\widehat{X}_{j}$ if such reconstruction point exists or
equal to zero otherwise. 222Note that some hyperplane arrangements induce less
that $\mathrm{r}(d,k)$ regions: we assume that there exists a natural
numbering of the possible $\mathrm{r}(d,k)$ regions. The matrix $M$ can be
thought of as one among a finite number of ways in which hyperplanes separate
the reconstruction points.

In this view, the local update maximizes the quantizer performance in a given
value of $M$. The global update, instead, allows “hyperplanes to jump over
centroids”, resulting in a different matrix $M$. This is illustrated in Fig.
4, where passing from step a. to step b. requires a global update (that is
changing the structure of the arrangement) but passing from step b. to c. and
c. to d. only require local updates (optimizing the configuration).

Algorithm 3 Quantizer design algorithm

Input: $s$ and $[T_{\max}]\in\mathbb{N}$, a hyperplane arrangement ${\cal
A}_{0}$, $d$ the dimension

Output: a hyperplane arrangement for quantization

for $t$ in $[T_{\max}]$ do

if $random(0,1)<\exp(-s.t)$ then

let ${\cal H}_{i,t}$ be the $i^{th}$ hyperplane of ${\cal A}_{t}$

let ${\cal V}_{i,t}$ be a vector of $d$ realisations of ${\cal
N}(0,\sigma_{t})$ (with $\sigma_{t}$ a parameter decreasing over $t$)

for all $i$, ${\cal H}_{i,t}\leftarrow{\cal H}_{i,t}+{\cal V}_{i,t}$

else

select one hyperplane ${\cal H}$ of ${\cal A}_{t}$

for $i$ in $[d]$ do

let $p_{i}$ be the $i^{th}$ coefficient of ${\cal H}$ and compute distortion
for alternate values of $p_{i}$

interpolate the function $f:p_{i}\mapsto distortion$

$p_{i}\leftarrow\min(f)_{p_{i}\in\mathbb{R}}$

end for

end if

end for

###  4.3 Discussion

A fundamental element in step (ii) of the design algorithm is the evaluation
of the distortion (the MSE in our case) for a given hyperplane arrangement and
set of reconstruction points. Given numerical precision limitations, our MSE
evaluation has to be an approximation method, that uses numerical integration
methods and particle filters.  More precisely, we generate random points
following the source distribution, until a certain number of points is reached
for every region of the hyperplane arrangement. For each random point
generated, the output of the quantizer is the reconstruction point of the
region in which that random point is located. Finally, the average of the
squared distance between the random points and their corresponding output
point is taken as the estimation of the MSE.

We use a similar approach to estimate the centroids of the regions formed by a
hyperplane arrangement, since that estimation requires numerical integration
too. To estimate the centroid positions, random points are also generated
following the source distribution. Each point is assigned to the region it is
in. We can then obtain the position of each centroid by computing the average
coordinates of all the points generated in each region.

Note that these two steps (estimation of the centroids and estimation of the
MSE) are done separately. In particular, they use distinct samples of the
input distribution.

With the alternance of steps (i) and (ii) of the algorithm, the hyperplane
configuration converges to an optimum. Multiple random restarts of the
algorithm occasionally lead to convergence to distinct local minima. These
minimal values are due either to a limitation in the precision of the
numerical integration or to a local minimum in the quantizer configuration. In
the next section, we shall further comment on this and other aspects of the
proposed optimization algorithm in view of numerical evaluation results.

##  5 Numerical evaluations

In this section, we numerically investigate the a few aspects of the proposed
algorithm. In Sec. 5.1 we will investigate the performance of the genetic
algorithm discussed in Sec. 4.1.2. In Sec. 5.2, we discuss the optimization
step in Sec. 4.2. Finally, in Sec. 5.3, we provide a performance comparison of
the proposed algorithm and the classic LBG algorithm in Sec. 2.2.

###  5.1 Genetic Initialization Performance

Let us return to the genetic initialization discussed in Sec. 4.1.2. The
nature of this genetic algorithm is such that the performance depends on the
number of comparator configurations in input to the algorithm (see line 1 of
Algo 2). The more are the available configurations, the bigger is the pool of
genes that can be used. On the other hand, adding too many configurations can
be computationally expensive.

Fig. 5 shows the result obtained by the genetic algorithm with an input of
dimension 3 and configurations of 3 or 4 hyperplanes in Fig. 5(a) and Fig.
5(b) respectively. Every run of the algorithm used a pool of 10
configurations. These configurations are then mixed following the steps
described in Sec. 4.1.2. The algorithm was run twice for 5(a) and twice again
for 5(b). Both (a) and (b) display these two runs of the algorithm in light
blue, and their average in dark blue.

$0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$22$$24$$26$$28$$30$$0.36$$0.38$$0.4$$0.42$$0.44$$0.46$$0.48$$0.5$iterationsdistortionaverage
of several executionsdistinct executions of the algorithm (a) Two random
restarts of the genetic algorithm in gray and their average in blue, on 30
iterations, 3 hyperplanes in 3 dimensions.

$0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$22$$24$$26$$28$$30$$0.3$$0.31$$0.32$$0.33$$0.34$$0.35$$0.36$$0.37$$0.38$$0.39$$0.4$iterationsdistortionaverage
of several executionsdistinct executions of the algorithm (b) Two random
restarts of the genetic algorithm in gray and their average in blue on 30
iterations, 4 hyperplanes in 3 dimensions.

Figure 5: Genetic algorithm evolution.

Note that in these two particular executions, in both Fig. 5(a) and Fig. 5(b),
that the distortion decreases very irregularly. This curve is composed of flat
sections and a few sudden drops: this suggests that improvements happen from
time to time when the right combination of two configurations is formed. For
this reason, added to the general propensity of genetic algorithms to converge
to local minima, this algorithm cannot be used for the whole optimization
process. However, its computational efficiency makes it a good pre-
optimization tool. Indeed, the distortion function is the often the most
computationally heavy part of the optimization. This implementation calls to
the distortion function $2\times n$ times per iteration (with n the number of
configurations), while the previous algorithm we described in Sec. 4 calls it
$40\times H\times d$ with $H$ the number of hyperplanes and $d$ the dimension.
Thus, this genetic algorithm can be used to quickly obtain a rough estimate of
the global minimum, before a more precise and more computation-costly
algorithm can finish the optimization.

###  5.2 Optimization Step Performance

Next, we numerically investigate the optimization step discussed in Sec. 4.2.
As discussed in this section, the proposed alternate optimization encounters a
number of local minima: Fig. 6 (showing the quantizers obtained by 2 distinct
random restarts of the optimization algorithm) demonstrates the possibility of
convergence to distinct local minima.While the arrangement in Fig. 6.a has 6
centroids, the one in Fig 6.b has 7, though both of them score similar
performance for a standard Gaussian source distribution. Regarding the
structure of the arrangement, however, these two configurations are rather
different, and the proposed algorithm is unable to go from the configuration
in Fig.6.a to the slightly better configuration 6.b.

Figure 6: 2 local minima obtained with the algorithm in Sec. 4. The center of
the Gaussian source distribution is marked by the grey dot, and the centroid
of each region is indicated by a cross.

###  5.3 LBG performance comparison

Finally, we come to perhaps the most interesting numerical evaluations: the
comparison between the proposed quantization approach and the classic LBG
algorithm. Such comparison indeed addresses the observation in Fig. 1 which is
the initial motivation behind the proposed approach.

This section presents the performance of the quantizer described in Sec. 3
with a configuration resulting from the algorithm in Sec. 4 for the case of
(i) standard Gaussian and (ii) unitary uniform distribution. We obtained
similar performance results in both cases. In all cases, we compared the
performance attained with the performance of the quantizer obtained by the
classic LBG algorithm.

$1$$1.2$$1.4$$1.6$$1.8$$2$$2.2$$2.4$$2.6$$2.8$$3$$3.2$$3.4$$3.6$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$Number
of hyperplanes or regions, divided by dimensionDistortionLBG algorithmproposed
algorithm Figure 7: Diagram showing the average distortion (Mean Squared-Error
here) values attained by the classic LBG algorithm (green) and the proposed
algorithm (blue) on a Gaussian input signal. The number $d$ of dimensions
ranges from 2 to 6, and the number of regions (for LBG) or hyperplanes (for
the proposed algorithm) ranges from $d$ to 8. The results are averaged over
several random restarts of the algorithms. For better readability, these
results are spread on a X-axis representing the ratio between number of
hyperplanes or regions, and dimension.

The code for this algorithm is available online: see Chataignon .

####  5.3.1 Gaussian distribution

$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$Number
of comparatorsDistortionLBG optimal quantizer (for comparison)LBG, same number
of comparatorsproposed algorithm (a) Quantizer performance for a Gaussian
input distribution in $\mathbb{R}^{2}$.

$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$0$$0.1$$0.2$$0.3$$0.4$$0.5$Number of
comparatorsDistortionLBG optimal quantizer (for comparison)LBG, same number of
comparatorsproposed algorithm (b) Quantizer performance for a uniform
unitarian input distribution in $\mathbb{R}^{2}$.

Figure 8:  The red curve shows the distortion of the proposed algorithm and
the dark blue curve shows the distortion from the quantizer of LBG algorithm,
in function of the number of comparators available. For comparison purpose,
the light blue curve shows the performance of a quantizer obtained with LBG
algorithm, but with as many regions as the quantizer obtained with the
proposed algorithm (using more comparators than the proposed algorithm). For
each plot, the values were averaged over numerous restarts of the algorithm.

As one can observe on Fig. 8(a), the quantizer designed by our algorithm
performs better than the quantizer design using LBG algorithm. For a quantizer
using 5 comparators and an input signal in $\mathbb{R}^{2}$, its distortion is
0.64 that of LBG quantizer.  Fig. 7 presents the performance of the output of
LBG algorithm, compared to that of the proposed algorithm. One can observe
there that the quantizer designed by the proposed algorithm performs better
than the quantizer designed by LBG algorithm in every configuration
represented on the figure.

As could be expected, the quantizer obtained with the proposed algorithm
performs slightly worse than the optimal quantizer obtained with LBG algorithm
using the same number of reconstruction points (ie the same output bit rate).
This shows that using the CLVQ framework allows for quantizer designs that
cost little degradation in performance in the classic quantization (where the
constraint is on the cardinality of the output), but generate important
performance gains in the view of comparison-limited quantization.

A result that may be surprising and seems counter-intuitive is that the
hyperplanes do not always form as many regions as possible. Curiously, the
hyperplanes sometimes converge to positions in which they form less regions
than the maximum possible, but with more regular shapes, which implies a
lesser probability of values at a long distance from the representation point,
thus a lower MSE. Fig. 6.b is an example of such configuration. These
configurations turn out to perform relatively well and sometimes better than
configurations with more regions.

####  5.3.2 Uniform distribution

Running the algorithm with a unitary uniform distribution as source
distribution gives similar results as with a Gaussian distribution. In this
case again, our algorithm performs better than LBG with the same number of
comparators. For a quantizer using 5 comparators and an input signal in
$\mathbb{R}^{2}$, its distortion is 0.72 that of LBG quantizer. It also
performs slightly worse than LBG algorithm with the same number of
reconstruction points.

Because of the square shape of the distribution support, the hyperplanes are
even more likely to form rectangular regions with a uniform distribution than
with the Gaussian distribution, as Fig. 4 illustrates.

####  5.3.3 Entropy maximization

As commented in Sec. 3.2, one is often interested in characterizing the
entropy of the quantizer output, as this indicates whether it is necessary to
perform a further compression of the quantizer output. Setting as the
objective as maximization of the entropy of the output, instead of
minimization of the MSE, does not require one to make any significant changes
to the optimization.

In our simulations, the result of the entropy maximization is rather similar
to that of the MSE minimization. The most notable difference in using entropy
as a measure is that in cases like Fig. 6, it favors configurations of type
Fig. 6.a rather than type Fig. 6.b.

##  6 Conclusion

In this paper, a novel paradigm for vector quantization is considered, in
which the performance of the vector quantizer is constrained by the number of
comparators needed to obtain the quantized signal, instead of being
constrained by the output bit rate as in the classic vector quantization
problem. We study the scenario where the vector quantizer is made of $k$
comparators that each receive a linear combination of the quantizer input and
a constant bias, and produce the sign of the received signal as output. We
focus on the task of finding the linear combinations and constant biases,
given $k$, the distribution of the quantizer input and a distortion measure,
so as to minimize the distortion between input and output of the quantizer.

We present an algorithm to solve this optimization problem and show the
performance it attains in the case of mean squared error distortion and
Gaussian and uniform i.i.d. source distributions. The proposed algorithm is
composed of an initialization step and a distinct optimization step. In the
initialization step, a set of possible initial configurations are produced
which span the set of possible solutions. Diverse approaches are presented for
this step: in particular we propose a genetic algorithm which can produce a
rather good initial configuration for a linear complexity. In the optimization
step, the algorithm iteratively optimized the position of the reconstruction
points and the comfiguration of linear combinations received by the
comparators. We perform numerical simulations for this algorithm and compare
the performance attained to that of the classic Linde-Buzo-Gray quantizer. We
show that one can obtain similar performance by using a smaller number of
comparators than the classic approach would require.

A number of research directions remain open from this new vector quantizer
architecture. In particular, we are investigating the optimal performance
attainable in the limit of infinitely long vector quantizer in which the
number of available comparators $k$ and bits available to represent the
quantizer inputs both grow to infinity at a given constant ratio $\alpha$.
This limit should result in a rather interesting generalization of the classic
distortion-rate function.

##  Appendix A – Some notion of matroid theory

Oriented matroids are defined as combinatorial objects that can contain
properties of linear dependence in vector spaces. Matroids have a natural way
of being mapped to some objects such as sets of points or sets of hyperplanes,
while retaining relations of alignment and coplanarity between their elements.
Oriented matroids also have these properties, but add order relations to it.
In consequence, convexity relations inside vector or hyperplane arrangements
are also contained in oriented matroids Bjorner et al. (1994). An oriented
matroid is defined as follows: for a linear functional
$\mathbf{x}^{T}\in\mathbb{R}^{m}$, define the _covector_ of $\mathbf{x}$

$\displaystyle C_{{\cal A}}(\mathbf{x})=\left[{\rm
sign}(\mathbf{a}_{1}^{T}\mathbf{x}-b_{1}),\ldots,{\rm
sign}(\mathbf{a}_{n}^{T}\mathbf{x}-b_{n})\right],$ (24)

and let us define the set of all possible covectors of ${\cal A}$ as

$\displaystyle{\cal L}_{{\cal A}}=\left\\{C_{{\cal A}}(\mathbf{x})\ :\
\mathbf{x}\in\mathbb{R}^{m}\right\\},$ (25)

where the ${\rm sign}$ takes value zero when
$\mathbf{a}_{i}^{T}\mathbf{x}=b_{i}$, so that the support of ${\rm sign}$ is
$\\{-1,0,+1\\}$. The set $M=-([n],{\cal L}_{{\cal A}})$ is defined as the
oriented matroid associated to ${\cal A}$. Using the oriented matroid
framework, one is able to state questions such as the what is the set of
covectors that can be generated through the arrangement of $n$ hyperplanes in
dimension $m$. Regrettably, Vakil showed Vakil (2006) using Mnëv’s
universality theorem Mnëv (1986) that most of the stimulating problems that
could be framed with oriented matroids have no computationally efficient way
to be solved.

##  Appendix B – Graph representation

In order to efficiently describe various combinations of hyperplane
arrangements, we have used graphs to represent them. The object we need to
work with has to be able to represent a hyperplane arrangement’s structure
rather than one particular arrangement. Thus, the object has to be invariant
to rotation, translation and scaling and to any hyperplane modification that
does not change the overall structure of the arrangement. A satisfying
representation is the one shown in Fig. 9.

Figure 9: The left hand side is an arrangement of 3 hyperplanes that form 7
regions, with the grey circle representing infinity. The right hand side is
the corresponding graph (7 vertices for the 7 regions, and one for the
sphere).

The interpretation for Fig. 9 is the following: regions of the hyperplane
arrangement (left side) are individually associated to nodes in the graph
(right side). If two regions in the arrangement share an edge, the related
nodes are connected on the graph (the color of the edge on the graph
corresponds to that of the hyperplane in the arrangement). In consequence, the
number of edges of a region is translated into the degree of the corresponding
node. Projectivization is used in order to accurately represent unbounded
regions. A finite hypersphere is drawn around the zone of interest, and the
space $\mathbb{R}^{m}$ is projected onto that hypersphere so that coordinates
at infinity are mapped to its edge. Then, the edge of the sphere is seen as an
edge of the infinite regions. This way, regions of infinite size can be fully
characterized. Another way to describe this representation is to define it as
the dual graph of the graph that has hyperplane intersections as vertices and
the hyperplanes themselves as edges. The same reasoning for using
projectivization to characterize regions can be applied in a similar way.

Note that graph representation in Fig. 9 is substantially a representation of
the oriented matroid in App. A. Indeed, given a label for each region, these
labels can be mapped to the set of all possible covectors ${\cal L}_{{\cal
A}}$ of the arrangement. In order to obtain such region labels, one can
proceed as follows. First, give an arbitrary orientation to the hyperplanes.
Each hyperplane defines two half-spaces, one of which is attributed a positive
label and the other one a negative label. Then, each region can be attributed
a unique label by concatenating the positive or negative labels given by each
hyperplane (depending on whether the region is on the positive or negative
side of the hyperplane). Thus, each region receives a unique label composed of
a number of binary variables equal to the number of hyperplanes in the
arrangement. Fig. 10 shows an example of such labeling process for an
arrangement of 3 hyperplanes.

Figure 10:  Illustration of the labeling process. In this arrangement, each of
the 3 hyperplanes is oriented: it has a positive side (indicated by an arrow)
and a negative side. Each region is identified in a unique manner by its
location relative to each hyperplane. If it is on the positive side (resp.
negative side) of the $k^{th}$ hyperplane, it has a "+" (resp. "-") in the
$k^{th}$ position of its label.

##  Appendix C – Progressive arrangement growth

By numerical experimentation, it is easy to observe that the optimization step
which relies on random initialization in Sec. 4.1.1 has a high risk of being
stuck in a local minimum. In order to alleviate this problem, we have
considered a second approach to linear comparator initialization that makes
use of the notions described in App. A and of the graph representation
described in App. B. More specifically, (i) we begin by considering a simple
graph with two vertices and one edge, representing two regions separated by
one hyperplane then, recursively, (ii) one generates a new set of graphs by
considering all the possible ways in which a new configuration can be
generated from the old one adding one more hyperplane.

The recursive repetition of step (ii) above results in an exponentially
growing number of configurations. Moreover, this process is susceptible to
produce duplicates or equivalent configurations. To detect and remove
duplicate configurations, one can consider verifying if the graph
representations of two configurations are isomorphic. Verifying graph
isomorphisms is a well studied problem that can be solved in logarithmic time
in the case of planar graphs (which is the case of representation graphs for
arrangements in 2 dimensions). However it is in sub-exponential time for non-
planar graphs (and thus for higher dimension arrangements). This presents an
obstacle to the problem of determining all possible hyperplane configurations.
More generally, there exist two other obstacles that prevent one from
generating all possible initial configurations and optimize each of them.
First, the task of generating all these configurations is computationally very
heavy. This is even more true in higher dimensions, as complexity grows super-
exponentially with the number of dimensions. Second, the task of converting a
hyperplane arrangement structure (that is a graph representation in our case)
into an actual hyperplane arrangement with numerical values is highly complex,
as is implied by Mnëv’s universality theorem Mnëv (1986). The high
computational complexity of these two steps cause the approach not to be
applicable in practice.

The work of S. Rini is supported by the MOST grant number 109 -2927-I-009
-512.

###### Acknowledgements.

We would like to give thanks to Arnau Padrol from Sorbonne Université for his
help on understanding oriented matroids. The authors declare no conflict of
interest. The following abbreviations are used in this manuscript:
MSE | Mean Squared Error
---|---
LBG | Linde-Buzo-Gray algorithm
MIMO | Massive Input Massive Output
GP | General Position
A2D | Analog-to-Digital
VQ | Vector Quantization
CLVQ | Comparison-Limited Vector Quantization
References yes

## References

  * Gray (1984) Gray, R.  Vector quantization.  IEEE Assp Magazine 1984, 1, 4–29.
  * Gersho and Gray (2012) Gersho, A.; Gray, R.M.  Vector quantization and signal compression; Vol. 159, Springer Science & Business Media, 2012.
  * Makhoul et al. (1985) Makhoul, J.; Roucos, S.; Gish, H.  Vector quantization in speech coding.  Proceedings of the IEEE 1985, 73, 1551–1588.
  * Paliwal and Atal (1993) Paliwal, K.K.; Atal, B.S.  Efficient vector quantization of LPC parameters at 24 bits/frame.  IEEE transactions on speech and audio processing 1993, 1, 3–14.
  * Yahampath and Rondeau (2007) Yahampath, P.; Rondeau, P.  Multiple-description predictive-vector quantization with applications to low bit-rate speech coding over networks.  IEEE transactions on audio, speech, and language processing 2007, 15, 749–755.
  * Gersho and Ramamurthi (1982) Gersho, A.; Ramamurthi, B.  Image coding using vector quantization.  ICASSP’82. IEEE International Conference on Acoustics, Speech, and Signal Processing. IEEE, 1982, Vol. 7, pp. 428–431.
  * Nasrabadi and King (1988) Nasrabadi, N.M.; King, R.A.  Image coding using vector quantization: A review.  IEEE Transactions on communications 1988, 36, 957–971.
  * Chiranjeevi and Jena (2018) Chiranjeevi, K.; Jena, U.R.  Image compression based on vector quantization using cuckoo search optimization technique.  Ain Shams Engineering Journal 2018, 9, 1417–1431.
  * Lee and Woods (1995) Lee, Y.Y.; Woods, J.W.  Motion vector quantization for video coding.  IEEE transactions on image processing: a publication of the IEEE Signal Processing Society 1995, 4, 378–382.
  * Regunathan et al. (2014) Regunathan, S.; Sun, S.; Tu, C.; Lin, C.L.  Adaptive quantization for enhancement layer video coding, 2014.  US Patent 8,897,359.
  * Soong et al. (1987) Soong, F.K.; Rosenberg, A.E.; Juang, B.H.; Rabiner, L.R.  Report: A vector quantization approach to speaker recognition.  AT&T technical journal 1987, 66, 14–26.
  * Hasan et al. (2004) Hasan, M.R.; Jamil, M.; Rahman, M.; others.  Speaker identification using mel frequency cepstral coefficients.  variations 2004, 1.
  * Wu and Chang (2005) Wu, H.C.; Chang, C.C.  A novel digital image watermarking scheme based on the vector quantization technique.  Computers & Security 2005, 24, 460–471.
  * Lu et al. (2005) Lu, Z.M.; Xu, D.G.; Sun, S.H.  Multipurpose image watermarking algorithm based on multistage vector quantization.  IEEE Transactions on Image Processing 2005, 14, 822–831.
  * Equitz (1989) Equitz, W.H.  A new vector quantization clustering algorithm.  IEEE transactions on acoustics, speech, and signal processing 1989, 37, 1568–1575.
  * Lughofer (2008) Lughofer, E.  Extensions of vector quantization for incremental clustering.  Pattern recognition 2008, 41, 995–1011.
  * Berger (2003) Berger, T.  Rate-distortion theory.  Wiley Encyclopedia of Telecommunications 2003.
  * Shannon (1948) Shannon, C.E.  A mathematical theory of communication.  The Bell system technical journal 1948, 27, 379–423.
  * Mezghani and Nossek (2007) Mezghani, A.; Nossek, J.A.  On ultra-wideband MIMO systems with 1-bit quantized outputs: Performance analysis and input optimization.  2007 IEEE International Symposium on Information Theory. IEEE, 2007, pp. 1286–1289.
  * Sundstrom et al. (2008) Sundstrom, T.; Murmann, B.; Svensson, C.  Power dissipation bounds for high-speed Nyquist analog-to-digital converters.  IEEE Transactions on Circuits and Systems I: Regular Papers 2008, 56, 509–518.
  * Mo and Heath (2014) Mo, J.; Heath, R.W.  High SNR capacity of millimeter wave MIMO systems with one-bit quantization.  2014 Information Theory and Applications Workshop (ITA). IEEE, 2014, pp. 1–5.
  * Rini et al. (2017) Rini, S.; Barletta, L.; Eldar, Y.C.; Erkip, E.  A general framework for MIMO receivers with low-resolution quantization.  2017 IEEE Information Theory Workshop (ITW). IEEE, 2017, pp. 599–603.
  * Khalili et al. (2018) Khalili, A.; Rini, S.; Barletta, L.; Erkip, E.; Eldar, Y.C.  On MIMO channel capacity with output quantization constraints.  2018 IEEE International Symposium on Information Theory (ISIT). IEEE, 2018, pp. 1355–1359.
  * Shlezinger et al. (2018) Shlezinger, N.; Eldar, Y.C.; Rodrigues, M.R.  Hardware-limited task-based quantization.  arXiv preprint arXiv:1807.08305 2018.
  * Shlezinger and Eldar (2019) Shlezinger, N.; Eldar, Y.C.  Deep task-based quantization.  arXiv preprint arXiv:1908.06845 2019.
  * Shohat et al. (2019) Shohat, M.; Tsintsadze, G.; Shlezinger, N.; Eldar, Y.C.  Deep Quantization for MIMO Channel Estimation.  ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2019, pp. 3912–3916.  doi:black10.1109/ICASSP.2019.8682704.
  * Khobahi et al. (2019) Khobahi, S.; Naimipour, N.; Soltanalian, M.; Eldar, Y.C.  Deep Signal Recovery with One-bit Quantization.  ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2019, pp. 2987–2991.  doi:black10.1109/ICASSP.2019.8683876.
  * Shlezinger et al. (2019) Shlezinger, N.; Eldar, Y.C.; Salamatian, S.; Médard, M.  Task-based quantization for recovering quadratic functions using principal inertia components 2019.
  * Shlezinger et al. (2019) Shlezinger, N.; Eldar, Y.C.; Rodrigues, M.R.D.  Asymptotic Task-Based Quantization With Application to Massive MIMO.  IEEE Transactions on Signal Processing 2019, 67, 3995–4012.  doi:black10.1109/TSP.2019.2923149.
  * Kamilov et al. (2012) Kamilov, U.S.; Bourquard, A.; Amini, A.; Unser, M.  One-bit measurements with adaptive thresholds.  IEEE Signal Processing Letters 2012, 19, 607–610.
  * Plan and Vershynin (2013) Plan, Y.; Vershynin, R.  One-Bit Compressed Sensing by Linear Programming.  Communications on Pure and Applied Mathematics 2013, 66, 1275–1297.
  * Khobahi et al. (2019) Khobahi, S.; Naimipour, N.; Soltanalian, M.; Eldar, Y.C.  Deep signal recovery with one-bit quantization.  ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019, pp. 2987–2991.
  * Max (1960) Max, S.  Quantizing for minimum distortion.  IEEE transactions on information theory 1960, 6, 7–12.
  * Lloyd (1982) Lloyd, S.  Least squares quantization in PCM.  IEEE transactions on information theory 1982, 28, 129–137.
  * Linde (1980) Linde, Buzo, G.  An Algorithm for Vector Quantizer Design.  IEEE transactions on Communications 1980, 28, 84–95.
  * Cover (1999) Cover, T.M.  Elements of information theory; John Wiley & Sons, 1999.
  * Conway and Sloane (1982) Conway, J.; Sloane, N.  Fast quantizing and decoding and algorithms for lattice quantizers and codes.  IEEE Transactions on Information Theory 1982, 28, 227–232.
  * Zamir and Feder (1992) Zamir, R.; Feder, M.  On universal quantization by randomized uniform/lattice quantizers.  IEEE Transactions on Information Theory 1992, 38, 428–436.
  * Stanley (2006) Stanley, R.P.  An introduction to Hyperplane Arrangements; 2006.
  * Dimca and PARUSINSKI (2017) Dimca, A.; PARUSINSKI, A.  Hyperplane arrangements; Springer, 2017.
  * Nelson and Gailly (1995) Nelson, M.; Gailly, J.L.  The data compression book 2nd edition.  M & T Books, New York, NY 1995.
  * Chou et al. (1989) Chou, P.A.; Lookabaugh, T.; Gray, R.M.  Entropy-constrained vector quantization.  IEEE Transactions on Acoustics, Speech, and Signal Processing 1989, 37, 31–42.
  * Orlik and Terao (2013) Orlik, P.; Terao, H.  Arrangements of hyperplanes; Vol. 300, Springer Science & Business Media, 2013.
  * (44) Chataignon, J.  Comparison-limited vector quantizer design algorithm.  doi:black10.5281/zenodo.3552797.
  * Bjorner et al. (1994) Bjorner, A.; Las Vergnas, M.; Sturmfels, B.; White, N.; Ziegler, G.  Oriented Matroids; 1994.
  * Vakil (2006) Vakil, R.  Murphy’s law in algebraic geometry: Badly-behaved deformation spaces.  Inventiones mathematicae 2006, 164, 569–590.
  * Mnëv (1986) Mnëv, N.  The topology of configuration varieties and convex polytopes varieties.  PhD thesis, Stanford University, 1986.