ABSTRACT
PILGRIM, RYAN ZACHARY. Source Coding Optimization for Distributed Average Consensus.
(Under the direction of Dror Baron.)
Consensus is a common method for computing a function of the data distributed among
the nodes of a network. Of particular interest is distributed average consensus, whereby the

arXiv:1710.01816v1 [cs.IT] 4 Oct 2017

nodes iteratively compute the sample average of the data stored at all the nodes of the network
using only near-neighbor communications. In real-world scenarios, these communications
must undergo quantization, which introduces distortion to the internode messages. In this
thesis, a model for the evolution of the network state statistics at each iteration is developed
under the assumptions of Gaussian data and additive quantization error. It is shown that
minimization of the communication load in terms of aggregate source coding rate can be
posed as a generalized geometric program, for which an equivalent convex optimization
can efficiently solve for the global minimum. Optimization procedures are developed for
rate-distortion-optimal vector quantization, uniform entropy-coded scalar quantization, and
fixed-rate uniform quantization. Numerical results demonstrate the performance of these
approaches. For small numbers of iterations, the fixed-rate optimizations are verified using
exhaustive search. Comparison to the prior art suggests competitive performance under certain
circumstances but strongly motivates the incorporation of more sophisticated coding strategies,
such as differential, predictive, or Wyner-Ziv coding.

© Copyright 2017 by Ryan Zachary Pilgrim
All Rights Reserved

Source Coding Optimization for Distributed Average Consensus

by
Ryan Zachary Pilgrim

A thesis submitted to the Graduate Faculty of
North Carolina State University
in partial fulfillment of the
requirements for the Degree of
Master of Science

Electrical Engineering
Raleigh, North Carolina
2017

APPROVED BY:

David Lalush

Michael Daniele

Dror Baron
Chair of Advisory Committee

DEDICATION
To everyone who has stood by me in tough times: my parents, my brothers, and my girlfriend,
Kristin. My success rests on a foundation of love and support from you all, and I am grateful.

ii

BIOGRAPHY
The author was born in Boulder, Colorado in 1991. His childhood was spent in Colorado,
Virginia, Mississippi, and North Carolina. After graduating from Hickory High school, he
attended the University of North Carolina at Asheville. He then transferred to North Carolina
State University to study Biomedical Engineering, graduating magna cum laude with a Bachelor
of Science in 2014.
His research interests span statistical signal processing, compressive sensing, and machine
learning. He also enjoys programming and embedded systems development. In his spare time,
he plays disc golf and the guitar.
The author completed an internship at Southwest Research Institute in 2015, where he
helped create a GPS spoofing detection system. He is currently a Project Engineer at Vadum,
Inc. in Raleigh, North Carolina.

iii

ACKNOWLEDGEMENTS
I would like to thank my advisor, Dror Baron, for his patient mentorship and excellent advice.
Thanks also to my collaborators Waheed Bajwa and Junan Zhu for their many helpful inputs
regarding consensus and the research process. I would furthermore like to acknowledge
Yanting Ma for helping to expand the main theorem of the thesis to the more general, variable
distortion setting.
Special thanks to Mehmet Ercan Yildiz and Anna Scaglione for graciously providing their
code for comparison, and to Yaoqing Yang and Pulkit Grover for helpful discussions about
their information theoretic work with Soummya Kar on the consensus problem.

iv

TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
Chapter 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . .
1.1 Summary of prior art . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 Static, dynamic, and dithered quantization schemes
1.1.3 Topology and weight tuning . . . . . . . . . . . . . .
1.1.4 Sequence filtering . . . . . . . . . . . . . . . . . . . . .
1.1.5 Wireless considerations . . . . . . . . . . . . . . . . .
1.1.6 Information theoretic approaches . . . . . . . . . . .
1.2 Motivation and contributions . . . . . . . . . . . . . . . . . .
1.3 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Notation and acronyms . . . . . . . . . . . . . . . . . . . . .
1.4.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.2 Acronyms . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

. . 1
.
2
.
2
.
3
.
5
.
6
.
6
.
7
.
8
.
9
.
10
.
10
.
12

Chapter 2 QUANTIZED DISTRIBUTED AVERAGE CONSENSUS
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Mathematical preliminaries . . . . . . . . . . . . . . . . . . . . .
2.2.1 Matrix and graph theory concepts . . . . . . . . . . . . .
2.3 The consensus problem . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Lossy consensus . . . . . . . . . . . . . . . . . . . . . . .
2.4 Source coding, quantization, and rate-distortion theory . . . . .
2.4.1 Coding and information theory . . . . . . . . . . . . . .
2.4.2 Rate-distortion theory . . . . . . . . . . . . . . . . . . . .
2.4.3 Uniform and nonuniform quantization . . . . . . . . . .
2.4.4 The quantization noise model and dithering . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
14
.
14
.
15
.
15
.
16
.
18
. . 21
. . 21
.
22
.
26
.
27

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
30
. . 31
.
33
.
36
.
37
.
37
.
39
.
45
.
47
.
50
. . 51
.
54

Chapter 3 CODING RATE OPTIMIZATION IN CONSENSUS
3.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . .
3.2 Abstractions for vector-valued data . . . . . . . . . . . . .
3.3 Assumptions and analytical results . . . . . . . . . . . . . .
3.3.1 Additive noise model . . . . . . . . . . . . . . . . .
3.3.2 Preservation of Gaussian distribution . . . . . . . .
3.3.3 State evolution for vector consensus . . . . . . . . .
3.4 Optimization via generalized geometric programming . .
3.4.1 Cost function formulation . . . . . . . . . . . . . . .
3.4.2 Basics of generalized geometric programming . . .
3.4.3 Proof that the optimization is a GGP . . . . . . . .
3.5 Efficient implementation of GGP optimization . . . . . . .

v

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

3.5.1
3.5.2

Equal-distortion simplification and implementation . . . . . . . . . . . .
Efficient search heuristic for fixed-rate quantizers . . . . . . . . . . . . .

Chapter 4 NUMERICAL RESULTS . . . . . . . . . . .
4.1 Structure of solutions . . . . . . . . . . . . . . . .
4.1.1 Nonmonotone behavior . . . . . . . . . .
4.2 Comparison to prior art . . . . . . . . . . . . . .
4.2.1 Implementation details . . . . . . . . . .
4.2.2 Discussion of comparison results . . . . .
4.2.3 Comparison to predictive coding scheme
4.2.4 Comparison to the ProgQ scheme . . . .
Chapter 5

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

55
56

.
58
.
60
.
62
.
70
. . 71
.
72
.
75
.
78

CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

vi

LIST OF TABLES

vii

LIST OF FIGURES
The RD function RpDq of a Gaussian memoryless source. The function is
a nonincreasing, convex function of the normalized distortion σD2 , where
σ2 is the variance of the Gaussian source. The RD pairs above the curve
are achievable, while those below the curve are not. . . . . . . . . . . . .
Figure 2.2 Operational RD relationships of practical quantizers: ECSQ, Lloyd-Max,
and fixed-rate uniform quantization. Note that all of these approach
an asymptote that is a constant offset from the Gaussian RD function
as R Ñ 8. The ECSQ performance corresponds to the uniform scalar
quantizer with a level at zero. The representation level placement may be
suboptimal, but it allows for coding rates below one bit per symbol. This
figure was inspired by a similar plot in lecture notes from a TU Berlin
source coding course [56]. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.1

Figure 3.1

Surface showing the evolution of MSEi ptq for the case of lossless consensus. Note that the value of MSEi ptq is not equal across all nodes for a
particular iteration, and the evolution of MSEi ptq is heavily dependent on
network topology. Also observe that the MSE decreases monotonically
at each node as t Ñ 8. The variance νi ptq evolves in a similar fashion. .

23

25

46

Optimal rates and MSE sequences from the solution of (3.38) (T “ 5,
ρc “ 0.35, σx2 “ 1, σn2 “ 0.5, m “ 20). Top: Optimal rate sequences for the
variable-distortion optimization problem. The rates are plotted against
iteration indices, and each line represents the rates used by a different
sensor. Bottom: MSE sequence corresponding to the above rate sequence.
The MSE values are plotted against node indices, and each line represents
a different iteration. In this case, the MSE values decrease monotonically
for all nodes. The lines represent values predicted by the state evolution
equations (3.23)–(3.26), and the overlaid circles represent the simulated
values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
Figure 4.2 Optimal rates and MSE sequences from the solution of (3.39) (T “ 5,
ρc “ 0.35, σx2 “ 1, σn2 “ 0.5, m “ 20). Top: Optimal rate sequences for the
constant-distortion optimization problem. The rates are plotted against
iteration indices, and each line represents the rates used by a different
sensor. Bottom: MSE sequence corresponding to the above rate sequence.
The MSE values are plotted against node indices, and each line represents
a different iteration. In this case, the MSE values decrease monotonically
at all nodes. The lines represent values predicted by the state evolution
equations (3.23)–(3.26), and the overlaid circles represent the simulated
values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Figure 4.1

viii

Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

Influence of connectivity on the optimal rates from the solution of (3.39)
(T “ 5, ρc “ 0.35, σx2 “ 1, σn2 “ 0.5, m “ 20). Top: Scatter plot of aggregate
rate used at each node versus node degree. Bottom: Surface of optimal
rate sequences at each node, where the node indices are sorted by node
degree in decreasing order. . . . . . . . . . . . . . . . . . . . . . . . . . .
Nonmonotone optimal rate sequences for the low-SNR, low-MSE setting
(T “ 5, ρc “ 0.35, σx2 “ 1, σn2 “ 20, m “ 20). Each line corresponds
to the rate sequence used by a different node. Top: Variable-distortion
optimization (3.38) result. Bottom: Constant-distortion optimization (3.39)
result. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Nonmonotone MSE sequences. Each line corresponds to the MSE sequence of a different node (T “ 5, ρc “ 0.35, σx2 “ 1, σn2 “ 0.5). Top:
MSE evolution for the variable-distortion problem (3.38). Bottom: MSE
evolution for the constant-distortion problem (3.39). . . . . . . . . . . .
Coefficients of the noise variance σn2 (solid blue curve) and distortion
Dp0q (dashed red curve) in the equation for the MSE (3.31) when the
states are initialized as a Gaussian random vector corrupted by AWGN
(ρc “ 0.35, m “ 20). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
RD trade-off curves for the proposed schemes versus order-one predictive
coding [21]. All schemes were simulated for random geometric graphs on
a unit torus with connectivity radius ρc “ 0.35. Top: number of iterations
T “ 5. Bottom: number of iterations T “ 7. . . . . . . . . . . . . . . . . .
RD trade-off curves for the proposed schemes versus order-one predictive
coding [21]. All schemes were simulated for random geometric graphs on
a unit torus with connectivity radius ρc “ 0.45. Top: number of iterations
T “ 5. Bottom: number of iterations T “ 7. . . . . . . . . . . . . . . . . .
RD trade-off curves for the proposed schemes versus ProgQ [20]. All
schemes were simulated for random geometric graphs on a unit torus
with connectivity radius ρc “ 0.35. Top: number of iterations T “ 5.
Bottom: number of iterations T “ 7. . . . . . . . . . . . . . . . . . . . . .
RD trade-off curves for the proposed schemes versus ProgQ [20]. All
schemes were simulated for random geometric graphs on a unit torus
with connectivity radius ρc “ 0.45. Top: number of iterations T “ 5.
Bottom: number of iterations T “ 7. . . . . . . . . . . . . . . . . . . . . .

ix

63

64

66

69

73

74

76

77

CHAPTER

1
INTRODUCTION

The wide availability of wireless sensors and large data sets in recent years has provided
significant motivation for the development of distributed computing methods. Large sensor
networks and big, distributed data sets pose a number of unique optimization opportunities, including energy management in battery-powered sensor networks [1, 2] and run time
reduction in cloud computing applications [3].
In many distributed computing problems, it is necessary to compute a function of data
that may be dispersed among a number of computing nodes, for instance in sensor networks,
where each agent observes a different measurement of a physical process, or large-scale
server farms, where the size of the data set requires distributed storage [4]. The class of
distributed algorithms considered in this thesis computes these functions with only local, nearneighbor interactions. This distributed approach offers several potential advantages, including
robustness to link and node failure and maintenance of privacy [5]. One popular approach
to distributed function computation is called consensus [5, 6], which has many variants.
Consensus protocols have found application in a wide variety of settings, including oscillator
synchronization, robot swarm control, load balancing in computer networks, distributed
sensor fusion, belief propagation [6], detection, estimation, target tracking [5], distributed

1

optimization, filtering, and environmental monitoring [7]. Consensus algorithms have also
found applications in computer vision [8], where they can be used for null-space and leastsquares estimation, distributed singular value decomposition, principal components and
generalized principal components analyses, point triangulation, and linear pose estimation.
Furthermore, consensus algorithms have recently found application in distributed eigenvector
estimation and dictionary learning [9].
Although it has a wide variety of engineering applications, consensus originated as a
model of social interaction in the management science and statistics communities [10] during
the 1960s, according to Olfati-Saber et al. [6]. Its application to the problem of distributed
networked computation was first studied by Borkar and Varaiya [11] and Tsitsiklis [12, 13],
where it was applied to distributed estimation, optimization, and decision making. Much of the
early work on consensus assumed that the nodes could communicate real-valued data to one
another [14]. In realistic scenarios, the nodes must communicate within bandwidth and energy
constraints, which can have a significant impact on the convergence of distributed averaging
algorithms. Although many papers have been published on quantized consensus in recent
years [4, 7, 14–28], a large portion consider trade-offs among run time, communication load,
and final accuracy in ways that do not take advantage of the tools of constrained optimization.
This work attempts to take a more structured approach to the reduction of communication
in consensus than past efforts by considering rate-distortion (RD) theory [29, 30] and convex
optimization techniques [31].

1.1
1.1.1

Summary of prior art
Overview

Many previous works have addressed consensus with limited coding rate. Due to the multidisciplinary nature of the subject, several approaches have been taken to address the issue of
transmitting limited-precision information. These include deterministic and dithered quan-

2

tization strategies using both static and dynamic quantization schemes. Because the classic
linear-update consensus algorithms suffer from divergence in the presence of quantized or
noisy exchanges [32], many of these works focus on proposing new algorithms and proving
their asymptotic convergence properties. A variety of methods has been employed to sidestep
the many issues quantization poses, including adapting the quantization range, tuning the
weight link sequence, using dithering or randomized quantization, and filtering the past and
present state values. This section provides a survey of these works before introducing and
distinguishing our approach.

1.1.2

Static, dynamic, and dithered quantization schemes

Early publications on the topic, such as the work of Xiao et al. [32], show that introducing
perturbations of constant variance (such as quantization error) into the traditional consensus
state update prevents convergence due to the limited precision of the quantizer. Similar
convergence issues preclude convergence when the quantization range is held constant [14].
Chamie et al. [16], like Frasca et al. [14], considered the case of nonadaptive, deterministically
quantized consensus, and they showed that in finite time, consensus is achieved in the sense
that the network state converges to one of the quantization levels.
To surmount the convergence issue associated with quantization, Aysal et al. [19] proposed
a randomized quantization scheme, and they showed it was equivalent to dithering. With this
probabilistic quantization scheme, Aysal et al. proved that consensus can be achieved almost
surely at one of the neighboring quantization values and analyzed its performance relative to
unquantized consensus. In the case of networks with quantization and random link failures
(i.e., time-varying topology), Kar and Moura [27] showed that by subjecting the communication
links to a persistence condition, dithered quantized consensus can converge. They also derived
probability bounds on consensus within a certain mean square error (MSE) range. Under
similar topological assumptions, Nedić et al. [25] assessed the performance degradation
resulting from quantization and provided tight polynomial bounds on the convergence time

3

for a wide class of distributed averaging algorithms.
Due to the difficulties associated with quantization error, many works address the incorporation of dynamic encoding/decoding strategies into consensus protocols. However, many
of these schemes do not explicitly consider the RD trade-off and offer certain heuristics to
optimize communication performance within their proposed frameworks.
Carli et al. [15, 17] assessed the performance of a “zoom-in, zoom-out” strategy originally
studied in the context of quantized feedback control system stabilization. In this scheme, the
quantization range grows in the case of saturation and shrinks otherwise. The authors demonstrated convergence for certain topologies, numbers of quantizer levels, and initial quantizer
range values. Interestingly, they showed that under certain circumstances, convergence speed
can be faster than in the case of ideal, unquantized transmission.
Similarly to Carli et al. [15, 17], Rego et al. [33] proposed an algorithm with progressively
shrinking quantization range and derived conditions on the design parameters to guarantee
bounded steady-state error. To demonstrate the efficacy of their approach, they simulated a
vehicle formation control problem.
Li et al. [34] proposed a control scheme capable of guaranteeing asymptotic consensus with
as little as one-bit data exchange per iteration. The authors quantified the convergence rate
in terms of network properties and coding rate. To achieve convergence while transmitting a
single bit at each node per iteration, the authors relied on differential encoding/decoding [35,
Sec. 7.2]. Fang and Li [18] proposed a scheme that adjusts its quantizer parameters by learning
from previous runs, and they showed convergence under the assumption that the quantization
error converges to zero. Similarly, Thanou et al. [20] proposed a differential encoding strategy
with exponentially decaying quantization range. Their work is based on an “average-case”
analysis that does not provide the same strict performance guarantees as Li et al. [34], but
allows for better typical performance. The numerical results demonstrate lower MSE than a
number of previous works (including Li et al. [34], Fang and Li [18], and Carli et al. [15, 17])
with equal communication load, and an error decay that matches the rate of the unquantized

4

algorithm in the limit of many iterations. Both of these works assume a constant coding rate
throughout all iterations and rely on differential encoding/decoding to achieve convergence.
Yildiz and Scaglione [21, 22, 36], unlike other authors, explicitly considered the RD tradeoff to achieve an asymptotic MSE value in consensus with Gaussian states. They proposed
schemes based on differential [22], predictive, and Wyner-Ziv coding [21, 36]. In Yildiz and
Scaglione [21], the authors exploited the correlation of the current network state with that
of previous iterations, and showed that bounded steady-state error is possible under their
schemes using shrinking coding rates. Modeling the quantization error as an additive noise,
they also provided a necessary and sufficient condition on the variance of the quantization
error to ensure convergence. The main focus of the work [21] is to show that convergence can be
achieved using asymptotically decreasing coding rates under predictive and Wyner-Ziv coding
schemes. In [22], Yildiz and Scaglione considered a simpler case of differential coding and
showed similar results for this approach. Yildiz and Scaglione [22] also imposed a parametric
form on the distortion sequence and examined the effect of varying the convergence rate on
the aggregate coding rate required to achieve a target asymptotic MSE.

1.1.3

Topology and weight tuning

Another approach to combat the effect of quantization and channel noise, and to improve
energy consumption, is to tune the network topology or estimation update rule. In many of
the most popular consensus algorithms, the update rule consists of a convex combination of
the current node state with those of its neighbors. By adapting which nodes are considered
neighbors or by altering the weights associated with each node in the neighborhood, it is
possible to improve convergence properties in the presence of quantization. In the absence of
quantization, Xiao and Boyd [37] showed that the fastest linear iterative update can be found
by semidefinite programming. They also listed a number of heuristics for weight design in
the case of incomplete knowledge of the network topology. In their later work, Xiao et al. [32]
considered the problem of noisy exchanges and showed that the weights corresponding to the

5

smallest steady-state error can be found via convex optimization.
Mosquera et al. [23] considered a greedy approach to updating the weight sequence by
minimizing the minimum MSE (MMSE) of the estimates at each node during each iteration
and proposed a modified scheme that only requires statistical knowledge about the topology.
The modified scheme approximates a random geometric network by a regular graph, for
which each node has the same number of neighbors.

1.1.4

Sequence filtering

To suppress the perturbations resulting from quantization, some authors have considered
the possibility of using the history of past state values. Zhu et al. [26, 38] considered the
problem of distributed parameter estimation with quantization error, proposed a scheme
to reduce randomness using a moving average, and bounded its almost sure performance.
Similarly, Thanou et al. [39] considered the problem of quantized distributed averaging and
demonstrated a technique to find the optimal polynomial filter coefficients to minimize the
effect of quantization error. Fang and Li [40] developed a sequence averaging approach with
convergence properties that improved over Frasca et al. [14]. In contrast to Zhu et al. and
Thanou et al., Fang and Li only computed the average in the final iteration. These approaches
both reduce the randomness introduced by quantization and accelerate the convergence of
consensus algorithms.

1.1.5

Wireless considerations

In cases where the network is wireless, certain aspects of the communication medium can
be exploited. Nokleby et al. [7] assessed the resource consumption of consensus under a
wireless path-loss model. They considered total transmit energy, elapsed time, and timebandwidth product and showed that by recursively forming geographic clusters, consensus
can be achieved with nearly order-optimal performance relative to all three metrics. For
comparison, they also assessed the performance of some popular gossip algorithms. Similarly,

6

Huang and Hua [28] designed an energy planning algorithm for progressive estimation and
consensus in multihop wireless sensor networks (WSNs). They formulated energy models
based on assumptions on the wireless channels, and proposed energy optimization approaches
for both types of estimation in the presence of quantization. However, Huang and Hua [28]
only considered the simple fixed-rate uniform quantizer and did not allow the coding rate to
vary over the iterations or nodes in their analysis of consensus.

1.1.6

Information theoretic approaches

Although much of the literature considers the design of specific protocols for consensus
averaging, a handful of works have explored the fundamental RD limits in the distributed
computation problem. Two of these derive bounds on computation time using information
theoretic inequalities [41, 42], but they do not consider RD theory [29, 30] in their analyses.
Recently, Su and El Gamal [43] considered the problem of computing the RD function for
distributed average consensus using peer-to-peer communication protocols. The authors
derived a closed form for the RD function for a two-node network, and derived upper and
lower bounds on the RD function for weighted-sum and gossip-based protocols in arbitrary
networks. However, the authors restricted their attention to a class of protocols that uses
communication between node pairs and assumes time-invariant normalized distortion. The
analysis of this thesis, by contrast, assumes communication among more than one agent
at each time step (corresponding to broadcast, rather than peer-to-peer, protocols), greater
flexibility in the selection of edge weights, and time-varying distortions. Furthermore, Su and
El Gamal did not present results on the tightness of their weighted-sum or gossip bounds,
nor did they demonstrate the achievability of these bounds except in the case of a large star
network with a centralized protocol.
In a similar vein, Yang et al. [44] considered RD bounds for two lossy in-network function
computation scenarios. The first of these is data aggregation, in which data is routed through
a tree network to a fusion center. Along the way, each node computes a partial result before

7

communicating its result to the next node in the path. The second scenario is that of consensus,
in which each node forms an estimate of the desired quantity. Although Yang et al. [44]
provided bounds on the RD relationship for consensus in tree networks and proved the
achievability of the derived bounds, their analysis is limited to the setting where the network
is tree-structured. Often it is beneficial to consider more flexible topologies, such as random
geometric graphs [45], which have been used to model WSNs [46]. In general, random
geometric graphs and their real-world WSN counterparts have loops.

1.2

Motivation and contributions

This thesis presents a framework for attaining an estimate of the network sample mean at
each node, within a desired average level of accuracy, with finite run time and minimal
total communication cost using either deterministic or dithered quantization. Our analysis is
informed by the results of RD theory [29, 30] and convex optimization [31], which allow a more
structured approach than some of the prior studies in the literature that neglect constrained
optimization. The proposed cost function, which originally appeared in the work of Zhu
and coauthors [47–51] is simple but capable of modeling diverse networks, from large-scale
server farms and cloud services [3] to battery-powered WSNs [1, 2]. Unlike much of the prior
art, which is focused on showing asymptotic convergence for particular protocols, this thesis
takes a constrained optimization approach similar to Huang and Hua [28] to optimizing the
communication scheme for finite run time, subject to an accuracy constraint. Unlike Huang
and Hua [28], however, we allow time- and node-varying coding rates. We explicitly assess
the trade-offs among run time, total communication load as measured by aggregate source
coding rate, and quantizer complexity to attain a final MSE across the network. The focus of
this work is on a Gaussian-distributed initial network state, but the results can also be used
to design fixed-rate uniform quantization schemes for other distributions. Additionally, for
the case of fixed-rate coding, we present a heuristic that reduces optimization complexity but
scales well compared to the exact problem. The accuracy of this heuristic is verified using

8

numerical experiments.
The advantages of our approach are (i) ignorance about the parametric form of the
distortions, which allows greater flexibility in the selection of the optimal rates, (ii) support for
different rates at each node and iteration of the algorithm, and (iii) optimization with respect
to exact MSE quantities for finite iteration count.
Point (i) distinguishes our approach from Yildiz and Scaglione [22], who restricted their
study to optimal distortion sequences that formed convergent series. Point (ii) contrasts this
thesis with the work of Huang and Hua [28] and Thanou et al. [20], who required the use of a
fixed-rate uniform quantizer with a rate that was constant over both nodes of the network and
iterations of the algorithm. Point (iii) differentiates this thesis from both Huang and Hua [28]
and Yildiz and Scaglione [21, 22]; Huang and Hua [28] optimized with respect to a bound on
the MSE, and Yildiz and Scaglione [21, 22] considered the asymptotic MSE as the number of
iterations goes to infinity.

1.3

Organization

The remainder of this thesis is organized as follows. Chapter 2 presents the previous results
in the literature, including the basics of consensus, quantization, rate-distortion theory, the
additive quantization noise model, and dithering. Therefore, Chapter 2 does not contain
any original results. In Chapter 3, we present the main contributions of this thesis, namely
the state-evolution model for the modified consensus iteration of Frasca et al. [14] and the
application of generalized geometric programming (GGP) to optimize the source coding in
consensus. Chapter 4 presents numerical results demonstrating the performance of the GGP
approach and compares its performance to some of the prior art. Finally, Chapter 5 includes
some concluding remarks on the contributions of this thesis, the merits and drawbacks of the
presented approach, and possible topics for future research.

9

1.4
1.4.1

Notation and acronyms
Notation

In this thesis, uppercase bold letters (e.g., A) will be used to denote matrices, and lowercase
bold letters (e.g., x) will be used to denote vectors. Vectors that vary with time are assigned
a time index, so that a time-varying x becomes xptq. Scalars that vary with time are denoted
similarly (e.g., xptq). A superscript on a square matrix (e.g., Ak , A P Rnˆn , k a positive integer)
denotes raising that matrix to the kth power (i.e, multiplying that matrix by itself k ´ 1 times).
Random variables (RVs) are not distinguished by notation to avoid confusing vectors and
matrices.
The following list enumerates frequently used notation and variables and their associated
meanings. The meaning of the following notation may not be clear until later, and it is provided
here for the reader’s reference.
• t¨uJ : transpose
• Qp¨q: quantization function
• I: identity matrix
• 0: the matrix or vector of all zeros
• 1: vector of all ones
• R: real numbers
• Z: integers
• Zě0 : nonnegative integers
• Zą0 : positive integers
• E r¨s: statistical expectation

10

• N pµ, Σq: multivariate Gaussian distribution with mean µ and covariance Σ
• k¨k p : ` p norm
• G : undirected graph
• V : vertex (node) set of a graph
• E : edge set of a graph
• Ni : neighborhood of node i
• X , Y : finite subsets of R
• Hpxq: entropy of the random variable x
• Ipx; yq: mutual information between the random variables x and y
• ppxq: probability mass function (PMF) of the discrete random variable x
• f pxq: probability distribution function (pdf) of the continuous random variable x
• Φ x pωq: characteristic function of the random variable x
• U pa, bq: uniform distribution with minimum a and maximum b
• r¨si : ith component of a vector
• r¨sij : pi, jqth component of a matrix
• µx ptq: mean of the random vector xptq
• Σx ptq: covariance matrix of the random vector xptq
• tr p¨q: the trace of a matrix
• |¨|: the cardinality of a set
• deg i: the degree of node i (i.e., |Ni |)

11

1.4.2

Acronyms

Here we define a number of acronyms that will be used throughout the work.
• AWGN: additive white Gaussian noise
• CLT: central limit theorem
• ECSQ: entropy-coded scalar (uniform) quantization/quantizer
• GP: geometric program/programming
• GGP: generalized geometric program/programming
• i.i.d.: independent and identically distributed
• LMMSE: linear minimum mean square error
• LSE: log-sum-exponential or log-sum-exp
• MMSE: minimum mean square error
• MSE: mean square error
• PMF: probability mass function
• pdf: probability density function
• RD: rate-distortion
• RV: random variable
• RVec: random vector
• SDR: signal-to-distortion ratio
• SNR: signal-to-noise ratio
• VQ: vector quantization/quantizer

12

• WSN: wireless sensor network

13

CHAPTER

2
QUANTIZED DISTRIBUTED AVERAGE CONSENSUS

This chapter presents much of the background required to understand consensus and source
coding. The results presented here are not original, and they are presented to make the thesis
self-contained. Our original results appear in Chapters 3 and 4.

2.1

Overview

One of the many attractive features of the variety of distributed average consensus algorithms
explored in this thesis is its linearity. At each iteration, every node takes a weighted average
of incoming messages from its neighbors, and in the absence of quantization errors, the
performance is elegantly described using concepts from spectral graph theory [6, 37].
Unfortunately, quantization introduces nonlinearity into the state update, which complicates analysis. Using the additive quantization noise model [52] and dithering [53], however, it
is possible to linearize the impact of quantization on the performance. The use of dithering
has been found in previous work [19] to greatly simplify the design of distributed averaging
algorithms, because the quantization error becomes uncorrelated with the quantizer input
signal.

14

Although dithering and high-resolution assumptions greatly simplify the analysis of
distributed average consensus, close attention must be paid to the design of the state update.
If quantization errors are introduced into the lossless algorithm, then the network state will
not converge to the true sample average [21, 32]. However, by making a simple modification
to the state update, Frasca et al. [14] guaranteed convergence to within one quantization bin.
In contrast, prior efforts [21, 22, 32] could only guarantee bounded asymptotic mean square
error (MSE) due to the drift from the sample average in the presence of quantization error.
In this chapter, we explore the fundamentals of distributed average consensus, algorithmic
modifications to account for quantization effects, rate-distortion (RD) theory, and the basics of
the additive noise model for quantization error.

2.2

Mathematical preliminaries

Before introducing the consensus problem, it is first necessary to introduce a number of
mathematical concepts related to linear algebra and graph theory.

2.2.1

Matrix and graph theory concepts

To represent the network of interest, it is necessary to model the nodes, which represent the
computing elements or agents. These can be wireless sensors [7, 28], servers [9], cameras [8], or
robots [6]. We assume that each node can only communicate with a subset of the other nodes
of the network, so it is also necessary to model the presence or absence of communication
links between them, which can be wireless channels or wired connections. These relationships
are modeled by a graph [37].
In this thesis, the communication links are bidirectional connections, so we model the
network as an undirected graph. This graph can be represented in a number of ways. One of

15

the simplest representations of the graph is a pair of sets,

G “ tV , E u,

(2.1)

where the graph G is comprised of a set of vertices (nodes) V and a set of edges E between
pairs of vertices [54, Sec. 1.2]. Because the communication links are bidirectional, each edge
ti, ju P E is represented as an unordered pair of vertices i and j [37].

2.3

The consensus problem

In the simplest case of the consensus problem, each node i P t1, . . . , mu has an initial scalar
quantity
zi p0q P R,

(2.2)

and the goal is to have all nodes of the network agree upon the sample mean of these
quantities by iteratively exchanging messages with their neighbors [37]. The quantities zi ptq
will be referred to as “states,” which in this thesis are assumed to be real-valued scalar random
variables (RVs) with known distribution. More formally, let the (discrete) iteration index be a
m are the initial values to be averaged
nonnegative integer,1 t P Zě0 . At t “ 0, the states tzi ptqui“1

by the consensus algorithm. For t ě 1, the state zi ptq represents the estimate of the sample
average
m
1 ÿ
s
z :“
zi p0q
m
i“1

at node i. The objective of consensus is for the state zi ptq to eventually equal the sample mean of
the initial states [37]. Mathematically, this is expressed as limtÑ8 zi ptq “ s
z, @i P t1, . . . , mu [37].
1 We denote the positive subset of a set S by S . The nonnegative subset is similarly denoted S . The integers
ą0
ě0
are denoted by Z, and the real numbers by R.

16

We leave the discussion of vector-valued states for the following chapter. In this thesis, we
restrict our attention to deterministic, synchronous-update consensus algorithms. We assume
the following: (i) the communication link topology of the network is fixed and does not change
with time, (ii) at each iteration, every node of the network exchanges messages with only its
neighbors, and (iii) the communication channels between nodes are noiseless.
Given the above assumptions on communication, one of the most popular algorithms for
consensus relies on linear updates [6, 37]. Each node is assigned an index i P t1, . . . , mu. Let
zi ptq denote the state of the ith node at iteration t. Each node updates its state by taking a
weighted sum of its own state with those of its neighbors [37],

zi pt ` 1q “ wii zi ptq `

ÿ

wij z j ptq,

(2.3)

jPNi

where wij ą 0 @i, j,

řm
k“1

wik “ 1, and Ni denotes the neighborhood of node i,

Ni :“ tj | ti, ju P E u.
The degree of node i is defined as

deg i :“ |Ni |,

where | ¨ | represents set cardinality. By this definition, deg i is the number of neighbors of
node i. The weights wij are designed such that [37]

lim zi ptq “ s
z.

tÑ8

If the state of each node of the network is collected in a vector

zptq :“ rz1 ptq, . . . , zm ptqsJ ,

17

(2.4)

and the averaging weights wij are collected in a matrix,

rWsij :“ wij ,

(2.5)

then the above update equation (2.3) can be written in matrix-vector form as [37]

zpt ` 1q “ Wzptq.

(2.6)

The design of the weight matrix W that yields the fastest convergence is a well studied
problem; the interested reader is referred to Xiao and Boyd [37]. To converge asymptotically,
that is,

lim zptq “

tÑ8

1 J
11 zp0q “ s
z 1,
m

the weight matrix must be doubly stochastic,

ř
i

wij “

ř
j

wij “ 1, and the modulus of its

largest eigenvalue must be less than unity [37].

2.3.1

Lossy consensus

By introducing quantization error into the internode messages, the simple linear iteration
above (2.6) is not guaranteed to converge. Instead, we use the modified iteration proposed
by Frasca et al. [14], which allows the sample average to be preserved in the presence of
quantization error.
Let Q : Rm Ñ X m represent quantization to a finite set of representation levels X m Ă Rm
(i.e., Qpzptqq “ rQ1 pz1 ptqq, . . . , Qm pzm ptqqsJ ). The subscripts on Q indicate that each node can
use a different quantizer in general. State updates for the case of consensus with quantized
messages have been studied by a number of authors, as discussed in Chapter 1. In this thesis,

18

we use the update proposed by Frasca et al. [14], which is

zi pt ` 1q “ zi ptq `

m
ÿ

`
˘
wij Q j pz j ptqq ´ Qi pzi ptqq ,

(2.7)

j“1

or in matrix-vector form,

zpt ` 1q “ zptq ` pW ´ Iq Qpzptqq.

The key advantage of this update is that the average

1
m

řm

i“1 zi ptq

(2.8)

of the states zi ptq is preserved

at each step t, despite the presence of quantization error [14]. Defining the quantization error

eptq :“ Qpzptqq ´ zptq,

(2.9)

the preservation of the sample mean can be seen by rewriting (2.8) as

zpt ` 1q “ zptq ` pW ´ Iq pzptq ` eptqq
“ Wzptq ` pW ´ Iq eptq,

(2.10)

and then taking the sample mean [14],
1 J
1
1
1 zpt ` 1q “ 1J Wzptq ` 1J pW ´ Iq eptq.
m
m
m
Because W is doubly stochastic, 1J W “ 1J , and the mean at iteration t ` 1 is equal to the

19

mean at iteration t [14]:
m
1
1 ÿ
zi pt ` 1q “ 1J zpt ` 1q
m
m
i“1

1 J
1
1 Wzptq ` 1J pW ´ Iq eptq
m
m
¯
1
1 ´ J
“ 1J zptq `
1 ´ 1J eptq
m
m looooomooooon

“

“0
m
1 ÿ
1
zi ptq.
“ 1J zptq “
m
m

(2.11)

i“1

In this thesis, we model the quantization error eptq as additive noise.
Note that the target state of consensus, termed the average consensus state, can be written,

z˚ :“ s
z1,

which means that zi ptq “ s
z, @i P t1, . . . , mu. We also define the average consensus operator

1
J
m 11 ,

so that [37]
1 J
11 zp0q “ z˚ .
m
Defining the error from the true sample mean,

eptq :“ zptq ´ z˚
“ zptq ´

1 J
11 zp0q,
m

and noting that the average is preserved over the iterations, (i.e.,
error eptq can be expressed as [14]
ˆ
eptq “

˙
1 J
I ´ 11 zptq.
m

20

(2.12)
1 J
m 1 zptq

“

1 J
m 1 zpt ` 1q),

the

2.4

Source coding, quantization, and rate-distortion theory

Digital systems rely on the ability to transmit and store information. To perform these
tasks reliably, the information takes the form of strings of symbols belonging to some finite
alphabet [30]. As a cornerstone of digital communication and storage, coding has been studied
quite extensively. In this section, we present some of the fundamental coding and quantization
concepts required to understand this thesis.

2.4.1

Coding and information theory

Understanding data compression requires a few concepts from information theory [30], initially
developed by Shannon in the 1940s [55]. This theory offers insight into the fundamental limiting
performance of communication and compression systems [30]. In this subsection, we focus on
the case where the source to be encoded, x P Rn , is mapped to a point x̂ in some finite set X n
by the quantization operator Qp¨q.
The entropy of a scalar RV x P X with probability mass function (PMF) ppxq is given by [30]
ÿ
Hpxq :“ ´ ppxq log2 ppxq.
xPX

Intuitively, the entropy of an RV is a measure of its uncertainty or information content [30].
Let x P Rn represent a long sequence x, the entries of which are independent and identically
distributed (i.i.d.) according to the probability mass function (PMF) ppxq. Then nHpxq is the
minimum expected binary sequence length required to describe x without error [30]. That
is, if we wish to describe the sequence x by a string of binary digits x̃ P t0, 1u M , then a code
exists such that the original sequence x can be noiselessly reconstructed from x̃ provided that
E rMs ě nHpxq [30]. To extend this concept to sequences of arbitrary length, we define the
coding rate per symbol as n1 E rMs, which is the average number of bits used to describe a single
source symbol [30].
Another important quantity in information theory is the mutual information between two

21

RVs. Let x and y be RVs with joint PMF ppx, yq and marginal PMFs ppxq and ppyq, respectively.2
Then the mutual information between x P X and y P Y is given by [30]

Ipx; yq :“

ÿ ÿ

ppx, yq log2

xPX yPY

ppx, yq
.
ppxqppyq

The utility of these information-theoretic quantities will become apparent during the discussion
of RD theory.

2.4.2

Rate-distortion theory

The previous discussion of coding considered only RVs that take on a finite set of values. If
we wish to digitally communicate or store a continuous source, it must first be quantized [35,
Ch. 1]. The two most basic elements of a quantizer are a set of representation levels, which are
used to approximate the unquantized signal, and a set of decision thresholds, which determine
the mapping from the input set to the output set [35, pp. 133–135]. If we imagine the source
data as an RV x P Rn , then the quantizer Q : Rn Ñ X n maps x to one of finitely many
representations x̂ P X n [35]. For real-valued sources, quantizers necessarily introduce a certain
expected distortion D into their representation of the input signal [35, pp. 144–145]. This
distortion can be quantified using a number of metrics, but for the purpose of this thesis, we
use the square error metric [35]
dpx, x̂q “ kx ´ x̂k22 ,
“
‰
so that the expected distortion per entry or per letter is given by D “ n1 E kx ´ x̂k22 . In general,
using a higher coding rate R results in a lower distortion D, with the drawback of greater
communication load. RD theory [29, 30] quantifies the best possible trade-off between coding
rate and distortion using the tools of information theory.
2 This is an abuse of notation: it would be more correct to write p pxq and p pyq, indicating that the two
x
y
functions are, in general, different. Here the subscripts are omitted for simplicity.

22

Gaussian RD Function

7

6

R(D) (bits/symbol)

5

4

3

2

1

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

D/σ 2
Figure 2.1 The RD function RpDq of a Gaussian memoryless source. The function is a nonincreasing,
convex function of the normalized distortion σD2 , where σ2 is the variance of the Gaussian source. The
RD pairs above the curve are achievable, while those below the curve are not.

23

Assume that we wish to encode a memoryless source, which is represented by the random
vector (RVec) x P Rn , xi i.i.d. „ f pxi q, @i P t1, . . . , nu. For the remainder of this subsection,
we drop the indices on xi due to the assumption that all entries of x are i.i.d. The minimum
coding rate R required to produce an expected distortion less than or equal to a particular
value D is given by the so-called RD function RpDq [29, 30],

RpDq :“

min

ppx̂|xq:Erdpx,x̂qsďD

Ipx; x̂q.

In words, the RD function is the minimum of the mutual information over all possible “test
channels” ppx̂|xq, subject to the constraint that the expected distortion per entry E rdpx, x̂qs is
less than a specified value D [30]. Operationally, the RD function is the minimum number
of bits per symbol required to describe a long i.i.d. source within the prescribed distortion
D [30]. The computation of a closed form for RpDq is difficult in general; however, RpDq can be
computed numerically [57–59]. When a particular quantizer is used, it will often have an RD
trade-off curve that differs from RpDq, which is a bound on the best possible performance [35].
In this thesis, we term such a trade-off curve for a particular practical quantizer an operational
RD relationship to avoid ambiguity. An example RD function for the Gaussian case is given
in Figure 2.1. Figure 2.2 shows the performance of practical quantizers compared to both
the Gaussian RD function and their respective operational RD relationships. Note that the
signal-to-distortion ratio (SDR),

SDR :“

var pxq
,
D

increases approximately exponentially in the high-rate limit for all practical quantizers in Figure 2.2. This exponential growth is a common feature of Gaussian operational RD relationships
that will be exploited in the following chapter.

24

Asymptotic Gaussian RD Relationships

5

4.5

4

R (bits/symbol)

3.5

3

2.5

2

ECSQ asymptote
ECSQ
Lloyd asymptote
Lloyd
Uniform asymptote
Uniform
Gaussian RD function

1.5

1

0.5

0
-5

0

5

10

15

20

25

30

10 log10 σ 2 /D (dB)
Figure 2.2 Operational RD relationships of practical quantizers: ECSQ, Lloyd-Max, and fixed-rate
uniform quantization. Note that all of these approach an asymptote that is a constant offset from the
Gaussian RD function as R Ñ 8. The ECSQ performance corresponds to the uniform scalar quantizer
with a level at zero. The representation level placement may be suboptimal, but it allows for coding
rates below one bit per symbol. This figure was inspired by a similar plot in lecture notes from a TU
Berlin source coding course [56].

25

2.4.3

Uniform and nonuniform quantization

From a conceptual standpoint, one of the simplest quantizers is the familiar uniform scalar
quantizer, for which both the representation levels and the decision thresholds are uniformly
spaced [35]. For a given input x, a uniform scalar quantizer simply rounds each xi to the
representation level x̂i nearest to xi [35]. This type of quantizer is frequently encountered
in analog-to-digital conversion and various digital signal processing systems. Despite its
simplicity, this quantizer performs surprisingly well in certain scenarios.
Consider a source x P Rn that is not uniformly distributed. For a uniform scalar quantizer,
each representation level x̂i will occur with different probability. In this case, we can choose
the representation levels based on the input statistics to obtain a low-distortion quantizer [35].
This idea was explored by Lloyd [60] and Max [61], who developed necessary and sufficient
conditions for quantizer optimality, given a fixed number of representation levels and a certain
source distribution. Lloyd and Max also developed algorithms for achieving an efficient
quantizer by iteratively updating the thresholds and representation levels [35]. Because the
quantizer is optimized for a particular number of representation levels, Lloyd-Max quantizers
are useful for fixed-length codes, where the encoded sequence x̃ P t0, 1unR . However, more
sophisticated coding techniques, such as Huffman coding [62], can be used to approach the
previously discussed entropy lower bound on code length if the length M of the binary code
sequence x̃ P t0, 1u M is allowed to vary with the input [35].
In the case of variable-length codes, we encounter a surprising result. For certain source
distributions, uniform scalar quantization followed by entropy coding, called entropy-coded
scalar quantization (ECSQ), is the optimal scalar quantizer asymptotically as R Ñ 8 [35, 63].
In fact, for the memoryless Gaussian source, ECSQ has been proved to approach the RD
limit within 0.255 bits as R Ñ 8 [35, 64]. For many memoryless sources, the performance of
ECSQ is similar [35, 64]: Farvardin and Modestino [65] showed less than 0.3 bits redundancy
over RpDq in the worst case for all the memoryless distributions they studied. Similarly, the
operational RD relationships for Lloyd-Max quantization and fixed-rate uniform quantization

26

approach asymptotes that are a constant offset from the RpDq function in the limit R Ñ 8 [35].
The asymptotic offset is illustrated in Figure 2.2.
The curious reader might wonder how the optimal quantizer can differ from the bound
if the bound is, in fact, achievable. Better performance can be achieved by quantizing long
sequences of the source jointly. Interestingly, this is advantageous even if the source symbols
are independent [30]. This technique is termed vector quantization (VQ) [35].
Because of the performance improvement associated with VQ, we also consider the
multidimensional extension of the scalar quantizer, which is the lattice quantizer [35]. These
quantizers approach the RD lower bound for memoryless Gaussian sources [66] as the block
dimension approaches infinity, while also permitting the application of dithering, which we
will discuss in the following sections.

2.4.4

The quantization noise model and dithering

One of the nicest properties of distributed averaging algorithms is their linearity [14]. At each
iteration, a weighted sum of the incoming messages is computed at each node [6, 37]. This
allows us to exploit the properties of linear operators, and, for Gaussian-distributed initial
states, we can describe the entire network state statistics using linear algebra [21]. However,
quantization is inherently nonlinear, which complicates the analysis [14]. Luckily, under certain
conditions, the quantization error can be modeled as an additive noise that is uncorrelated
with the source [52, 53, 67].
In undithered quantization systems, the conditions for i.i.d. uniform quantization error
are rather restrictive [53]. In most cases, these conditions are not met exactly, and the additive
quantization noise model is used as an approximation to the true behavior that simplifies
system design [52]. In general, the approximation resulting from assuming the additive
quantization noise model improves with higher coding rates [52, Sec. 9.5].
More precisely, let x P Rn represent a vector to be quantized. The quantizer Qp¨q : R Ñ X n
maps the values of x to a set of representation levels x̂ P X n [35]. Defining the quantization

27

error,

e :“ x̂ ´ x,

the quantization noise model can, for our purposes, be summarized by the following three
relationships: (i) the quantization error is zero-mean,
E res “ 0,

(ii) the quantization error is uncorrelated with the source,
”
ı
E exJ “ 0,

and (iii) the quantization error entries resi and res j , i ‰ j, are uncorrelated,
”
ı
b2
I,
E eeJ “
12
where b is the quantizer bin size and I is the identity matrix [52].
In situations that discourage or forbid the use of high rates, a technique called dithering
can be used to cause the quantization error to be uncorrelated with the source [53]. Dithering
is the process of intentionally adding noise to the source signal x P Rn to randomize the
quantization error, which is otherwise a deterministic function of the quantizer input [53].
If the i.i.d. noise w P Rn satisfies certain technical conditions, and the receiver subtracts w
from the reconstructed signal Qpxq, then the previously discussed additive quantization noise
model becomes exact in the sense that [53]
”
ı
E pQpx ` wq ´ px ` wqq xJ “ 0

28

and
”
ı
b2
I.
E pQpx ` wq ´ px ` wqqpQpx ` wq ´ px ` wqqJ “
12
Therefore, the quantization error can be modeled as an additive white noise with variance
b2
12 ,

which is uncorrelated with the input signal x [53]. The subtraction of the dither can be

achieved in practice using a particular pseudorandom number generator with a seed that is
known to both the transmitter and receiver.
Perhaps the simplest dither signal that causes the quantization error to be independent of
the source and spectrally white is the uniform distribution U p´ 2b , 2b q [53]. Uniformly distributed
dither w „ U p´ 2b , 2b q is used throughout this thesis whenever dithering is mentioned.

29

CHAPTER

3
CODING RATE OPTIMIZATION IN CONSENSUS

Depending on the system of interest, communication may be more or less expensive relative
to run time [47]. For instance, consider two extremes of distributed computing: wireless sensor
networks (WSNs) [1, 2] and multiserver networks [3]. WSNs have limited battery life and must
communicate over wireless channels, which consumes a large amount of power. In many WSN
applications, communication is expensive due to the low-power constraints at each node [47].
By contrast, multiserver networks, such as those used in cloud services [3] have less stringent
communication requirements, and the run time is more significant to the user [47]. In these
networks, the nodes communicate over wired connections and do not have battery power
constraints. To produce an optimization that is useful in diverse distributed settings, then, we
pursue a general and simple approach that models the relative costs of communication and
computation.
The contents of this chapter constitute the main analytical results of this thesis. The
contributions of this chapter are twofold: (i) a set of state evolution equations for the network’s
sufficient statistics are derived under the assumption of initial Gaussian distribution and
additive white quantization error, and (ii) it is proven that the problem of minimizing the
aggregate rate required to achieve a certain mean square error (MSE) can be posed as a

30

generalized geometric program (GGP).
The first contribution is minor, since similar approaches to computing the network state
statistics were derived in Yildiz and Scaglione [21] and Huang and Hua [28]; however, our
derivation allows for the explicit computation of MSE expressions in closed form using the
iteration (2.7) introduced in Frasca et al. [14]. Yildiz and Scaglione [21] only derived a closed
form for the asymptotic MSE as the iteration index t Ñ 8, and Huang and Hua optimized
with respect to an MSE bound. Our analysis allows an exact optimization to be performed in
the finite-iteration (total iterations T ă 8) case.
The second contribution allows the problem to be formulated in a well-understood framework. The GGP model can be converted into convex form and efficiently solved using common
numerical techniques [31]. Furthermore, due to the convexity of the equivalent transformed
optimization problem, the solver is guaranteed to find the globally optimal solution [31]. In
addition to the GGP model, we present a heuristic approach for optimization of fixed-rate
source coding schemes. This model reduces the size of the optimization problem, which still
provides the same solutions as exhaustive search.

3.1

Problem formulation

At every iteration t P t0, . . . , T ´ 1u of the consensus process, each node i P t1, . . . , mu uses a
coding rate Ri ptq to encode its state for transmission to the neighboring nodes.1 In general,
Ri ptq can vary across both nodes and iterations. In practice, we want to terminate the consensus
process after a finite number of iterations T. We simplify notation by defining the coding rate
vector
R :“ rR1 p0q, . . . , Rm p0q, . . . , R1 pT ´ 1q, . . . Rm pT ´ 1qsJ .
1 The

choice of zero as the initial iteration is arbitrary.

31

(3.1)

We denote the expected distortions per entry incurred by transmitting states using the coding
rates R by the distortion vector
D :“ rD1 p0q, . . . , Dm p0q, . . . , D1 pT ´ 1q, . . . Dm pT ´ 1qsJ .

(3.2)

One key quantity we use to determine the cost of running the consensus process is the
aggregate coding rate [47, 51]:

Ragg :“

T´1
m
ÿÿ

Ri ptq,

(3.3)

t“0 i“1

which represents the total coding rate used over the T iterations of the consensus algorithm by
all m nodes of the network. To address the trade-off between communication and computation,
we borrow the following composite cost function from Zhu and coauthors [47–51] that models
both communication and run-time costs:

Ctotal pR, Tq :“ K1 Ragg ` K2 T,

(3.4)

where K1 and K2 represent the costs of communication and computation, respectively. Since
the iteration (2.8) consists only of a quantization and a weighted sum at each node, the total
run-time cost is modeled as linear in the number of iterations. This cost function can be used
to model both battery-powered sensor networks and server networks by varying the relative
values of K1 and K2 [47, 51].
In this chapter, we present minimization strategies of the subcost

Ccomm pR, Tq :“ K1 Ragg

(3.5)

for fixed- and variable-length codes for (i) Gaussian-distributed sources using a variety
of quantizers, and (ii) arbitrarily distributed sources using a specific fixed-rate uniform
quantizer. For both proposed schemes, the overall solution to the minimization of (3.4) can

32

be found by fixing T and minimizing the communication cost Ccomm pR, Tq, subject to the
constraint that R achieves a target mean square error (MSE) value. Repeating this procedure
for T P tTmin , . . . , Tmax u is equivalent to minimizing over T and R jointly.
To efficiently encode the data stored across the network, it is necessary to understand
the evolution of its distribution over time. In general, this is not tractable due to the need to
track dependences among source elements and transformations of source densities. Instead,
we propose an efficient optimization scheme for entropy-coded uniform scalar quantization
(ECSQ) [35] of stationary Gaussian states and rate-distortion-optimal (RD-optimal) vector
quantization (VQ) [29, 30] of memoryless Gaussian-distributed initial states.2 We also use the
resulting strategies to design a coding rate optimization heuristic for a particular dithered
fixed-rate uniform quantizer [52, 53], for which the expected distortion is the same regardless
of source distribution [53].
In this thesis, we assume that the node-to-node communication is taking place in a broadcast
manner [36], which means that the total rate expended by a transmitting node is independent
of the size of its neighborhood. This also simplifies the analysis of the introduced distortion,
because each neighbor of a particular node will receive the same quantized message.

3.2

Abstractions for vector-valued data

The concepts presented in this section are included for clarity. We will later present simplifications for the case that each node has a stationary Gaussian source. Although not used directly
in our optimization model, the mathematics presented in this section is useful for computing
the source statistics across the network in a general setting.
To model the vector-valued data stored at each node during the averaging process, a useful
abstraction is introduced, based on similar approaches by Huang and Hua [28]. Let each node
i P t1, . . . , mu at iteration t P t0, . . . , T ´ 1u of consensus have state ζ i ptq P R L . We represent the
2 For

the scalar quantization schemes, we assume that the state of each node is a sample from a stationary,
ergodic Gaussian random process. The scalar quantizer acts elementwise, so we expect the RD performance of
correlated sources with scalar quantization and entropy coding to be the same as the memoryless RD performance.

33

jth element of the state ζ i ptq as
j

ζ i ptq :“ rζ i ptqs j .

(3.6)

”
ı
These states are Gaussian random vectors (RVecs) that may be correlated, so that E ζ ik ptqζ lj ptq ‰
0 in general.
To model the state of the entire network, we define the network state supervector ζptq P RmL
as follows:
”
1
ζptq :“ ζ 11 ptq ¨ ¨ ¨ ζ m
ptq

ˇ
ˇ
ˇ
ıJ
ˇ 2
ˇ
ˇ L
L
2
ζ
ptq
¨
¨
¨
ζ
ptq
¨
¨
¨
ζ
ptq
¨
¨
¨
ζ
ptq
.
ˇ 1
ˇ
ˇ 1
m
m

(3.7)

During consensus, each partition of ζptq in (3.7) is averaged independently (i.e., the vector
states are averaged over the network nodes elementwise). To represent this elementwise action,
we use the direct sum operator for matrices [68]:
»

fi
A1

—
—
—
n
à
—
Ai :“ —
—
i“1
—
–

0

0
..
.

¨¨¨

A2 ¨ ¨ ¨
.. . .
.
.

0

0

0
0
..
.

ffi
ffi
ffi
ffi
ffi .
ffi
ffi
fl

¨ ¨ ¨ An

À
n . Note that
Therefore, the operator
creates a block-diagonal matrix from the matrices tAi ui“1
Àn
i“1 Ai is shorthand for A1 ‘ A2 ‘ ¨ ¨ ¨ ‘ An . This is an alternative to the Kronecker product
notation and the matrix structure introduced by Huang and Hua [28].
Defining the weighted averaging matrix Ω P RmLˆmL :
Ω :“

L
à

W

(3.8)

l“1

as the operator equivalent to the weight matrix W (2.5) for vector consensus, the state update

34

equivalent to (2.8) for vector consensus is

ζpt ` 1q “ ζptq ` pΩ ´ IqQpζptqq,

(3.9)

where
ˇ
ˇ
ˇ
”
ıJ
ˇ
ˇ
ˇ
1
2
L
Qpζptqq :“ Q1 pζ 11 ptqq ¨ ¨ ¨ Qm pζ m
ptqq ˇ Q1 pζ 12 ptqq ¨ ¨ ¨ Qm pζ m
ptqq ˇ ¨ ¨ ¨ ˇ Q1 pζ 1L ptqq ¨ ¨ ¨ Qm pζ m
ptqq .
As in the scalar consensus case, we define the quantization error,

δptq :“ Qpζptqq ´ ζptq.

(3.10)

This is the vector consensus analog to eptq (2.9). Using the preceding definitions, the state
update (3.9) can be rewritten as
ζpt ` 1q “ Ωζptq ` pΩ ´ Iqδptq.

(3.11)

The goal of vector consensus is to iteratively drive each ζ i ptq to the sample mean of the initial
states ζ i p0q. Define s
ζ P R L as the sample mean (over nodes) of the initial states,
m
1 ÿ
s
ζ :“
ζ i p0q.
m
i“1

The goal of consensus is thus

lim ζ i ptq “ s
ζ, @i P t1, . . . , mu.

tÑ8

We term the state supervector ζ ˚ P RmL for which ζ i ptq “ s
ζ, @i P t1, . . . , mu, the vector average

35

consensus state. This state has the following form:
ffJ

«
ζ˚ “

ˇ “ ‰
“ ‰
“ ‰ ˇˇ “ ‰
“ ‰ ˇˇ
“ ‰
ˇ
s
ζ 1¨¨¨ s
ζ 1 ˇ s
ζ 2¨¨¨ s
ζ 2 ˇ ¨¨¨ ˇ s
ζ L ¨¨¨ s
ζ L

.

The vector average consensus operator is defined as

M :“

L
à
1 J
11 ,
m
l“1

(3.12)

so that
Mζp0q “ ζ ˚ .

The state estimation error, which is the error in the current estimate of the true average, is
defined as

ηptq :“ ζptq ´ Mζp0q.

(3.13)

Because Mζpt ` 1q “ Mζptq “ Mζp0q [this follows from (2.11)], the error (3.13) can also be
expressed as

ηptq “ pI ´ Mq ζptq.

3.3

Assumptions and analytical results

In this section, we present certain assumptions and results on the distribution of ζptq (3.7),
the network graph structure (i.e., G “ tV , E u (2.1)), and the quantization error δptq (3.10).
j

We also present relationships among the MSE, the variance of ζ i ptq (3.6), and the distortions
introduced, D (3.2).

36

3.3.1

Additive noise model

Following Frasca et al. [14] and Yildiz and Scaglione [21, 22], we adopt the additive model for
quantization noise [52]. This model assumes that the errors introduced by the quantization
process are uncorrelated with the quantizer input. A formal statement of the corresponding
assumptions is provided.
We assume that the quantization error δptq is spectrally white:
Assumption 1 The quantization error element δik ptq corresponding to ζ ik ptq (3.6) (i.e., quantization
error at the kth element of the ith node at iteration t) is uncorrelated with all other quantization error
”
ı
ı
“
‰ ”
elements: if j ‰ i, l ‰ k, or τ ‰ t, E δik ptqδjl pτq “ E δik ptq E δjl pτq “ 0.
We further assume that the quantization error is independent of the source to be quantized:
Assumption 2 The quantization error δptq is orthogonal to the state supervector ζptq with respect to
”
ı
node index (i.e., E ζptqδJ ptq “ 0).
These assumptions are well motivated for uniform quantizers with high coding rate [52], or
for uniform quantizers using dither at arbitrary rates [53]. In general, probability-densityoptimized (pdf-optimized) quantizers will have quantization error that is linearly correlated
with the input signal (i.e., E rxpQpxq ´ xqs “ ´D, where x is a scalar RV to be quantized,
Qpxq is the result of quantization, and D is the expected distortion) [35, p. 182, 69]. Although
dither for nonuniform quantization is not as well studied as for uniform quantization, the
problem of correlated quantization noise may be alleviated by dithering using a compander
as shown in a recent paper by Aykol and Rose [70]. For quantizers with nonuniform level
spacing, Assumptions 1 and 2 are not rigorously justified, but numerical results are included
to assess performance.

3.3.2

Preservation of Gaussian distribution

In this work, we assume that the initial state supervector ζptq (3.7) has a multidimensional
Gaussian distribution. Under certain conditions on the connectivity of the network, the

37

distribution of the state supervector ζptq will be Gaussian for all t P t0, . . . , T ´ 1u.
Lemma 1 Assume the technical conditions of Assumption 1 hold. Further assume that the degree of
each node is bounded from below by a monotone increasing function gpmq of the network size m. Then
the state supervector ζptq, t P t0, . . . , T ´ 1u, will be Gaussian distributed.
Proof: The proof follows from a simple application of the central limit theorem (CLT). At the
beginning of iteration zero, the random vector ζp0q is jointly Gaussian by assumption. The
state supervector is updated as (3.11),
ζpt ` 1q “ Ωζptq ` pΩ ´ Iq δptq.
If the elements of ζptq are jointly Gaussian, then the term Ωζptq will also be Gaussian
distributed, because it is a linear transformation of jointly Gaussian RVs. As m Ñ 8,

|Ni | ě gpmq Ñ 8, @i. Define the “aggregate quantization error” as
δagg ptq :“ pΩ ´ Iqδptq.
“
‰
Because the degree of each node deg i Ñ 8, @i P t1, . . . , mu, each element δagg ptq k of δagg ptq
is a sum of many independent RVs. By the CLT and our assumption of zero-mean independent quantization error, δagg ptq converges in distribution to N p0, pΩ ´ IqΣδ ptqpΩ ´ Iqq,
”
ı
where Σδ ptq :“ E δptqδJ ptq is the covariance matrix of the quantization error δptq. Because
both terms of (3.11) are Gaussian distributed, if ζptq is Gaussian, then ζpt ` 1q will be Gaussian. Therefore, by our assumption of initial Gaussian state, the state ζptq will be Gaussian

˝

distributed for all t P t0, . . . , T ´ 1u.

Note: For high-dimensional VQ, the quantization noise approaches a white Gaussian
process [66], and the above node degree bounding assumptions are unnecessary.

38

3.3.3

State evolution for vector consensus

Because the network state supervector ζptq (3.7) is Gaussian for all t, the sufficient statistics of
ζptq are a mean vector and covariance matrix. Likewise, because the error ηptq (3.13) is a linear
transformation of ζptq, it has a multidimensional Gaussian distribution. By deriving formulas
for the evolution of these statistics, we can minimize the cost function (3.5) to find the optimal
source coding rates Ropt .
We define the mean of ζptq (3.7) as
µζ ptq :“ E rζptqs,

(3.14)

µη ptq :“ E rηptqs.

(3.15)

and the mean of ηptq (3.13) as

The covariance matrices of ζptq and ηptq are defined as
”
ı
Σζ ptq :“ E ζptqζ J ptq ´ µζ ptqµJ
ζ ptq,

(3.16)

”
ı
Ση ptq :“ E ηptqηJ ptq ´ µη ptqµJ
η ptq,

(3.17)

and

respectively.
Lemma 2 The quantities (3.14)–(3.17) can be computed as follows,
µζ ptq “ Ωt µζ p0q,

39

(3.18)

µη ptq “ pI ´ Mq µζ p0q,

Σζ ptq “ Ωt Σζ p0qΩt `

t´1
ÿ

Ωt´s´1 pΩ ´ IqΣδ psqpΩ ´ IqΩt´s´1 ,

(3.19)

(3.20)

s“0

and
Ση ptq “ pI ´ Mq Σζ ptq pI ´ Mq ,

(3.21)

”
ı
where Σδ ptq :“ E δptqδJ ptq is the quantization error (3.10) covariance matrix.
Proof: The state update used in vector consensus is (3.9)

ζpt ` 1q “ ζptq ` pΩ ´ IqQpζptqq,

where ζptq is the state supervector (3.7). The mean µζ ptq can be written in terms of the weight
matrix Ω (3.8) and the quantization error δptq (3.10) as
µζ ptq “ E rΩζpt ´ 1qs ` E rpΩ ´ Iqδpt ´ 1qs.
Because Ω and pΩ ´ Iq are constants, this expression is equivalent to
µζ ptq “ ΩE rζpt ´ 1qs ` pΩ ´ IqE rδpt ´ 1qs.
Recalling that the mean of the quantization error is assumed to be zero,
µζ ptq “ Ωµζ pt ´ 1q.

40

By recursion on the above equation,
µζ ptq “ Ωt µζ p0q,
which is (3.18).
The covariance Σζ ptq can be written
ı
”
Σζ ptq “E pΩζpt ´ 1q ` pΩ ´ Iqδpt ´ 1qq pΩζpt ´ 1q ` pΩ ´ Iqδpt ´ 1qqJ
ı
”
´ E rΩζpt ´ 1qsE pΩζpt ´ 1qqJ .
Noting that Ω is a constant, and that Ω “ ΩJ , expanding the above expression gives3
”
ı
”
ı
Σζ ptq “E Ωζpt ´ 1qζ J pt ´ 1qΩ ` E pΩ ´ Iqδpt ´ 1qζ J pt ´ 1qΩ
ı
”
ı
”
` E Ωζpt ´ 1qδJ pt ´ 1qpΩ ´ Iq ` E pΩ ´ Iqδpt ´ 1qδJ pt ´ 1qpΩ ´ Iq
”
ı
´ ΩE rζpt ´ 1qsE ζ J pt ´ 1q Ω.

By Assumption 2, the second and third terms above are zero, so that
”
ı
Σζ ptq “ ΩE ζpt ´ 1qζ J pt ´ 1q Ω ` pΩ ´ IqΣδ pt ´ 1qpΩ ´ Iq ´ Ωµζ pt ´ 1qµJ
ζ pt ´ 1qΩ.
Grouping terms,
´ ”
ı
¯
Σζ ptq “ Ω E ζpt ´ 1qζ J pt ´ 1q ´ µζ pt ´ 1qµJ
pt
´
1q
Ω ` pΩ ´ IqΣδ pt ´ 1qpΩ ´ Iq.
ζ
Therefore,
Σζ ptq “ ΩΣζ pt ´ 1qΩ ` pΩ ´ IqΣδ pt ´ 1qpΩ ´ Iq.
3 The weight matrix W is symmetric for the optimal weights found by Xiao and coauthors [32, 37] and other
simple weighting strategies, such as Metropolis-Hastings [32] and max-degree [37].

41

Performing recursion on the above equation gives

Σζ ptq “ Ω Σζ p0qΩ `
t

t

t´1
ÿ

Ωt´s´1 pΩ ´ IqΣδ psqpΩ ´ IqΩt´s´1 ,

s“0

which is (3.20).
Because the error from the true sample mean is ηptq “ pI ´ Mq ζptq (from (3.12) and (3.13))
and pI ´ Mq “ pI ´ MqJ , it follows that

µη ptq “ pI ´ Mq µζ ptq,
which is (3.19), and
Ση ptq “ pI ´ Mq Σζ ptq pI ´ Mq ,

which is (3.21). Similar steps are taken to obtain the simplified forms of these expressions for

˝

scalar consensus.

The equations (3.18)–(3.21) can be simplified for optimization modeling in the case that
the states ζ i p0q are stationary Gaussian processes. In this scenario, the distortion incurred
by uniform scalar quantization at node i and iteration t can be modeled as the distortion
from quantization of an independent and identically distributed (i.i.d.) source under certain
conditions [52, Sec. 9.5]. In particular, the RD performance can be described well by the
marginal variance of the state ζ i ptq,
νi pD, tq “ E

„´

„´
¯2 
”
ı2
¯2 
”
ı2
ζ i1 ptq
´ E ζ i1 ptq “ ¨ ¨ ¨ “ E ζ iL ptq
´ E ζ iL ptq ,

“
‰
provided that the correlations E ζ ik ptqζ il ptq , k ‰ l, are not extreme and the quantizer bin size b
is not too coarse [52, Sec. 9.5]. For this scenario, we can consider the scalar states zi ptq P R (2.2),
i P t1, . . . , mu and the state vector zptq P Rm (2.4) in place of ζ i ptq P R L , i P t1, . . . , mu and
ζptq P RmL , respectively, where zi ptq is distributed according to the marginal pdf f pζ i1 ptqq,

42

@i P t1, . . . , mu.
We define the means of zptq (2.4) and eptq (2.12),
µz ptq :“ E rzptqs
and
µe ptq :“ E reptqs,
and the covariances of zptq and eptq,
”
ı
Σz ptq :“ E zptqzJ ptq ´ µz ptqµJ
z ptq

(3.22)

and
”
ı
Σe ptq :“ E eptqeJ ptq ´ µe ptqµJ
e ptq.
The equations (3.18)–(3.21) reduce to
µz ptq “ Wt µz p0q,

ˆ
µe ptq “

Σz ptq “ Wt Σz p0qWt `

t´1
ÿ

(3.23)

˙
1 J
I ´ 11 Wt µe p0q,
m

(3.24)

Wt´s´1 pW ´ IqΣe psqpW ´ IqWt´s´1 ,

(3.25)

s“0

43

and
ˆ
Σe ptq “

˙
ˆ
˙
1 J
1 J
I ´ 11 Σz ptq I ´ 11 .
m
m

(3.26)

Given the above state-evolution equations (3.23)–(3.26), we can determine the variance νi pD, tq
and MSE, MSEi pD, tq, at each node i P t1, . . . , mu. The source variance at node i and iteration t
is defined as
“
‰
νi pD, tq :“ E z2i ptq ´ E rzi ptqs2 .

(3.27)

Let us discuss our state estimation error (2.12) metrics. We define the node MSE, which
represents the MSE at node i and iteration t, as
“
‰
MSEi pD, tq :“ E ei2 ptq ,

(3.28)

where ei ptq is the ith element of eptq (2.12). We now turn to quantifying the estimation error
over the entire network. For this purpose, we define the network MSE, which is the arithmetic
m ,
mean of the node MSEs, tMSEi pD, tqui“1

MSEpD, tq :“

m
1 ÿ
MSEi pD, tq.
m

(3.29)

i“1

Using these definitions, we present the following mathematical relationships, which we term
the state evolution equations. These equations allow us to perform the optimization of the coding
rate vector R using the cost function (3.5).
Lemma 3 Given the above equations (3.25) and (3.26) for the state and error covariances, the variance
at node i at iteration t (3.27) is given by

νi pD, tq “ rΣz ptqsii ,

44

(3.30)

the node MSE (3.28) is given by
”
ı
MSEi pD, tq “ Σe ptq ` µe ptqµJ
ptq
,
e
ii

(3.31)

and the network MSE (3.29) is given by

MSEpD, tq “

¯
1 ´
tr Σe ptq ` µe ptqµJ
ptq
,
e
m

(3.32)

where tr p¨q denotes the trace of a matrix.
The derivations of equations (3.23)–(3.26) take almost identical steps to those in the proof of
Lemma 2. An example of the MSEi ptq predicted by (3.31) is given in Figure 3.1.

3.4

Optimization via generalized geometric programming

The key insight of this work is the ability to pose the optimization of the coding rate vector R as
a generalized geometric program (GGP), whose global optimum can be found efficiently [71].
Using variable-length codes [62], the coding rates
Ri ptq P Rą0 , @i P t1, . . . , mu, @t P t0, . . . , T ´ 1u,
in contrast to the fixed-length coding case where Ri ptq P Zą0 , @i P t1, . . . , mu, @t P t0, . . . , T ´ 1u.
For the case of variable-length codes, we propose a high-rate approximation to the operational RD relationship to find an efficient rate vector R capable of achieving a target value
of MSEpD, Tq (3.32).4 This modeling strategy will later be used to design an optimization
mT .
procedure for fixed-length coding, where R P Zą0
4 Note that the final network MSE corresponds to iteration index T, and not T ´ 1. This is because MSEpD, tq is
the network MSE after the end of t iterations. The quantity MSEpD, t ´ 1q would thus correspond to the end of the
pt ´ 1qth iteration, or the beginning of the tth iteration.

45

0.25

MSEi (t)

0.2

0.15

0.1

0.05

0
1
2
3

0
5

4

Iteration index t

10
5

15
20

Node index i

Figure 3.1 Surface showing the evolution of MSEi ptq for the case of lossless consensus. Note that the
value of MSEi ptq is not equal across all nodes for a particular iteration, and the evolution of MSEi ptq is
heavily dependent on network topology. Also observe that the MSE decreases monotonically at each
node as t Ñ 8. The variance νi ptq evolves in a similar fashion.

46

3.4.1

Cost function formulation

For the quantizers of interest operating on Gaussian sources, the operational RD relationship
(the reader is reminded that the this term was defined in Section 2.4.2) in the high-rate regime
has the form
$
ˆ 2˙
σ
1
’
’
& log2
` Rc ,
D
RpDq « 2
’
’
%
0,

D P p0, σ2 Dmax s
,

(3.33)

otherwise

where σ2 represents the variance of the data to be encoded and Dmax is some constant [63]. In
some cases, such as infinite-dimensional VQ with Gaussian memoryless sources and dithered
scalar uniform quantization [66], the relationship (3.33) holds for all rates. Note that (3.33) is
equal to the Gaussian RD function
$
ˆ 2˙
1
σ
’
’
& log2
,
D
RG pDq “ 2
’
’
%0,

D P p0, σ2 s
,
otherwise

with a modified domain and offset. In the following discussion, it is easier to consider the
inverse to the operational RD relationship, denoted DpRq.
For uniform scalar quantization of a memoryless Gaussian source followed by block
entropy coding, the true operational RD relationship is approximated well in the high-rate
regime by the Gaussian RD function plus 0.255 bits [35, p. 302, 64]. More generally, the
operational distortion-rate relationship for ECSQ at high rates is determined by

DpRq «

1 2hpxq ´2R
2
2
,
12

(3.34)

where hpxq is the differential entropy of the source x [64]. Similar results hold for Lloyd-Max
quantization [64, 72] and fixed-rate uniform quantization [64]. For Lloyd-Max quantization of
memoryless Gaussian sources at high resolution, the operational RD relationship is approxi-

47

mated by
1
DpRq «
12

3

„ż 8
f

1{3

pxqdx

2´2R ,

(3.35)

´8

where f pxq is the probability distribution of the source x [64, 72]. For Gaussian densities,
both (3.34) and (3.35) scale with the source variance σ2 . This scaling can be proven by
straightforward computation.
As in Section 2.4.2, we denote the representation levels of the N-codeword quantizer as
N . The maximum representation level of the quantizer is denoted x̂
tx̂i ui“1
max , and the minimum

representation level is denoted x̂min . For subtractively dithered fixed-rate uniform quantization
of any source at all rates with ∆ :“ σ1 px̂max ´ x̂min q (neglecting overload distortion) [64],
DpRq “

∆2 σ2 ´2R
2
.
12

In our case, the source variance is νi pD, tq, defined in (3.30), which is a function of the
initial state vector covariance Σz p0q (3.22) and the distortion vector D (3.2). It evolves as
described by (3.30). The operational RD relationship at all nodes i P t1, . . . , mu and iterations
t P t0, . . . , T ´ 1u can be expressed as follows. Any Gaussian scalar quantizer will have a
maximum possible distortion for a particular source. Because the distortion is proportional
to the source variance νi pD, tq in the Gaussian case, we consider the maximum distortion
when νi pD, tq “ 1. Letting Dmax denote this maximum possible distortion for a unit-variance
Gaussian source,
$
ˆ
˙
νi pD, tq
1
’
’
& log2
` Rc ,
Di ptq
Ri pD, tq “ 2
’
’
%0,

48

Di ptq P p0, νi pD, tqDmax s
.
Di ptq ą νi pD, tqDmax

This operational RD relationship is equivalent to
"
ˆ
˙ *
νi pD, tq
1
, c ` Rc
Ri pD, tq “ max log2
2
Di ptq
ˆ
"
*˙
1
νi pD, tq
“ log2 max
,k
` Rc , Di ptq ą 0,
2
Di ptq

(3.36)

where c and k are terms (constant in D and t) that depend on Dmax . Concretely, k “ 2´2Rc ą 0.
Given a number of iterations T, the goal is to minimize the aggregate coding rate (3.3),
subject to a constraint on the final MSE, MSEpD, Tq ď MSE˚ (3.32). For a particular number
of iterations T and target MSE˚ , the minimum number of iterations required to achieve that
MSE, Tmin “ arg minT tT | MSElossless pD, Tq ă MSE˚ u, can be readily obtained using the stateevolution equations (3.23)–(3.26). More formally, using the operational RD relationship (3.36),
the optimization problem is

minimize
D

T´1
m
ÿÿ
t“0 i“1

ˆ
"
*˙
νi pD, tq
1
log2 max
,k
` Rc ,
2
Di ptq

subject to MSEpD, Tq ď MSE˚ ,
Di ptq ą 0, @i, t.
Note that the above optimization can be rewritten as
˜
minimize
D

ln

T´1
m
źź
t“0 i“1

"

νi pD, tq
max
,k
Di ptq

subject to MSEpD, Tq ď MSE˚ ,

*¸
,
(3.37)

Di ptq ą 0, @i, t,
where the constants

1
2 lnp2q

and Rc are omitted, because they make no difference to the optimum

point of the minimization. We will now introduce GGP and show that the optimization (3.37)
reduces to such a problem.

49

3.4.2

Basics of generalized geometric programming

The following information can be found in the book Convex Optimization by Boyd and Vandenberghe [31]. For this tutorial, we stay close to the authors’ original notation. In the language of
geometric programming (GP), a function of the form
f pxq “ cx1a1 x2a2 ¨ ¨ ¨ xnan , c ą 0, xi ą 0 @i, ai P R @i,
is called a monomial [31, 71]. Similarly, a function of the form

f pxq “

k
ÿ

gi px1 , . . . , xn q,

i“1

where the gi px1 , . . . , xn q are monomials, is termed a posynomial (not to be confused with the
word polynomial from algebra) [71]. The word posynomial is intended to convey that the
coefficients and values of the variables are constrained to be positive [71]. Boyd and Vandenberghe [31] state a certain number of closure properties:

Posynomials are closed under addition, multiplication, and nonnegative scaling.
Monomials are closed under multiplication and division. If a posynomial is multiplied by a monomial, the result is a posynomial; similarly, a posynomial can be
divided by a monomial, with the result a posynomial [31, p. 161].

By introducing dummy variables and additional constraints, more general cost and constraint functions, called generalized posynomials, can be converted into posynomials and
used to formulate GPs [71]. A generalized posynomial is a function formed from posynomials
using addition, multiplication, the maximum operation, and positive exponentiation [71].
Generalized posynomials are closed under addition, multiplication, division by monomials,
and the maximum operation [31, 71]. A standard inequality-constrained generalized GGP is of

50

the form
minimize Cpx1 , . . . , xn q,
x1 ,...,xn

subject to

f i px1 , . . . , xn q ď 1, @i P t1, . . . , n f u,
gi px1 , . . . , xn q “ 1, @i P t1, . . . , n g u,

where the cost Cpx1 , . . . , xn q and all the inequality constraints f i px1 , . . . , xn q are generalized
posynomials, and all the equality constraints gi px1 , . . . , xn q are monomials [31]. By a logarithmic
change of variables and objective function transformation, GGPs can be cast in convex form
and efficiently solved using numerical techniques [71]. Note that there is an implicit constraint
that xi ą 0 @i P t1, . . . , nu [71].
We now turn to the problem of posing the optimization (3.37) as a GGP. Boyd and
Vandenberghe [31] state that applying a monotonically increasing function to the objective
results in an equivalent problem. Therefore, by applying the exponential function to the cost
function (3.37), we obtain the equivalent problem,

minimize
D

T´1
m
źź
t“0 i“1

"

*
νi pD, tq
max
,k ,
Di ptq

subject to MSEpD, Tq ď MSE˚ ,

(3.38)

Di ptq ą 0, @i, t.

3.4.3

Proof that the optimization is a GGP

Until now, the analysis has resembled the prior art. The state evolution equations (3.23)–(3.26)
were derived for a different update equation than ours by Yildiz and Scaglione [21]. Huang
and Hua [28] extended analysis of consensus to vectors, although they also used a different
state update. However, neither of these works optimize with respect to exact MSE values (3.28),
and Huang and Hua did not account for the possibilities of different coding rates at each node
and iteration or entropy coding.
The main original contribution of this work is to pose the coding rate optimization problem

51

with node- and iteration-dependent rates as a GGP, which can solve quickly and reliably for
the global optimum coding rate vector Ropt . The following theorem proves that the optimization (3.38) can be cast as such a problem.

Theorem The optimization problem (3.38) is a GGP.
Proof: To show that (3.38) is a GGP, our argument takes the following steps: (i) show that
i pD,tq
is a posynomial, (iii) show that maxt νD
, ku is a
i ptq
śT´1 śm
νi pD,tq
generalized posynomial, (iv) show that t“0
i“1 maxt Di ptq , ku is a generalized posynomial,

νi pD, tq is a posynomial, (ii) show that

and (v) show that

1
MSEpD, tq
MSE˚

νi pD,tq
Di ptq

is a posynomial.

(i) To show that the variance νi pD, tq is a posynomial, it is necessary to express νi pD, tq in
terms of the individual distortions Di ptq. From the state-evolution equation (3.30),
«
νi pD, tq “ W Σz p0qW `
t

t

t´1
ÿ

ff
W

t´s´1

pW ´ IqΣe psqpW ´ IqW

t´s´1

s“0

.
ii

Recall that Σe ptq is a diagonal matrix by Assumption 1 and (2.9),
$
’
’
& Di ptq, i “ j
rΣe ptqsij “

.

’
’
%0, otherwise
x sq “ Wt´s´1 pW ´ Iq, and denote the ith column of Wpt,
x sq as w
p i pt, sq P Rm . Each term
Let Wpt,
x sqΣe psqW
xJ pt, sq in the sum is then equal to
Wpt,
x sqΣe psqW
xJ pt, sq “
Wpt,

m
ÿ

p k pt, sqw
pJ
Dk psqw
k pt, sq.

k“1

Therefore,
t´1
m
”
ı
ÿÿ
“
‰
p k pt, sqw
pJ
νi pD, tq “ Wt Σz p0qWt ii `
Dk psq w
pt,
sq
.
k
s“0 k“1

52

ii

Because the diagonals of the outer product of a real vector with itself are always nonnegative,
“
‰
each coefficient in the above summation is nonnegative. The first term, Wt Σz p0qWt ii is
“
‰
p k pt, sqw
pJ
a monomial. It is clear that even if one or more of the w
k pt, sq ii terms from the
summation is zero, νi pD, tq can still be considered a posynomial by omitting those terms from
the summation. Therefore, νi pD, tq is a posynomial.
(ii) Di ptq is a monomial, and a posynomial divided by a monomial is a posynomial.
Therefore,

νi pD,tq
Di ptq

is a posynomial.

(iii) The maximum of a posynomial and a monomial is a generalized posynomial. Because
i pD,tq
k is positive, it is trivially a monomial. Therefore, maxt νD
, ku is a generalized posynomial.
i ptq
ś
i pD,tq
(iv) Generalized posynomials are closed under multiplication. Therefore, i,t maxt νD
, ku
i ptq

is a generalized posynomial.
`
˘
| sq “ I ´ 1 11J Wt´s´1 pW ´ Iq and let w
q i pt, sq P Rm denote the ith column
(v) Define Wpt,
m
| sq. The MSE at each node (3.28) is given by
of Wpt,
«ˆ

˙
ˆ
˙
1 J
1 J
t
t
MSEi pD, tq “
I ´ 11 W Σe p0qW I ´ 11
m
m
ff
t´1
ÿ
| sqΣe psqW
|J pt, sq ` µ ptqµJ ptq .
`
Wpt,
e

e

s“0

ii

As in the proof of (i), the above equation reduces to
«ˆ
MSEi pD, tq “

˙
ˆ
˙
1 J
1 J
t
t
I ´ 11 W Σe p0qW I ´ 11
m
m

t´1
m
ÿÿ

`

ff
J
q k pt, sqw
qJ
Dk psqw
k pt, sq ` µe ptqµe ptq

s“0 k“1

.
ii

By the same argument as (i), MSEi pD, tq is a posynomial. Likewise, the network MSE,
řm
MSEpD, tq “ m1 i“1
MSEi pD, tq is a posynomial.
Because the objective and constraint functions are posynomials, and the variables Di ptq are
constrained to be positive, the optimization (3.38) is a GGP.

53

˝

Two aspects of the above proof should be highlighted. Because the constraints are allowed
to be generalized posynomials, one could also optimize with respect to a constraint on the
maximum node MSE, for example
m
max tMSEi pD, Tqui“1
ď MSE˚ ,
i

or constraints on each of the node MSE values,

MSEi pD, Tq ď MSE˚i , @i P t1, . . . , mu,
where MSEi pD, tq is defined in (3.31). Also, in its most general form, the optimization allows
each node to use a different rate or distortion. In the interest of designing a distributed
protocol (or for computational efficiency), one may wish to constrain the rates or distortions at
each node to be the same. The constraint that all distortions be the same is a straightforward
modification of (3.38) and is also provably a GGP.

3.5

Efficient implementation of GGP optimization

Solving the exact optimization problem (3.38) naively requires explicit representation of all mT
distortions Di ptq and all the coefficients of the log-sum-exp (LSE) model5 required to compute
MSEi pD, tq and νi pD, tq from D. The result of this explicit representation is large time and
memory complexity. In this section, we explore a simplification of the optimization problem
to combat these issues, and we present a heuristic to apply the GGP optimization method to
fixed-rate coding schemes.
5 GGPs

are converted to convex LSE form for solution [73, 74].

54

3.5.1

Equal-distortion simplification and implementation

According to Boyd and coauthors [31, 71], GGPs can be converted into convex optimization
problems by a change of variables and a logarithmic transformation of the objective function.
For a detailed treatment in the case of GP, the reader is referred to the above references. GGPs
can be converted in a similar fashion; however, the conversion is somewhat more complicated.
In some software packages, such as CVX [73, 74], GGPs can be specified directly, and they are
automatically converted into corresponding convex programs [71]. This approach guarantees
that the global optimum can be found with little additional effort by the programmer.
In practice, for large networks, (i.e., m ą 20 and T ě 6) the memory and time requirements
of the optimization (3.38) seem to grow very quickly. Furthermore, using CVX [73, 74],
all compatible solvers tested displayed poor stability for large networks. Attempts to use
alternative modeling frameworks, such as GPkit [75], CVXPY [76], and YALMIP [77], failed
due to the high memory requirements of the model. Manually building a GGP model is
memory- and time-intensive; if explicit model representation can be avoided, however, it
may be possible to apply other convex optimization methods without these scaling issues.
In particular, avoiding the explicit representation of posynomial coefficients and exponents
would probably greatly improve performance for large problem sizes m and T. To provide a
program that is more easily solved in practice, we make two simplifications. First, we constrain
the distortions to be equal at each node, which is well motivated for the end goal of designing
a truly distributed protocol. Next, the program can be cast as

minimize
D

T´1
m
źź

"
max

t“0 i“1

*
νi ptq
,k ,
Dptq

subject to MSEpD, Tq ď MSE˚ ,

(3.39)

Dptq ą 0, @t.
In the following chapter, the results of solving this simplified problem (3.39) are compared

55

to the solutions of the exact program (3.38) and the prior art [20, 21]. Surprisingly, the
above approximation provides competitive results for random geometric networks [45], with
significant reduction in memory and run-time requirements.

3.5.2

Efficient search heuristic for fixed-rate quantizers

Using fixed-length coding with uniform scalar quantization, the rates to be used are integervalued [64], and arbitrary distortions are not achievable for a given source variance. To
account for this constraint, we suggest a search heuristic. Given a particular T and MSE˚ ,
mT (3.1). Clearly, if each node is allowed to
there exist a finite number of rate vectors R P Zą0

use different rates, then the size of the search space becomes too large to be practical, and we
constrain our attention to the case that rates are constant over nodes for a particular iteration:
Ri ptq “ R j ptq “ Rptq, @i, j P t1, . . . , mu, t P t0, . . . , T ´ 1u. Even in the node-constant rate case,
the number of possible sequences such that mini,t Ri ptq “ 1 and maxi,t Ri ptq “ Rmax grows
combinatorially in T.
To find an efficient (though perhaps suboptimal) sequence, we propose to use the output
of the GGP (3.38) as a starting point for the optimization. Because constraining the rates to be
equal at each node is not straightforward, we instead constrain the distortions to be equal at
each node as in (3.39). The optimal rates Ri ptq from the solution of (3.39) are then averaged
T´1
into a single sequence tRGGP ptqut“0
,

RGGP ptq :“

m
1 ÿ
Ri ptq.
m
i“1

A neighborhood is formed as follows. To the original sequence

RGGP :“ rRGGP p0q, . . . , RGGP pT ´ 1qsJ ,
we add ` P t´1, `1uT . Each possible ` is added, and the result is stored, until all 2T possible
transformations of RGGP have been computed. Because the minimum fixed rate is one bit per

56

symbol, all rates less than one bit per symbol from the previous step are re-assigned to one
bit per symbol. Next, we apply the operators A P tt¨u, r¨suT to the resulting sequences, where,
for example, rt¨u, r¨s, t¨us rx, y, zsJ :“ rtxu, rys, tzusJ . Finally, we eliminate duplicate sequences.
This results in a “trellis” containing the integers about each RGGP ptq. The maximum possible
total number of sequences in the search space is thus 4T . Brute-force search is conducted on
this subset of sequences to find the one that satisfies the MSE constraint with the minimum
řT´1
aggregate rate Ragg “ m t“0
Rptq. In the case that multiple sequences achieve the MSE
constraint with equal aggregate rate, the one that produces the lowest MSE is selected.

57

CHAPTER

4
NUMERICAL RESULTS

In this chapter, we present numerical results that provide further insight into the optimal
rate sequences and their relationship to the target mean square error (MSE) and the network
topology. We further compare the performance of the proposed generalized geometric program
(GGP) and heuristic optimizations to the prior art. To test the effectiveness of the proposed
approach, we used MATLAB [78] for simulation. We used the CVX toolbox [73, 74] to solve the
GGP presented in Chapter 3, which converts the GGP model to a convex problem automatically.
In the results presented, the networks were generated by random geometric graph models [45], and each node state ζ i p0q was initialized with the same independent and identically
distributed (i.i.d.) variance-σx2 zero-mean Gaussian vector x P R L corrupted by a different
variance-σn2 , zero-mean Gaussian noise ni P R L , i P t1, . . . , mu. More precisely,
ζ i p0q “ x ` ni , @i P t1, . . . , mu,

(4.1)

where x „ N p0, σx2 Iq and ni „ N p0, σn2 Iq, @i P t1, . . . , mu. The random geometric graphs were
generated on the unit torus (i.e., edge effects were neglected by “wrapping” the top and
bottom and left and right edges of r0, 1s2 ) [45, 79, 80].

58

Figure 4.1 Optimal rates and MSE sequences from the solution of (3.38) (T “ 5, ρc “ 0.35, σx2 “ 1,
σn2 “ 0.5, m “ 20). Top: Optimal rate sequences for the variable-distortion optimization problem.
The rates are plotted against iteration indices, and each line represents the rates used by a different
sensor. Bottom: MSE sequence corresponding to the above rate sequence. The MSE values are plotted
against node indices, and each line represents a different iteration. In this case, the MSE values
decrease monotonically for all nodes. The lines represent values predicted by the state evolution
equations (3.23)–(3.26), and the overlaid circles represent the simulated values.

59

In the case of very low rates using entropy-coded uniform scalar quantization (ECSQ),
nothing would be communicated (all elements transmitted/received will be zeros). To prevent
this behavior, the maximum normalized distortion Dmax allowed was set such that 1% of the
elements were nonzero.

4.1

Structure of solutions

To verify the state-evolution equations (3.23)–(3.26) and better understand the structure of the
solutions to the optimization problems (3.38) and (3.39), we present some anecdotal simulation
results. Each of these results is taken from a different single instantiation of the optimization
problem.
T´1
The optimal rate sequences tRi ptqut“0
, i P t1, . . . , mu for both the variable-distortion

(3.38) and constant-distortion (3.39) problems typically exhibit monotonically nondecreasing
structure, with an increasing rate of change toward the final iterations. In the constant´
¯
i pD,tq
i pD,tq
distortion case, the rates Ri ptq « 12 log2 νD
in the
` Rc are similar because the ratios νD
i ptq
i ptq
Gaussian operational rate-distortion (RD) relationship (3.33) are similar across the network.
Examples of the optimal rate sequences and corresponding MSE state evolution are provided
for both variants of the optimization problem (3.38) and (3.39) in Figures 4.1 and 4.2.
The reader may wonder why the rate sequences of Figures 4.1 and 4.2 look different. In
m,T´1
Figure 4.1, the distortions tDi ptqui“1,t“0
corresponding to the optimization of (3.38) are allowed

to vary across both the node index i P t1, . . . , mu and the iteration index t P t0, . . . , T ´ 1u. This
m,T´1
optimization is more flexible with respect to the choice of rates tRi ptqui“1,t“0
corresponding
m,T´1
to tDi ptqui“1,t“0
. In Figure 4.2, however, the distortion is constrained to be the same the
m,T´1
network for a particular iteration (i.e., tDi ptqui“1,t“0
varies only across t P t1, . . . , T ´ 1u). This

corresponds to the solution of (3.39).
The pattern of these rate sequences is intuitive and mirrors the results of Zhu and coauthors’ study of multiprocessor approximate message passing [47, 48, 51]. As the estimate of
the sample mean at each node increases in precision, higher-resolution messages must be

60

Figure 4.2 Optimal rates and MSE sequences from the solution of (3.39) (T “ 5, ρc “ 0.35, σx2 “ 1,
σn2 “ 0.5, m “ 20). Top: Optimal rate sequences for the constant-distortion optimization problem. The
rates are plotted against iteration indices, and each line represents the rates used by a different sensor.
Bottom: MSE sequence corresponding to the above rate sequence. The MSE values are plotted against
node indices, and each line represents a different iteration. In this case, the MSE values decrease
monotonically at all nodes. The lines represent values predicted by the state evolution equations (3.23)–
(3.26), and the overlaid circles represent the simulated values.

61

exchanged among processors to achieve increasing estimate quality. In the case of coding
without side information, this implies the use of larger coding rates in the later iterations of
the algorithm.
The results of the variable-distortion problem (3.38) also exhibit structure based on the
connectivity of each node. Figure 4.3 contains a scatter plot of the aggregate rate versus the
node degree, which is a proxy for the “connectedness” of each node of the network. There
is a negative correlation between the node degree and the aggregate rate used by each node.
Figure 4.3 also contains a surface demonstrating the changing shape of the rate sequence as
the node degree decreases. In general, the experimental results suggest that nodes with higher
degree tend to have flatter optimal rate sequences, whereas the nodes with lower degree ramp
up their coding rates more aggressively from beginning to end.
Intuitively, this pattern can be explained by considering the overall goal of consensus and
the impact of quantization on the information content [44] of the network’s sample mean
estimates. The more poorly connected nodes must use greater rates at each iteration to diffuse
their information throughout the network, due to topological bottleneck effects [41]. The
more highly connected nodes, however, can save coding rate by exploiting their topological
advantage, as the information content of their messages will be diffused more quickly.
In all cases simulated, the differences in aggregate rate between the exact optimization
problem (3.38) and the constant-distortion problem (3.39) were small, on the order of 5%. For
the sequences plotted in Figures 4.1 and 4.2, the difference was 4.45%. This similarity may
be an effect of the highly regular random geometric graph structure. For random geometric
networks on squares rather than tori, the difference may be more significant due to the
irregularity arising from edge effects.

4.1.1

Nonmonotone behavior

Intuitively, it seems reasonable to expect the MSE to decrease monotonically, and the rates to
increase monotonically, since greater precision is required to encode the increasingly precise

62

Figure 4.3 Influence of connectivity on the optimal rates from the solution of (3.39) (T “ 5, ρc “ 0.35,
σx2 “ 1, σn2 “ 0.5, m “ 20). Top: Scatter plot of aggregate rate used at each node versus node degree.
Bottom: Surface of optimal rate sequences at each node, where the node indices are sorted by node
degree in decreasing order.

63

Figure 4.4 Nonmonotone optimal rate sequences for the low-SNR, low-MSE setting (T “ 5, ρc “ 0.35,
σx2 “ 1, σn2 “ 20, m “ 20). Each line corresponds to the rate sequence used by a different node. Top:
Variable-distortion optimization (3.38) result. Bottom: Constant-distortion optimization (3.39) result.

64

estimates of the sample average available to each node. However, in certain edge cases, neither
of these patterns holds. In this subsection, we discuss these behaviors and present some
hypothesized reasons for these nonintuitive results.
Figure 4.4 demonstrates a case in which the optimal rate sequences are nonmonotone. In
this case, the signal-to-noise ratio (SNR),

SNR :“ 10 log10

σx2
,
σn2

is low (-13.01 dB), which introduces a trade-off for the rates Ri ptq corresponding to the first
few iterations of the consensus algorithm. Recall that the simulation initializes each node to
have access to the same zero-mean Gaussian signal plus its own realization of additive white
Gaussian noise (AWGN). In this case, the covariance matrix (3.25) has the form
Σz pt “ 0q “ σx2 11J ` σn2 I.

(4.2)

If the signal variance is set to zero (corresponding to the case of averaging independent
AWGN processes), then the diagonals of the covariance matrix (3.25) will converge to

σn2
m

as

the iteration index t Ñ 8, which is significantly less than the original variance σn2 for large
networks.
Therefore, for low SNR, the impact of quantization error on the variance becomes significant.
The variance state-evolution equation is (3.25)
«
νi ptq “ Wt Σz p0qWt `

t´1
ÿ

ff
Wt´s´1 pW ´ IqΣe psqpW ´ IqWt´s´1

s“0

.
ii

It can be seen from the above equation that both the original covariance matrix and the
distortions introduced in past iterations play a role in the variance at any iteration t P Zą0 .

65

Figure 4.5 Nonmonotone MSE sequences. Each line corresponds to the MSE sequence of a different
node (T “ 5, ρc “ 0.35, σx2 “ 1, σn2 “ 0.5). Top: MSE evolution for the variable-distortion problem (3.38).
Bottom: MSE evolution for the constant-distortion problem (3.39).

66

Furthermore, recall the family of operational RD relationships we study (3.33),
ˆ
"
*˙
νi pD, tq
1
,k
` Rc ,
Ri pD, tq « log2 max
2
Di ptq
which indicates that, holding the distortion Di ptq constant, larger coding rates are required if
the source has higher variance νi ptq. From (3.30) and (4.2), we can see that introducing larger
distortions will increase the source variance for all future iterations, and thus increase the
communication cost. Therefore, using higher coding rates (introducing lower distortions) in
the initial iterations may be preferable in settings where the source variance makes a significant
contribution to the overall cost, as in the low-SNR setting explored here.
It seems intuitive to expect the MSE to decrease monotonically in iterations as t Ñ 8,
but this has been observed not to be true in general. If the target MSE˚ is much higher
than the MSEpD, Tq (3.32) achievable in the lossless case (i.e., Ri ptq Ñ 8, @i P t1, . . . , mu,
t P t0, . . . , T ´ 1u), then the MSE at each node might increase initially, followed by a rapid
decrease in the later iterations. This behavior is illustrated in Figure 4.5. We hypothesize that
this nonmonotonicity is an effect of the different treatment of states and quantization errors in
the state update. Recall the iteration we use (2.8), borrowed from Frasca et al. [14]:

ζpt ` 1q “ ζptq ` pΩ ´ Iq Qpζptqq.

This equation can be rewritten as (2.10),
ζpt ` 1q “ Ωζptq ` pΩ ´ Iq δptq,
where the reader is reminded that ζptq P RmL is the state supervector (3.7) that contains
the vector-valued states ζ i ptq at each node i P t1, . . . , mu, and δptq P RmL is the quantization
error (3.10) corresponding to each element of ζptq. Because the state ζptq and quantization
error δptq are incorporated into the state update differently, they have different impacts on the

67

MSE. In our case, this difference appears to be favorable in the sense that it allows the impact
of distortions to decay quickly.
Mathematically, the reason for the nonmonotone optimal MSE behavior may be discerned
from the state evolution equations. The reader is reminded of the MSE state evolution relationship (3.31),
˙
ˆ
˙
1 J
1 J
t
t
MSEi ptq “
I ´ 11 W Σe p0qW I ´ 11
m
m
ˆ
˙
ˆ
˙ff
t´1
ÿ
1 J
1 J
t´s´1
t´s´1
pW ´ IqΣe psqpW ´ IqW
I ´ 11 W
`
I ´ 11
,
m
m
s“0
„ˆ

(4.3)

ii

where the µe ptqµJ
e ptq term is omitted because the simulated signal and noise are zero-mean.
To simplify the explanation of the nonmonotone MSE behavior, we will consider the
constant-distortion optimization, for which the above (4.3) reduces to
„ˆ

˙
ˆ
˙
1 J
1 J
t
t
MSEi ptq “
I ´ 11 W Σe p0qW I ´ 11
m
m
ˆ
˙
ˆ
˙ff
t´1
ÿ
1 J
1
2
`
Dpsq I ´ 11 Wt´s´1 pW ´ Iq Wt´s´1 I ´ 11J
.
m
m
s“0
ii

In the Gaussian-signal-plus-noise setting (4.1) mentioned at the beginning of this chapter,
the initial error covariance matrix (3.26) is Σe p0q “ σn2 pI ´ m1 11J q. Because I ´ m1 11J is idempo`
˘
`
˘
tent,1 the first term in the above equation reduces to σn2 I ´ m1 11J W2t I ´ m1 11J , and we can
assess the different contributions of noise variance σn2 and distortion Dptq to the network MSE,
MSEpD, tq, by examining the traces of the matrices
ˆ

1 An

˙
ˆ
˙
1 J
1 J
2t
I ´ 11 W I ´ 11
m
m

idempotent matrix A is one for which A2 “ A [81].

68

Figure 4.6 Coefficients of the noise variance σn2 (solid blue curve) and distortion Dp0q (dashed red
curve) in the equation for the MSE (3.31) when the states are initialized as a Gaussian random vector
corrupted by AWGN (ρc “ 0.35, m “ 20).

69

and
ˆ

˙
˙
ˆ
1 J
1 J
2 t
t
I ´ 11 W pW ´ Iq W I ´ 11 ,
m
m

for t P t0, . . . , T ´ 1u. The reader is reminded that the trace of a matrix A P Rnˆn , written
řn
tr pAq, is the sum of its diagonal elements (i.e., tr pAq :“ i“1
aii ).
These diagonal components are the coefficients of the noise variance σn2 and the distortions
Di ptq in the equation for the MSE. Figure 4.6 plots the diagonal values for the above matrices,
@t P t0, . . . , T ´ 1u. The coefficients of Dptq decrease sharply after a single iteration, meaning
that the distortion contribution to the MSE vanishes quickly compared to the contribution of
the initial state variance.

4.2

Comparison to prior art

To compare our work to the prior art [20, 21], we generated 32 random geometric networks [45]
with connectivity radius2 ρc P t0.35, 0.45u on a two-dimensional unit torus. For each of these
networks, consensus was run on 1000 realizations of the initial states, which were length-10000
i.i.d. Gaussian vectors zi p0q “ x ` ni , @i P t1, . . . , mu, x „ N p0, Iq, ni „ N p0, 0.5Iq (SNR = 2,
which is 3.01 dB).
We simulated ProgQ [20] and order-one predictive coding [21], using initial rates Ri p0q P
t4, . . . , 7u @i and Ri p0q P t3, . . . , 6u @i, respectively. The measured final MSE (3.31) for these
schemes were set as the target values for the proposed GGP and heuristic optimizations.
For all schemes, MSEpD, Tq and Ragg were computed. These values were averaged over all
32 realizations of each (ρc , Ri p0q, T) setting, and the resulting averages were plotted against
each other. Assuming an additive noise model for the quantization error means that the
MSE for quantized communication will always be higher than the MSE of the unquantized
algorithm. We therefore introduce two terms to define the MSE performance relative to the
2 A random geometric network is one for which each node V P V is associated with a geographic coordinate v .
i
i
Nodes Vi , Vj are connected if kvi ´ v j k2 ď ρc [46].

70

ideal, unquantized algorithm. To compensate for the effect of network topology on the MSE,
we define the lossless MSE,
ˇ
ˇ
MSElossless ptq :“ MSEpD, tqˇ

D“0

.

We then define the excess MSE (EMSE) as

EMSE :“ 10 log10

MSEpD, Tq
,
MSElossless pTq

which represents the increase in MSE over lossless consensus resulting from distortion. The
EMSE was used for the generation of RD trade-off curves.

4.2.1

Implementation details

To generate the weight matrix W, the max-degree heuristic [37, 82] was used, so that

W “ I ´ αL “ I ´

0.9
m L,
maxi tdeg iui“1

(4.4)

where deg i is the degree of node i and L is the graph Laplacian [54, Sec. 10.9]. This weighting
scheme was selected for simplicity, following the implementation of Yildiz and Scaglione [21].
The factor of 0.9 is used to ensure numerical stability, since the coefficient α of L in (4.4) must
m q to guarantee convergence [37]. Based on (4.4), W is symmetric
satisfy α P p0, 1{ maxi tdeg iui“1

and satisfies all the requirements of the matrix properties discussed in Chapter 2.
The implementation of Yildiz and Scaglione’s order-one predictive coding [21], which was
graciously provided by the authors, was modified to use fixed-rate uniform quantization but
allow for the rate to vary with iteration and node indices. This capability was implemented by
running two rate update recursions—one to keep track of the ideal (real-valued) rates given by
the quantization noise variance recursion [21], and another to perform the predictive coding
using rates that were rounded to the nearest integral value. Interested readers can access the

71

code, which will be posted online, for further details.
The optimization modeler for the GGP and its associated solvers suffered from poor
stability at high iteration count, and poor scaling in the network size and number of iterations
for the exact program (3.38). Therefore, only the fixed-distortion problem (3.39) was solved for
comparison to the prior art.
The bin size of all fixed-rate uniform quantizers was set to 12 times the standard deviation
of the Gaussian data to prevent clipping.

4.2.2

Discussion of comparison results

To adequately discuss the RD results, it is first necessary to comment on some properties of
each of the schemes presented. The ProgQ algorithm of Thanou et al. [20] uses a time- and
node-invariant fixed-rate uniform quantizer (i.e., Ri ptq “ R, @i P t1, . . . , mu, t P t0, . . . , T ´ 1u),
whereas Yildiz and Scaglione [21] allow the use of different rates at each node and iteration.
Thanou et al. [20] use the same state update as ours (2.8), but Yildiz and Scaglione [21]
use a different update that is incapable of truly converging in the presence of quantization
error. The final asymptotic MSE for the predictive scheme depends on the sum of distortions
Di ptq, t P Zě0 . If these distortions are chosen to form a convergent series, then the MSE
will converge to a nonzero, but bounded, value. Because of this limitation, the predictive
scheme [21] is heavily dependent on the starting rates Ri p0q.
For each fixed-rate uniform quantizer sequence in the prior art comparisons, the result of
the heuristic from Sec. 3.4 was compared to the optimum obtained by an exhaustive search over
řT´1
all possible rate sequences with a maximum aggregate rate value per node of t“0
RGGP ptq.
In every case, the heuristic was able to find the optimal result. These results seem to justify our
heuristic search. Despite the reduced search complexity (exponential instead of combinatorial)
of this heuristic, the simulated performance was no worse. Therefore, at the least, this heuristic
seems to be a promising approach whenever fixed-rate coding is used. The RD performances
of the proposed optimization schemes are compared to the predictive coding scheme of Yildiz

72

Figure 4.7 RD trade-off curves for the proposed schemes versus order-one predictive coding [21].
All schemes were simulated for random geometric graphs on a unit torus with connectivity radius
ρc “ 0.35. Top: number of iterations T “ 5. Bottom: number of iterations T “ 7.

73

Figure 4.8 RD trade-off curves for the proposed schemes versus order-one predictive coding [21].
All schemes were simulated for random geometric graphs on a unit torus with connectivity radius
ρc “ 0.45. Top: number of iterations T “ 5. Bottom: number of iterations T “ 7.

74

and Scaglione [21] in Figures 4.7 and 4.8 and to the ProgQ algorithm of Thanou et al. [20] in
Figures 4.9 and 4.10.
These RD trade-off curves (Figures 4.7–4.10) contain the performance of both ECSQ and
fixed-rate uniform quantization. For both of our schemes, the predicted RD performance (as
computed by the state evolution equations (3.23)–(3.26)) is compared to the actual performance.
The distortion (measured by the EMSE) is given by the y-axis, and the aggregate rate Ragg
by the x-axis. The predicted performance is denoted by a dashed line, while the simulated
performance is represented as a solid line. The legend disambiguates the many curves that are
plotted. Alongside our approaches, we plot the RD performance of one of the comparators. A
curve closer to the bottom left corner of these figures indicates better performance, meaning
lower aggregate rate Ragg to achieve the same EMSE, or alternatively lower EMSE for a
particular Ragg investment.
In some cases, the predicted performance and measured performance of ECSQ do not
match. Because the ECSQ used in the simulations is not dithered, the additive quantization
model only holds approximately. As the aggregate rate increases, the performance improves,
and better adherence to predicted performance can be accomplished using dithering.
In addition to the ECSQ and fixed-rate uniform quantizers, Lloyd-Max quantization [60,
61] was studied. However, its empirical MSE performance did not match predictions, because
its quantization error δptq is correlated with the source ζptq. Due to this correlation, the use of
Lloyd-Max quantization in our framework seems inappropriate.

4.2.3

Comparison to predictive coding scheme

If we neglect the effect of entropy coding on the RD curves in Figures 4.7 and 4.8, then we
can more easily assess the impact of predictive coding on the aggregate rate. In each of these
cases, the predictive coding scheme requires lower aggregate rate than the proposed fixed-rate
uniform quantization strategy. However, when ECSQ is used in place of fixed-rate uniform
quantization, the results improve significantly.

75

Figure 4.9 RD trade-off curves for the proposed schemes versus ProgQ [20]. All schemes were simulated
for random geometric graphs on a unit torus with connectivity radius ρc “ 0.35. Top: number of
iterations T “ 5. Bottom: number of iterations T “ 7.

76

Figure 4.10 RD trade-off curves for the proposed schemes versus ProgQ [20]. All schemes were
simulated for random geometric graphs on a unit torus with connectivity radius ρc “ 0.45. Top: number
of iterations T “ 5. Bottom: number of iterations T “ 7.

77

The performance advantage of the predictive scheme stems from its use of linear minimum
MSE (LMMSE) prediction [35] to reduce the source variance prior to quantization. As the
connectivity parameter ρc and total iteration count T increase, there is no significant change
in performance relative to the proposed schemes. This lack of improvement may be due to
the previously discussed nonzero asymptotic MSE. Clearly, the advantages of the predictive
quantization strategy combined with entropy coding and the modified state update of Frasca
et al. [14] would produce better results, and it is worthwhile to explore whether LMMSE
prediction can be incorporated into our framework.

4.2.4

Comparison to the ProgQ scheme

Unlike the predictive scheme [21], ProgQ [20] converges to the true sample average in the
limit t Ñ 8. In the low-connectivity (ρc “ 0.35) setting, the ProgQ scheme is outdone by our
fixed-rate uniform quantizer. Increasing the connectivity to ρc “ 0.45 also increases the speed
of convergence. Because the ProgQ scheme uses differential coding, faster convergence leads
to lower source variance, and hence lower distortion for a given coding rate (3.36). Therefore,
the performance of the ProgQ scheme improves relative to our approaches with higher ρc .
Furthermore, from Figure 4.10, it is clear that this advantage becomes more significant as the
iteration count T increases. Even with constant coding rates, the ProgQ scheme can converge in
the limit; however, it is clear from Figure 4.2 that the optimal sequences for our schemes grow
with t. Therefore, asymptotically, we expect that ProgQ will outdo our proposed schemes,
despite our constrained optimization approach.

78

CHAPTER

5
CONCLUSION

In this thesis, we presented a new approach to reducing the communication cost in distributed
average consensus using lossy compression and generalized geometric programming (GGP)
modeling in the case of initial Gaussian states with entropy-coded scalar uniform quantization
(ECSQ), or the case of an arbitrary distribution using fixed-rate uniform quantization. The
main results in Chapter 3 of this thesis are (i) the proof that the problem of minimizing
the coding rate for a particular number of consensus iterations T, subject to a constraint
on the mean square error (MSE), can be posed as a GGP, and (ii) to a lesser extent, the
state evolution equations for consensus using Frasca’s modified iteration [14]. This analysis
allows the minimization of communication cost to be tackled with a structured and principled
approach that may be promising for future research.
Numerical results in Chapter 4 verify the state-evolution equations and hint at properties
of the optimal rate sequences in the absence of side information at the encoder and decoder.
Additionally, in the case that fixed-rate coding is used, a heuristic is presented in Section 3.5.2
that reduces the optimization complexity substantially. The main advantage of the proposed
approach is the ability to specify fine trade-offs among accuracy, run time, and communication
cost as measured by the aggregate coding rate.

79

From the comparisons to prior art in Chapter 4, it can be seen that even simple incorporation
of side information into the coding strategy can result in significant performance improvements
in terms of aggregate rate. Using differential coding, as in the approach of Thanou et al. [20], is
simple. Exploring its effect on performance would make an interesting future research topic.
Unfortunately, the optimization (3.38) does not scale well, and it may be worthwhile to
study the properties of the optimal rate sequences to design some theoretically motivated
heuristics, or to investigate alternative optimization implementations that do not rely on
explicit GGP modeling. In the case of uniform quantization with fixed-rate coding, numerical
results showcased the effectiveness of our proposed heuristic approach. For each simulation
conducted, the reduced complexity search heuristic found the same result as an exhaustive
search, with a significant reduction in run time.
On the applied front, it would be interesting to explore potential uses of this optimized form
of consensus in existing distributed algorithms. One potential path forward is to use consensus
to perform the pseudodata fusion step in row-partitioned multiprocessor approximate message
passing [47–51]. Because of the nonhomogeneous MSE evolution at each node, this becomes a
nontrivial problem. The distributed nature of signal acquisition in sensor networks may make
approximate message passing in the multimeasurement vector framework [83] appropriate.
Other potential areas for application include dictionary learning [9] and computer vision [8].
Another possible application of this work is in information theory. With relaxations of some
of the state-update assumptions, it may be possible to numerically compute the operational
rate distortion (RD) relationship for synchronous-time consensus for Gaussian data using our
optimization (3.38). This would represent an upper bound on the operational RD relationship
in consensus, and it might offer insight into the continuing information-theoretic investigation
of source coding for distributed average consensus.

80

BIBLIOGRAPHY
[1] G. J. Pottie and W. J. Kaiser, “Wireless integrated network sensors,” Commun. ACM, vol.
43, no. 5, pp. 51–58, May 2000.
[2] D. Estrin, D. Culler, K. Pister, and G. Sukhatme, “Connecting the physical world with
pervasive networks,” IEEE Pervasive Comput., vol. 1, no. 1, pp. 59–69, Jan. 2002.
[3]

Amazon EC2, https://aws.amazon.com/ec2/.

[4] N. Noorshams and M. J. Wainwright, “Non-asymptotic analysis of an optimal algorithm
for network-constrained averaging with noisy links,” IEEE J. Sel. Topics Signal Process.,
vol. 5, no. 4, pp. 833–844, Aug. 2011.
[5] A. G. Dimakis, S. Kar, J. M. F. Moura, M. G. Rabbat, and A. Scaglione, “Gossip algorithms
for distributed signal processing,” Proc. IEEE, vol. 98, no. 11, pp. 1847–1864, Nov. 2010.
[6] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and cooperation in networked
multi-agent systems,” Proc. IEEE, vol. 95, no. 1, pp. 215–233, Jan. 2007.
[7] M. Nokleby, W. U. Bajwa, R. Calderbank, and B. Aazhang, “Toward resource-optimal
consensus over the wireless medium,” IEEE J. Sel. Topics Signal Process., vol. 7, no. 2,
pp. 284–295, Apr. 2013.
[8] T. Tron and R. Vidal, “Distributed computer vision algorithms through distributed averaging,” in Proc. IEEE Conf. Comput. Vision and Pattern Recognition (CVPR), Providence,
RI, Jun. 2011, pp. 57–63.
[9] H. Raja and W. U. Bajwa, “Cloud K-SVD: A collaborative dictionary learning algorithm
for big, distributed data,” IEEE Trans. Signal Process., vol. 64, no. 1, pp. 173–188, Jan.
2016.
[10]

M. H. DeGroot, “Reaching a consensus,” J. Amer. Statist. Assoc., vol. 69, no. 345, pp. 118–
121, 1974.

[11]

V. Borkar and P. Varaiya, “Asymptotic agreement in distributed estimation,” IEEE Trans.
Autom. Control, vol. AC-27, no. 3, pp. 650–655, Jun. 1982.

[12]

J. N. Tsitsiklis, “Problems in decentralized decision making and computation,” PhD
thesis, Massachussetts Inst. Technol., Cambridge, MA, Nov. 1984.

[13]

J. Tsitsiklis, D. Bertsekas, and M. Athans, “Distributed asynchronous deterministic and
stochastic gradient optimization algorithms,” IEEE Trans. Autom. Control, vol. 31, no. 9,
pp. 803–812, Sep. 1986.

81

[14]

P. Frasca, R. Carli, F. Fagnani, and S. Zampieri, “Average consensus on networks with
quantized communication,” Int. J. Robust Nonlinear Control, vol. 19, no. 16, pp. 1787–1816,
Nov. 2008.

[15]

R. Carli, F. Bullo, and S. Zampieri, “Quantized average consensus via dynamic coding/decoding schemes,” Int. J. Robust Nonlinear Control, vol. 20, no. 2, pp. 156–175, May
2009.

[16]

M. E. Chamie, J. Liu, and T. Başar, “Design and analysis of distributed averaging with
quantized communication,” in Proc. 53rd IEEE Conf. Decision Control, Dec. 2014, pp. 3860–
3865.

[17]

R. Carli, F. Fagnani, P. Frasca, and S. Zampieri, “Efficient quantization for average
consensus,” Arxiv preprint arXiv:0903.1337, Mar. 2009.

[18]

J. Fang and H. Li, “An adaptive quantization scheme for distributed consensus,” in Proc.
IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), Apr. 2009, pp. 2777–2780.

[19]

T. C. Aysal, M. J. Coates, and M. G. Rabbat, “Distributed average consensus with dithered
quantization,” IEEE Trans. Signal Process., vol. 56, no. 10, pp. 4905–4918, Oct. 2008.

[20]

D. Thanou, E. Kokiopoulou, Y. Pu, and P. Frossard, “Distributed average consensus
with quantization refinement,” IEEE Trans. Signal Process., vol. 61, no. 1, pp. 194–205, Jan.
2013.

[21]

M. E. Yildiz and A. Scaglione, “Coding with side information for rate-constrained
consensus,” IEEE Trans. Signal Process., vol. 56, no. 8, pp. 3753–3764, Aug. 2008.

[22]

——, “Limiting rate behavior and rate allocation strategies for average consensus problems with bounded convergence,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process.
(ICASSP), Apr. 2008, pp. 2717–2720.

[23]

C. Mosquera, R. López-Valcarce, and S. K. Jayaweera, “Stepsize sequence design for
distributed average consensus,” IEEE Signal Process. Lett., vol. 17, no. 2, pp. 169–172, Feb.
2010.

[24]

R. Rajagopal and M. J. Wainwright, “Network-based consensus averaging with general
noisy channels,” IEEE Trans. Signal Process., vol. 59, no. 1, pp. 373–385, Jan. 2011.

[25]

A. Nedić, A. Olshevsky, A. Ozdaglar, and J. N. Tsitsiklis, “On distributed averaging
algorithms and quantization effects,” IEEE Trans. Autom. Control, vol. 54, no. 11, pp. 2506–
2517, Nov. 2009.

[26]

S. Zhu, S. Liu, J. Xu, Y. C. Soh, and L. Xie, “Averaging based distributed estimation
algorithm for rate-constrained sensor networks with additive quantization model,” in
Proc. IEEE Ann. Conf. Decision Control (CDC), Osaka, Japan, Dec. 2015, pp. 6245–6250.

82

[27]

S. Kar and J. M. F. Moura, “Distributed consensus algorithms in sensor networks:
Quantized data and random link failures,” IEEE Trans. Signal Process., vol. 58, no. 3,
pp. 1383–1400, Mar. 2010.

[28]

Y. Huang and Y. Hua, “On energy for progressive and consensus estimation in multihop
sensor networks,” IEEE Trans. Signal Process., vol. 59, no. 8, pp. 3863–3875, Aug. 2011.

[29]

T. Berger, Rate distortion theory: Mathematical basis for data compression. Englewood Cliffs,
NJ: Prentice-Hall, 1971.

[30]

T. M. Cover and J. A. Thomas, Elements of information theory. New York, NY, USA:
Wiley-Interscience, 1991.

[31]

S. Boyd and L. Vandenberghe, Convex optimization. Cambridge, UK: Cambridge University Press, 2004.

[32]

L. Xiao, S. Boyd, and S.-J. Kim, “Distributed average consensus with least-mean-square
deviation,” J. Parallel Distributed Comput., vol. 67, no. 1, pp. 33–46, Jan. 2007.

[33]

F. F. C. Rego, Y. Pu, A. Alessandretti, A. P. Aguiar, and C. N. Jones, “A consensus
algorithm for networks with process noise and quantization error,” in Proc. 53rd Annu.
Allerton Conf. on Commun., Control, and Comput., Allerton House, UIUC, IL, Oct. 2015,
pp. 488–495.

[34]

T. Li, M. Fu, L. Xie, and J.-F. Zhang, “Distributed consensus with limited communication
data rate,” IEEE Trans. Autom. Control, vol. 56, no. 2, pp. 279–292, Feb. 2011.

[35]

A. Gersho and R. M. Gray, Vector quantization and signal compression. Norwell, MA:
Kluwer, 1993.

[36]

M. E. Yildiz and A. Scaglione, “Differential nested lattice encoding for consensus
problems,” in Proc. 6th Int. Symp. Inform. Process. Sensor Networks (IPSN), Cambridge,
MA, Apr. 2007, pp. 89–98.

[37]

L. Xiao and S. Boyd, “Fast linear iterations for distributed averaging,” Syst. Control Lett.,
vol. 53, no. 1, pp. 65–78, Sep. 2004.

[38]

S. Zhu, Y. C. Soh, and L. Xie, “Distributed parameter estimation with quantized communication via running average,” IEEE Trans. Signal Process., vol. 63, no. 17, pp. 4634–4646,
Sep. 2015.

[39]

D. Thanou, H. Park, E. Kokiopoulou, and P. Frossard, “Polynomial filter design for
quantized consensus,” in Proc. 18th European Signal Process. Conf. (EUSIPCO), Aalborg,
Denmark, Aug. 2010, pp. 184–188.

83

[40]

J. Fang and H. Li, “Distributed consensus with quantized data via sequence averaging,”
IEEE Trans. Signal Process., vol. 58, no. 2, pp. 944–948, Feb. 2010.

[41]

O. Ayaso, D. Shah, and M. A. Dahleh, “Information theoretic bounds for distributed
computation over networks of point-to-point channels,” IEEE Trans. Inf. Theory, vol. 56,
no. 12, pp. 6020–6039, Dec. 2010.

[42]

A. Xu and M. Raginsky, “A new information-theoretic lower bound for distributed
function computation,” in Proc. IEEE Int. Symp. Inform. Theory (ISIT), Honolulu, HI, Jun.
2014, pp. 2227–2231.

[43]

H.-I. Su and A. El Gamal, “Distributed lossy averaging,” IEEE Trans. Inf. Theory, vol. 56,
no. 7, pp. 3422–3437, Jul. 2010.

[44]

Y. Yang, P. Grover, and S. Kar, “Rate distortion for lossy in-network function computation:
Information dissipation and sequential reverse water-filling,” IEEE Trans. Inf. Theory, vol.
PP, no. 99, pp. 1–29, May 2017.

[45]

M. Penrose, Random geometric graphs. New York, NY: Oxford University Press, 2003.

[46]

S. Boyd, A. Ghosh, B. Prabhakar, and D. Shah, “Mixing times for random walks on
geometric random graphs,” in ALENEX/ANALCO, Vancouver, British Columbia, Canada,
Jan. 2005, pp. 240–249.

[47]

J. Zhu, D. Baron, and A. Beirami, “Optimal trade-offs in multi-processor approximate
message passing,” Arxiv preprint arXiv:1601.03790, Nov. 2016.

[48]

J. Zhu, A. Beirami, and D. Baron, “Performance trade-offs in multi-processor approximate message passing,” in Proc. IEEE Int. Symp. Inform. Theory (ISIT), Barcelona, Spain,
Jul. 2016, pp. 680–684.

[49]

J. Zhu, “Statistical physics and information theory perspectives on linear inverse problems,” PhD thesis, North Carolina State University, Raleigh, NC, Jan. 2017.

[50]

P. Han, J. Zhu, R. Niu, and D. Baron, “Multi-processor approximate message passing using lossy compression,” in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP),
Shanghai, China, Mar. 2016, pp. 6240–6244.

[51]

J. Zhu, R. Pilgrim, and D. Baron, “An overview of multi-processor approximate message
passing,” in Proc. 51st IEEE Conf. Inform. Sci. Syst., Baltimore, MD, Mar. 2017.

[52]

B. Widrow and I. Kollár, Quantization noise: Roundoff error in digital computation, signal
processing, control, and communications. New York, NY: Cambridge University Press, 2008.

[53]

S. P. Lipshitz, R. A. Wannamaker, and J. Vanderkooy, “Quantization and dither: A
theoretical survey,” J. Audio Eng. Soc., vol. 40, no. 5, pp. 355–375, 1992.

84

[54]

A. Brandstädt, T. Nishizeki, K. Thulasiraman, and S. Arumugam, Eds., Handbook of graph
theory, combinatorial optimization, and algorithms, ser. Chapman and Hall/CRC Computer
and Information Science Series. Boca Raton, FL: Chapman and Hall/CRC, 2016, vol. 34.

[55]

C. E. Shannon, “A mathematical theory of communication,” Bell Syst. Tech. J., pp. 379–
423,623–656, 1948.

[56]

T. Wiegand and H. Schwarz. (Dec. 2012). Quantization, [Online]. Available: https :
//www.ic.tu-berlin.de/fileadmin/fg121/Source-Coding_WS12/05_QuantizationWS12.pdf.

[57]

S. Arimoto, “An algorithm for computing the capacity of an arbitrary discrete memoryless channel,” IEEE Trans. Inf. Theory, vol. 18, no. 1, pp. 14–20, Jan. 1972.

[58]

R. E. Blahut, “Computation of channel capacity and rate-distortion functions,” IEEE
Trans. Inf. Theory, vol. 18, no. 4, pp. 460–473, Jul. 1972.

[59]

K. Rose, “A mapping approach to rate-distortion computation and analysis,” IEEE Trans.
Inf. Theory, vol. 40, no. 6, pp. 1939–1952, Nov. 1994.

[60]

S. P. Lloyd, “Least squares quantization in PCM,” IEEE Trans. Inf. Theory, vol. 28, no. 2,
pp. 129–137, Mar. 1982.

[61]

J. Max, “Quantization for minimum distortion,” IRE Trans. Inf. Theory, vol. 6, no. 1,
pp. 7–12, Mar. 1960.

[62]

D. A. Huffman, “A method for the construction of minimum-redundancy codes,” Proc.
Inst. Radio Eng., vol. 9, no. 40, pp. 1098–1101, Sep. 1952.

[63]

H. Gish and J. Pierce, “Asymptotically efficient quantizing,” IEEE Trans. Inf. Theory, vol.
14, no. 5, pp. 676–683, Sep. 1968.

[64]

R. M. Gray and D. L. Neuhoff, “Quantization,” IEEE Trans. Inf. Theory, vol. IT-44,
pp. 2325–2383, Oct. 1998.

[65]

N. Farvardin and J. W. Modestino, “Optimum quantizer performance for a class of
non-gaussian memoryless sources,” IEEE Trans. Inf. Theory, vol. IT-30, no. 3, pp. 485–497,
May 1984.

[66]

R. Zamir and M. Feder, “On lattice quantization noise,” IEEE Trans. Inf. Theory, vol. 42,
no. 4, pp. 1152–1159, Jul. 1996.

[67]

A. B. Sripad and D. L. Snyder, “A necessary and sufficient condition for quantization
errors to be uniform and white,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-25,
no. 5, pp. 442–448, Oct. 1977.

85

[68]

E. W. Weisstein. (2017). Matrix direct sum. From MathWorld—a Wolfram Web Resource,
[Online]. Available: http://mathworld.wolfram.com/MatrixDirectSum.html.

[69]

P. Kabal, “Correlation properties of quantization noise,” in Proc. IEEE Int. Conf. Acoust.,
Speech, Signal Process. (ICASSP), Prague, Czech Republic, May 2011, pp. 5244–5247.

[70]

E. Aykol and K. Rose, “On constrained randomized quantization,” IEEE Trans. Signal
Process., vol. 61, no. 13, pp. 3291–3302, Jul. 2013.

[71]

S. Boyd, S.-J. Kim, L. Vandenberghe, and A. Hassibi, “A tutorial on geometric programming,” Optimization Eng., vol. 8, no. 67, pp. 67–127, Mar. 2007.

[72]

P. F. Panter and W. Dite, “Quantization distortion in pulse-count modulation with
nonuniform spacing of levels,” Proc. IRE, vol. 39, no. 1, pp. 44–48, Jan. 1951.

[73]

M. Grant and S. Boyd, CVX: MATLAB software for disciplined convex programming, version
2.1, http://cvxr.com/cvx, Mar. 2014.

[74]

——, “Graph implementations for nonsmooth convex programs,” in Recent Advances in
Learning and Control, ser. Lecture Notes in Control and Information Sciences, V. Blondel,
S. Boyd, and H. Kimura, Eds., http://stanford.edu/~boyd/graph_dcp.html, SpringerVerlag Limited, 2008, pp. 95–110.

[75]

E. Burnell and W. Hoburg, GPkit software for geometric programming, https://github.
com/hoburg/gpkit, Version 0.5.2, 2017.

[76]

S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling language for convex
optimization,” J. Mach. Learning Research, vol. 17, no. 83, pp. 1–5, 2016.

[77]

J. Löfberg, “YALMIP: A toolbox for modeling and optimization in MATLAB,” in Proc.
CACSD Conf., Taipei, Taiwan, 2004.

[78]

MATLAB, release 2017a. Natick, MA: The MathWorks, Inc., 2017.

[79]

E. W. Weisstein. (2017). Torus. From MathWorld—a Wolfram Web Resource, [Online].
Available: http://mathworld.wolfram.com/Torus.html.

[80]

——, (2017). Square torus. From MathWorld—a Wolfram Web Resource, [Online]. Available: http://mathworld.wolfram.com/SquareTorus.html.

[81]

K. Hoffman and R. Kunze, Linear algebra. Englewood Cliffs, NJ: Prentice Hall, Inc., 1971.

[82]

L. Xiao, S. Boyd, and S. Lall, “A scheme for robust distributed sensor fusion based on
average consensus,” in Proc. 4th Int. Symp. Inform. Process. Sensor Networks (IPSN), Apr.
2005, pp. 63–70.

86

[83]

J. Zhu, D. Baron, and F. Krzakala, “Performance limits for noisy multimeasurement
vector problems,” IEEE Trans. Signal Process., vol. 65, no. 9, pp. 2444–2454, May 2017.

87

