# Deep Equilibrium Models for Video Snapshot Compressive Imaging

Yaping Zhao
Westlake University
<EMAIL_ADDRESS>Siming Zheng
Computer Network Information Center
<EMAIL_ADDRESS>Xin Yuan
Westlake University
<EMAIL_ADDRESS>Corresponding author.

###### Abstract

The ability of snapshot compressive imaging (SCI) systems to efficiently
capture high-dimensional (HD) data has led to an inverse problem, which
consists of recovering the HD signal from the compressed and noisy
measurement. While reconstruction algorithms grow fast to solve it with the
recent advances of deep learning, the fundamental issue of accurate and stable
recovery remains. To this end, we propose deep equilibrium models (DEQ) for
video SCI, fusing data-driven regularization and stable convergence in a
theoretically sound manner. Each equilibrium model implicitly learns a
nonexpansive operator and analytically computes the fixed point, thus enabling
unlimited iterative steps and infinite network depth with only a constant
memory requirement in training and testing. Specifically, we demonstrate how
DEQ can be applied to two existing models for video SCI reconstruction:
recurrent neural networks (RNN) and Plug-and-Play (PnP) algorithms. On a
variety of datasets and real data, both quantitative and qualitative
evaluations of our results demonstrate the effectiveness and stability of our
proposed method. The code and models will be released to the public.

##  1 Introduction

Figure 1: Our proposed deep equilibrium models (DEQ) for SCI can lead to
stable recovery as $K$ increases, where $K$ denotes the iteration number
during the corresponding optimization progress. We test our model under two
different frameworks, _i.e_., RNN [7] and PnP-GAP [34], the fidelity and
stability of our model can be obviously observed.

Aiming at the efficient and effective acquisition of high-dimensional (HD)
visual signal, snapshot compressive imaging (SCI) systems have benefited from
the advent of novel optical designs to sample the HD data as two-dimensional
(2D) measurements. Considering the video SCI system, the 2D measurement of a
video, i.e., a three-dimensional (3D) data-cube leads to an inverse problem.
The goal of such an inverse problem is to recover a video from a collection of
noisy snapshots, which could be modeled as [33]:

$\displaystyle\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{x}+{\boldsymbol{e}},$
(1)

where $\boldsymbol{y}\in{\mathbb{R}}^{n}$ is the 2D measurement with $n$
equaling the number of each video frame’s pixels,
$\boldsymbol{\Phi}\in{\mathbb{R}}^{n\times nB}$ is the sensing matrix,
$\boldsymbol{x}\in{\mathbb{R}}^{nB}$ is the 3D data (by vectorizing each frame
and stacking them), and ${\boldsymbol{e}}$ is the measurement noise; here $B$
denotes that every $B$ video frames are collapsed into a single 2D
measurement. Though algorithms have been fully developed to reconstruct the
video from its snapshot measurement in recent years, the fundamental issue
remains: this inverse problem is inherently ill-posed, which makes the
recovery of the signal $x$ inaccurate and unstable for noise-affected data
$\boldsymbol{y}$ [12].

The rapid advancement of deep learning and artificial intelligence have
empowered a new wave of revolutionary solutions towards these previously
intractable problems. For instance, BIRNAT[7] employed recurrent neural
networks (RNNs) to reconstruct the video frames in a sequential manner and
explore the temporal correlation within the video SCI signal. Inspired by
particular optimization algorithms, GAP-net[16], DUN-3DUnet[28] designed deep
unfolding structures, which consist of a fixed number of architecturally
identical blocks. The heart of RNN and deep unfolding are deep neural
networks, which have posed new challenges due to their ever-growing depth and
huge training memory occupation. To overcome these difficulties, inspired by
[11], a recent work (RevSCI) [6] utilized reversible convolutional neural
networks to develop a memory-efficient structure. However, all of these
aforementioned algorithms inevitably suffer growing memory occupation with
increasing layer depth, and thus models need to be painstakingly designed.

Inspired by the plug-and-play (PnP) framework [23, 22] which has been proposed
for inverse problems with provable convergence [5, 21], PnP-FFDNet [34] and
PnP-FastDVDNet [35] bridged the gap between deep learning and conventional
optimization algorithms with the plug-and-play (PnP) framework, utilizing a
pre-trained denoiser as the proximal operator. While enjoying the advantages
of both data-driven regularization and flexible iterative optimization steps,
those algorithms still have hyperparameters to be tuned. Nevertheless, an
accurate result must be guaranteed with a proper parameter setting. Due to the
intrinsic unstable characteristic of the iterative recovery, even some
complicated strategy needs to be employed [27]. As we illustrate in Fig. 1 and
Fig. 2, the hyperparameters are unavoidable to be handcrafted to achieve
satisfactory performance in traditional algorithms.

An important and interesting research topic in deep learning is to train
arbitrary deep networks, in which the deep equilibrium models (DEQ) [3] stands
up as the leading method. A recent work [10] leverages DEQ to solve the
inverse problems in imaging, which corresponds to the potentially infinite
number of iteration steps in the PnP reconstruction scheme.

To accommodate the state-of-the-art SCI architectures and to enable low-memory
stable reconstruction, this paper sets about utilizing DEQ for solving the
inverse problem of video SCI. Specifically, we applied DEQ to two existing
models for video SCI reconstruction: RNN and PnP. Therefore, the former one is
equivalent to an infinite-depth network using only constant memory; the latter
one is tuning-free, and directly solves for the fixed point during the
iterative optimization process. On a variety of simulations and real datasets,
both quantitative and qualitative evaluations of our results demonstrate the
effectiveness of our proposed method. As shown in Fig. 2, our reconstruction
converges to stable results along with the increasing iterations during
optimization.

Figure 2: The quantitative comparison of different frameworks with or without
our proposed DEQ for SCI. The convergence trends of different algorithms
demonstrate that our model’s results can converge to a higher level.

In a nutshell, we aim to address the following two challenges which the SCI
reconstruction are facing while using deep neural network and iterative
optimization algorithms:

  * •

How deep should the model be? Can it be infinite?

  * •

Is there a tuning-free framework to be used? If yes, how to use it for SCI
reconstruction?

By employing the most recent development of DEQ, we demonstrate that the
answers to all the above questions are positive. Our specific contributions
are as follows:

  * 1)

We propose deep equilibrium models for video SCI, which fuses data-driven
regularization and stable convergence in a theoretically sound manner.

  * 2)

Each equilibrium model analytically computes the fixed point, thus enabling
unlimited iterative steps and infinite network depth with only a constant
memory requirement in training and testing.

  * 3)

We derive convergence theory for each equilibrium model, to ensure the
implicit operators in our models are nonexpansive.

  * 4)

On a variety of simulations and real datasets, both quantitative and
qualitative evaluations of our results demonstrate the effectiveness and
stability of our proposed method.

##  2 Related Work

###  2.1 Snapshot Compressive Imaging

The underlying principle of SCI is to compress the 3D data cube into a 2D
measurement by hardware, and then reconstruct the desired signal by
algorithms. Considering video SCI, it compresses the data-cube across the
temporal dimension, and thus enables a low-speed camera to capture high-speed
scenes. For instance, Llull et al. [15] proposed the coded aperture
compressive temporal imaging (CACTI) system, which decomposes the 3D cube into
its constituent 2D frames and imposes 2D masks for modulation.

Given the masks and measurements, plenty of algorithms including conventional
optimization[14, 29, 30, 32], end-to-end deep learning [19], deep unfolding
[16, 28] and plug-and-play [34, 35] are proposed for reconstruction. To solve
the ill-posed problem in Eq. (1), additional regularization is usually needed
to ensure accurate and stable recovery with respect to noise perturbation. To
this end, these algorithms obtain the estimated value $\hat{\boldsymbol{x}}$
of $\boldsymbol{x}$ by solving the following problem:

$\displaystyle\hat{\boldsymbol{x}}=\textstyle\mathop{\arg\min}_{\boldsymbol{x}}\frac{1}{2}||\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{x}||_{2}^{2}+R(\boldsymbol{x}),$
(2)

where $||\boldsymbol{y}-\mathbf{\Phi}\boldsymbol{x}||_{2}^{2}$ is the fidelity
term and $R(\boldsymbol{x})$ is the regularization term.

By introducing an auxiliary parameter $\boldsymbol{v}$ , the unconstrained
optimization in Eq. (2) can be converted into:

$\displaystyle(\boldsymbol{x},\boldsymbol{v})=\mathop{\arg\min}_{\boldsymbol{x},\boldsymbol{v}}\frac{1}{2}||\boldsymbol{y}-\Phi\boldsymbol{x}||_{2}^{2}+R(\boldsymbol{v}),~{}s.t.~{}\boldsymbol{x}=\boldsymbol{v}.$
(3)

Using the alternating direction method of multipliers (ADMM) [4] and
introducing another parameter ${\boldsymbol{u}}$, Eq. (3) could be divided
into the following sequence of sub-problems:

$\displaystyle\boldsymbol{x}^{(k+1)}$
$\displaystyle=\mathop{\arg\min}_{\boldsymbol{x}}\frac{1}{2}\|\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{x}\|_{2}^{2}+\frac{\rho}{2}\|\boldsymbol{x}-(\boldsymbol{v}^{(k)}-\frac{1}{\rho}{\boldsymbol{u}}^{(k)})\|^{2}_{2},$
(4) $\displaystyle\boldsymbol{v}^{(k+1)}$
$\displaystyle=\textstyle\mathop{\arg\min}_{\boldsymbol{v}}\mu
R(\boldsymbol{v})+\frac{\rho}{2}||\boldsymbol{v}-(\boldsymbol{x}^{(k)}+\frac{1}{\rho}{\boldsymbol{u}}^{(k)})||_{2}^{2},$
(5) $\displaystyle{\boldsymbol{u}}^{(k+1)}$
$\displaystyle=\textstyle{\boldsymbol{u}}^{(k)}+\rho(\boldsymbol{x}^{(k+1)}-\boldsymbol{v}^{(k+1)}),$
(6)

where the superscript k denotes the iteration number; $\rho$ is the penalty
parameter and $\mu$ is the regularization weight. Since Eq. (5) can be
regarded as a denoising process of $\boldsymbol{v}$, implicitly we have:

$\displaystyle\boldsymbol{v}^{(k+1)}=\textstyle\mathcal{D}^{(k+1)}(\boldsymbol{x}^{(k+1)}+\frac{1}{\rho}{\boldsymbol{u}}^{(k)}),$
(7)

where $\mathcal{D}$ is a denoiser.

On the other hand, generalized alternating projection (GAP) [13] can be used
as a (little bit) lower computational workload algorithm with the following
two steps:

$\displaystyle\boldsymbol{x}^{(k+1)}$
$\displaystyle=\boldsymbol{v}^{(k)}+{\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}(\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{v}^{(k)}),$
(8) $\displaystyle\boldsymbol{v}^{(k+1)}$
$\displaystyle=\mathcal{D}^{(k+1)}(\boldsymbol{x}^{(k+1)}).$ (9)

Eq. (8) can be solved efficiently due to the special structure of
$\boldsymbol{\Phi}$ in SCI [12].

###  2.2 Deep Unfolding

Inspired by optimization algorithms such as ADMM [4] and GAP [13], deep
unfolding methods [16, 28] are proposed to solve inverse problems in SCI,
which consist of a fixed number of architecturally identical blocks. Each of
those blocks represents a single iterative step in conventional optimization
algorithms. Though deep unfolding successfully assimilate the advantages of
the iterative optimization algorithms and could be trained in an end-to-end
manner, the fixed number of network blocks in deep unfolding is needed to be
kept small for two reasons: $i$) these systems should be concise to keep a
high inference speed for real-time reconstruction; $ii$) it is challenging to
train deep unfolding networks for numerous stages due to memory limitations.

###  2.3 Plug-and-Play

The latest trend is to bridge the gap between deep learning and optimization
with the PnP framework. Yuan et al. [35] proposed PnP-ADMM framework and PnP-
GAP framework, using a pre-trained denoiser as the proximal operator in Eq.
(5) and Eq. (9), respectively. In contrast to deep unfolding, PnP relieves
itself from the limited memory by integrating a flexible denoising module into
the iterative optimization process. Nevertheless, it suffers manual parameter
tuning in addition to the time-consuming reconstruction process. That is, its
performance is highly sensitive to the internal parameter selection, including
but not limited to the penalty parameter, the denoising level, and the
terminal step number. Moreover, the optimal parameter setting differs image-
by-image, depending on the modulation process, noise level, noise type, and
the unknown image itself.

###  2.4 Memory-Efficient Deep Networks

Since the important factor that limits the development of deep learning and
deep unfolding for SCI is limited memory on hardware devices used for
training, to address this issue, RevSCI [6] developed a memory-efficient
network for large-scale video SCI. Using reversible neural networks, where
each layer’s input can be calculated from the layer’s activation during back-
propagation, which means the activation during training is not needed to be
stored for sake of saving memory. Nevertheless, it still suffers growing
memory occupation along with the increasing depth of the network. In contrast,
DEQ reduces memory consumption to a constant (i.e., independent of network
depth) by directly differentiating through the equilibrium point and thus
circumvents the construction and maintenance of layers.

###  2.5 Deep Equilibrium Models

Motivated by the surprisingly recent works [1, 9, 8] that employ the same
transformation in each layer and still achieve competitive results with the
state-of-the-art, Bai et al. [2] proposed a new approach to model this process
and directly computed the fixed point. To leverage ideas from DEQ, Gilton et
al. [10] proposed DEQ for inverse problems in imaging, which corresponds to a
potentially infinite number of iteration steps in the PnP reconstruction
scheme. In this paper, we present a novel approach for video SCI using DEQ,
taking both the PnP and the RNN framework into considerations.

##  3 Method

Given measurement $\boldsymbol{y}\in{\mathbb{R}}^{n}$ with compression rate
$B$ and sensing matrix $\boldsymbol{\Phi}\in{\mathbb{R}}^{n\times nB}$ as
input, we consider an optimization iteration or neural network as:

$\displaystyle\boldsymbol{x}^{(k+1)}=f_{\theta}(\boldsymbol{x}^{(k)};\boldsymbol{y},\boldsymbol{\Phi}),\quad
k=0,1,\dots,\infty,$ (10)

where $\theta$ denotes the weights of embedded neural networks;
$\boldsymbol{x}^{(k)}\in{\mathbb{R}}^{nB}$ is the output of the $k^{th}$
iterative step or hidden layer, and
$\boldsymbol{x}^{(0)}=\mathbf{\boldsymbol{\Phi}}^{\top}\boldsymbol{y}$;
$f_{\theta}(\cdot\ ;\boldsymbol{y},\boldsymbol{\Phi})$ is an iteration map
${\mathbb{R}}^{nB}\rightarrow{\mathbb{R}}^{nB}$ towards a stable equilibrium:

$\displaystyle\lim_{k\rightarrow+\infty}\boldsymbol{x}^{(k)}=\lim_{k\rightarrow+\infty}f_{\theta}(\boldsymbol{x}^{(k)};\boldsymbol{y},\boldsymbol{\Phi})\equiv\hat{\boldsymbol{x}}=f_{\theta}(\hat{\boldsymbol{x}};\boldsymbol{y},\boldsymbol{\Phi}),$
(11)

where $\hat{\boldsymbol{x}}\in{\mathbb{R}}^{nB}$ denotes the fixed point and
reconstruction result.

In Sec. 3.1, we first design different $f_{\theta}$ for SCI, in terms of the
implicit infinite-depth RNN architecture and infinitely iterative PnP
framework. Following [10], we utilize Anderson acceleration [24] to compute
the fixed point of $f_{\theta}$ efficiently. For gradient calculation, we
optimize the network wights $\theta$ by approximating the inverse Jacobian,
described in Sec. 3.2. Convergence of this scheme for specific $f_{\theta}$
designs is discussed in Sec. 3.3.

###  3.1 Forward Pass

Unlike the conventional optimization method where the terminal step number is
manually chosen or a network where the output is the activation from the
limited layers, the result of DEQ is the equilibrium point itself. Therefore,
the forward evaluation could be any procedure that solves for this equilibrium
point. Considering SCI reconstruction, we design novel iterative models that
converge to equilibrium.

####  3.1.1 Recurrent Neural Networks

Figure 3: Illustration of our proposed DEQ for SCI under the framework of
recurrent neural network (RNN), _i.e_., DE-RNN.

To achieve integration of DEQ and RNN for video SCI, we have:

$\displaystyle\boldsymbol{x}^{(k+1)}={\rm
RNN}_{\theta}(\boldsymbol{x}^{(k)},\boldsymbol{y},\boldsymbol{\Phi}),$ (12)

where ${\rm RNN}(\cdot\ )$ is a trainable RNN network learning to iteratively
reconstruct effective and stable data. And as shown in Fig 3, the
corresponding iteration map is:

$\displaystyle
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})={\rm
RNN}_{\theta}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{\Phi}).$ (13)

####  3.1.2 Generalized Alternating Projection

Figure 4: Illustration of our proposed DEQ for SCI under the framework of
generalized alternating projection (GAP), _i.e_., DE-GAP.

Regarding the optimization iterations in the GAP method, represented in Eq.
(8)-(9), we iteratively update $\boldsymbol{x}$ with:

$\displaystyle\boldsymbol{x}^{(k+1)}=\textstyle\mathcal{D}^{(k+1)}_{\theta}\left[\boldsymbol{x}^{(k)}+{\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}(\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{x}^{(k)})\right].$
(14)

Therefore, as illustrated in Fig. 4, the iteration map is:

$\displaystyle
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})=\textstyle\mathcal{D}_{\theta}(\boldsymbol{x}+{\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}(\boldsymbol{y}-\boldsymbol{\Phi}\boldsymbol{x})).$
(15)

####  3.1.3 Anderson Acceleration

To enforce fixed-point iterations converge more quickly, we make full use of
the ability to accelerate inference with standard fixed-point accelerators,
e.g., Anderson accelerator. Anderson acceleration utilizes previous iterations
to seek promising directions to move forward. Under the setting of Anderson
accelerator, we identify a vector
$\boldsymbol{\alpha}^{(k)}\in\mathbb{R}^{s}$, for $\delta>0$:

$\displaystyle\begin{aligned}
\boldsymbol{x}^{(k+1)}=&\textstyle(1-\delta)\sum_{i=0}^{s-1}\mathbf{\boldsymbol{\alpha}}_{i}^{(k)}\boldsymbol{x}^{(k-i)}\\\
&+\textstyle\delta\sum_{i=0}^{s-1}\mathbf{\boldsymbol{\alpha}}_{i}^{(k)}f_{\theta}(\boldsymbol{x}^{(k-i)};\boldsymbol{y},\boldsymbol{\Phi}),\end{aligned}$
(16)

where the vector $\mathbf{\boldsymbol{\alpha}}_{i}^{(k)}$ is the solution to
the optimization problem:

$\displaystyle\textstyle\operatorname*{arg\,min}_{\mathbf{\boldsymbol{\alpha}}}||{\bf
A}\boldsymbol{\alpha}||_{2}^{2},\quad
s.t.\quad\mathbf{1}^{\top}\boldsymbol{\alpha}=1,$ (17)

where ${\bf A}$ is a matrix whose $i$-th column is the vectorized residual
$f_{\theta}(\boldsymbol{x}^{(k-i)};\boldsymbol{y},\boldsymbol{\Phi})-\boldsymbol{x}^{(k-i)}$,
with $i=0,\dots,s-1$. When $s$ is small (e.g., $s=3$), the optimization
problem in Eq. (17) introduces trivial computation.

###  3.2 Backward Pass

While previous work often utilizes Newton’s method to achieve the equilibrium
and then backpropagate through all the Newton iterations, following [10], we
alternatively adopt another method with high efficiency and constant memory
requirement.

####  3.2.1 Loss Function

To optimize network parameters $\theta$, stochastic gradient descent is used
to minimize a loss function as below:

$\displaystyle\textstyle\theta^{*}=\operatorname*{arg\,min}_{\theta}\frac{1}{m}\sum_{i=1}^{m}\ell(f_{\theta}(\hat{\boldsymbol{x}}_{i};\boldsymbol{y}_{i},\boldsymbol{\Phi}_{i}),\boldsymbol{x}^{\star}_{i}),$
(18)

where $m$ is the number of training samples; $\ell(\cdot,\cdot)$ is a given
loss function, $\boldsymbol{x}^{\star}_{i}$ is the ground truth 3D data of the
$i$-th training sample, $\boldsymbol{y}_{i}$ is the paired measurement,
$\boldsymbol{\Phi}_{i}$ denotes the sensing matrix, and
$f_{\theta}(\hat{\boldsymbol{x}}_{i};\boldsymbol{y}_{i},\boldsymbol{\Phi}_{i})$
denotes the reconstruction result given as the fixed point
$\hat{\boldsymbol{x}}$ of the iteration map $f_{\theta}(\cdot\
;\boldsymbol{y},\boldsymbol{\Phi})$, as derived from Eq. (11). The mean-
squared error (MSE) loss is used for our video SCI reconstruction:

$\displaystyle\textstyle\ell(\hat{\boldsymbol{x}},\boldsymbol{x}^{\star})=\frac{1}{2}||\hat{\boldsymbol{x}}-\boldsymbol{x}^{\star}||_{2}^{2}.$
(19)

Since the reconstruction result is a fixed point of the iteration map
$f_{\theta}(\cdot\ ;\boldsymbol{y},\boldsymbol{\Phi})$, gradient calculation
of this loss term could be designed to avoid large memory demand. Following
[10], we calculate the gradient of the loss term, which takes the network
parameters $\theta$ into consideration.

####  3.2.2 Gradient Calculation

Let $\ell$ be an abbreviation of
$\ell(\hat{\boldsymbol{x}},\boldsymbol{x}^{\star})$ in Eq. (19), then the loss
gradient is:

$\displaystyle\begin{aligned}
\textstyle\frac{\partial\ell}{\partial\theta}&={\left(\frac{\partial\hat{\boldsymbol{x}}}{\partial\theta}\right)}^{\top}\frac{\partial\ell}{\partial\hat{\boldsymbol{x}}}&={\left(\frac{\partial\hat{\boldsymbol{x}}}{\partial\theta}\right)}^{\top}{\left(\hat{\boldsymbol{x}}-\boldsymbol{x}^{\star}\right)},\end{aligned}$
(20)

where $\frac{\partial\hat{\boldsymbol{x}}}{\partial\theta}$ is the Jacobian of
$\hat{\boldsymbol{x}}$ evaluated at $\theta$, and
$\frac{\partial\ell}{\partial\hat{\boldsymbol{x}}}$ is the gradient of $\ell$
evaluated at $\boldsymbol{x}^{\star}$.

Then to compute the Jacobian
$\frac{\partial\hat{\boldsymbol{x}}}{\partial\theta}$, we recall the fixed
point equation
$\hat{\boldsymbol{x}}=f_{\theta}(\hat{\boldsymbol{x}};\boldsymbol{y},\boldsymbol{\Phi})$
in Eq. (11). By implicitly differentiating both sides of this fixed point
equation, the Jacobian $\frac{\partial\hat{\boldsymbol{x}}}{\partial\theta}$
is solved as:

$\displaystyle\textstyle\frac{\partial\hat{\boldsymbol{x}}}{\partial\theta}=\left[\boldsymbol{I}-\left.\frac{\partial
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\boldsymbol{x}}\right|_{\boldsymbol{x}=\hat{\boldsymbol{x}}}\right]^{-1}\frac{\partial
f_{\theta}(\hat{\boldsymbol{x}};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\theta},$
(21)

which could be plugged into Eq. (20) and thus get:

$\displaystyle\begin{aligned}
\textstyle\frac{\partial\ell}{\partial\theta}=\left[\frac{\partial
f_{\theta}(\hat{\boldsymbol{x}};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\theta}\right]^{\top}\\!{\left[\boldsymbol{I}-\left.\frac{\partial
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\boldsymbol{x}}\right|_{\boldsymbol{x}=\hat{\boldsymbol{x}}}\right]^{-\top}}{(\hat{\boldsymbol{x}}-\boldsymbol{x}^{\star})},\end{aligned}$
(22)

where -⊤ denotes the inversion followed by transpose. As this method converted
gradient calculation to the problem of calculating an inverse Jacobian-vector
product, it avoids the backpropagation through many iterations of
$f_{\theta}(\hat{\boldsymbol{x}};\boldsymbol{y},\boldsymbol{\Phi})$. To
approximate the inverse Jacobian-vector product, we define the vector
$\boldsymbol{a}^{(\infty)}$ as:

$\displaystyle\textstyle\boldsymbol{a}^{(\infty)}={\left[\boldsymbol{I}-\left.\frac{\partial
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\boldsymbol{x}}\right|_{\boldsymbol{x}=\hat{\boldsymbol{x}}}\right]^{-\top}}\\!{(\hat{\boldsymbol{x}}-\boldsymbol{x}^{\star})}.$
(23)

Following [10], it is noted that $\boldsymbol{a}^{(\infty)}$ is a fixed point
of the equation:

$\displaystyle\begin{aligned}
\boldsymbol{a}^{(k+1)}=\textstyle{\left[\left.\frac{\partial
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\boldsymbol{x}}\right|_{\boldsymbol{x}=\hat{\boldsymbol{x}}}\right]^{-\top}}\boldsymbol{a}^{(k)}+{(\hat{\boldsymbol{x}}-\boldsymbol{x}^{\star})},\\\
k=0,1,\dots,\infty.\end{aligned}$ (24)

Therefore, the same algorithm used to calculate the fixed point
$\hat{\boldsymbol{x}}$ could also be used to calculate
$\boldsymbol{a}^{(\infty)}$ . The limit of fixed-point iterations for solving
Eq. (24) with initial iterate $\boldsymbol{a}^{(0)}=\mathbf{0}$ is denoted
equivalently to the Neumann series:

$\displaystyle\boldsymbol{a}^{(\infty)}=\textstyle\sum_{p=0}^{\infty}\left\\{{\left[\left.\frac{\partial
f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\boldsymbol{x}}\right|_{\boldsymbol{x}=\hat{\boldsymbol{x}}}\right]^{\top}}\right\\}^{p}\\!{(\hat{\boldsymbol{x}}-\boldsymbol{x}^{\star})}.$
(25)

To quickly calculate the vector-Jacobian products in Eq. (24) and Eq. (25), a
lot of auto-differentiation tools (e.g., autograd packages in Pytorch[17])
could be utilized. After the accurate approximation of
$\boldsymbol{a}^{(\infty)}$ is calculated, the gradient in Eq. (20) is given
by:

$\displaystyle\textstyle\frac{\partial\ell}{\partial\theta}={\left(\frac{\partial
f_{\theta}(\hat{\boldsymbol{x}};\boldsymbol{y},\boldsymbol{\Phi})}{\partial\theta}\right)}^{\top}\boldsymbol{a}^{(\infty)}.$
(26)

###  3.3 Convergence Theory

Given the iteration map $f_{\theta}(\cdot\
;\boldsymbol{y},\boldsymbol{\Phi}):{\mathbb{R}}^{nB}\rightarrow{\mathbb{R}}^{nB}$,
in this section, we discuss conditions that guarantee the convergence of the
proposed deep equilibrium models
$\boldsymbol{x}^{(k+1)}=f_{\theta}(\boldsymbol{x}^{(k)};\boldsymbol{y},\boldsymbol{\Phi})$
to a fixed-point $\hat{\boldsymbol{x}}$ as $k\rightarrow\infty$.

###### Theorem 1

(Convergence of DE-RNN). If there exists a constant $0\leq c<1$ satisfies
that:

$\displaystyle||{\rm
RNN}_{\theta}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{\Phi})-{\rm
RNN}_{\theta}(\boldsymbol{x}^{\prime},\boldsymbol{y},\boldsymbol{\Phi})||\leq
c||\boldsymbol{x}-\boldsymbol{x}^{\prime}||,$ (27)

then the DE-RNN iteration map
$f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})$ is nonexpansive.

Following [21], we assume that for all
$\boldsymbol{x},\boldsymbol{x}^{\prime}\in{\mathbb{R}}^{nB}$, there exists a
$\varepsilon>0$ such that the denoiser $\mathcal{D}$ satisfies:

$\displaystyle\|(\mathcal{D}_{\theta}-{\bf
I})(\boldsymbol{x})-(\mathcal{D}_{\theta}-{\bf
I})(\boldsymbol{x}^{\prime})\|\leq\varepsilon||\boldsymbol{x}-\boldsymbol{x}^{\prime}||,$
(28)

where $(\mathcal{D}_{\theta}-{\bf
I})(\boldsymbol{x}):=\mathcal{D}_{\theta}(\boldsymbol{x})-\boldsymbol{x}$,
that is, we assume the map $\mathcal{D}_{\theta}-{\bf I}$ is
$\varepsilon$-Lipschitz.

###### Theorem 2

(Convergence of DE-GAP). Under the assumption in Eq. (28), the DE-GAP
iteration map $f_{\theta}(\cdot;\boldsymbol{y},\boldsymbol{\Phi})$ defined in
Eq. (15) satisfies:

$\displaystyle\textstyle\|f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})-f_{\theta}(\boldsymbol{x}^{\prime};\boldsymbol{y},\boldsymbol{\Phi})\|\leq\eta\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\|$
(29)

for all $\boldsymbol{x},\boldsymbol{x}^{\prime}\in{\mathbb{R}}^{nB}$. The
coefficient $\eta$ is less than $1$ and is related to the eigenvalues of
${\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}\boldsymbol{\Phi}$,
in which case the DE-GAP iteration map
$f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})$ is contractive.

Proof. Based on the assumption in Eq. (28) and the DE-GAP iteration map
$f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})$ in Eq. (15), the
Jacobian of $f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})$ with
respect to $\boldsymbol{x}\in{\mathbb{R}}^{nB}$, denoted as
$\partial_{\boldsymbol{x}}f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})$,
is given by:

$\displaystyle\partial_{\boldsymbol{x}}f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})$
$\displaystyle=\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})({\bf
I}-{\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}\boldsymbol{\Phi}),$
(30)

where $\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}\in{\mathbb{R}}^{nB\times
nB}$ is the Jacobian of
$\mathcal{D_{\theta}}:{\mathbb{R}}^{nB}\rightarrow{\mathbb{R}}^{nB}$ with
respect to $\boldsymbol{x}\in{\mathbb{R}}^{nB}$. To prove
$f_{\theta}(\cdot;\boldsymbol{y},\boldsymbol{\Phi})$ is nonexpansive it
suffices to show
$||\partial_{\boldsymbol{x}}f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})||<1$
for all $\boldsymbol{x}\in{\mathbb{R}}^{nB}$, where $||\cdot||$ denotes the
spectral norm.

Following the derivation in [35], we define ${{\bf
Q}}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\top}$, which is a diagonal matrix. In
the following, we have:

$\displaystyle||\partial_{\boldsymbol{x}}f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})||$
$\displaystyle=$ $\displaystyle\
||\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})({\bf
I}-{\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}\boldsymbol{\Phi})||$
$\displaystyle=$ $\displaystyle\
||\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x}){\boldsymbol{\Phi}}^{\top}{{\bf
Q}}\boldsymbol{\Phi}||$ $\displaystyle=$ $\displaystyle\
||\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf I}+{\bf
I}-\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x}){\boldsymbol{\Phi}}^{\top}{{\bf
Q}}\boldsymbol{\Phi}||$ $\displaystyle=$ $\displaystyle\
||(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf I})+{\bf
I}-(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf I}+{\bf
I}){\boldsymbol{\Phi}}^{\top}{{\bf Q}}\boldsymbol{\Phi}||$ $\displaystyle=$
$\displaystyle\
||(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf I})+{\bf
I}-(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf
I}){\boldsymbol{\Phi}}^{\top}{{\bf
Q}}\boldsymbol{\Phi}-{\boldsymbol{\Phi}}^{\top}{{\bf Q}}\boldsymbol{\Phi}||$
$\displaystyle=$ $\displaystyle\
||(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf I})({\bf
I}-{\boldsymbol{\Phi}}^{\top}{{\bf Q}}\boldsymbol{\Phi})+{\bf
I}-{\boldsymbol{\Phi}}^{\top}{{\bf Q}}\boldsymbol{\Phi}||$ $\displaystyle=$
$\displaystyle\ ||({\bf
I}+[(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf
I})])({\bf I}-{\boldsymbol{\Phi}}^{\top}{{\bf Q}}\boldsymbol{\Phi})||$
$\displaystyle\leq$ $\displaystyle||({\bf
I}+[(\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf
I})])||\cdot||({\bf I}-{\boldsymbol{\Phi}}^{\top}{{\bf
Q}}\boldsymbol{\Phi})||$ $\displaystyle\leq$
$\displaystyle\textstyle(1+\varepsilon)\max_{i}|1-\lambda_{i}|,$ (31)

where $\lambda_{i}$ are eigenvalues of
${\boldsymbol{\Phi}}^{\top}(\boldsymbol{\Phi}{\boldsymbol{\Phi}}^{\top})^{-1}\boldsymbol{\Phi}$;
and the inequality Eq. (31) is based on the assumption that the map
$(\mathcal{D}_{\theta}-{\bf
I})(\boldsymbol{x}):=\mathcal{D}_{\theta}(\boldsymbol{x})-\boldsymbol{x}$ is
$\varepsilon$-Lipschitz. Therefore the spectral norm of its Jacobian
$\partial_{\boldsymbol{x}}\mathcal{D_{\theta}}(\boldsymbol{x})-{\bf I}$ is
bounded by $\eta,\eta=(1+\varepsilon)\max_{i}|1-\lambda_{i}|$.

Finally, under our assumption $\varepsilon>0$, we can achieve:

$\displaystyle\|\partial_{\boldsymbol{x}}f_{\theta}(\boldsymbol{x};\boldsymbol{y},\boldsymbol{\Phi})\|\leq(1+\varepsilon)\max_{i}|1-\lambda_{i}|.$
(32)

This demonstrates $f_{\theta}$ is $\eta$-Lipschitz with
$\eta=(1+\varepsilon)\max_{i}|1-\lambda_{i}|$.

Table 1:  The results in terms of PSNR (dB) and SSIM by different algorithms
on classical six datasets for video SCI reconstruction. Compared methods
include GAP-net [16], GAP-TV [32], E2E-CNN [19] and PnP-FFDnet [34].

Methods | Kobe | Traffic | Runner | Drop | Vehicle | Aerial | Average
---|---|---|---|---|---|---|---
GAP-net-AE-S9 | 24.20, 0.570 | 21.13, 0.685 | 29.18, 0.886 | 32.21, 0.907 | 24.19, 0.769 | 24.41, 0.744 | 25.89, 0.760
GAP-TV | 26.46, 0.885 | 20.89, 0.715 | 28.52, 0.909 | 34.63, 0.970 | 24.82 0.838 | 25.05, 0.828 | 26.73, 0.858
E2E-CNN | 29.02, 0.861 | 23.45, 0.838 | 34.43, 0.958 | 36.77, 0.974 | 26.40, 0.886 | 27.52, 0.882 | 29.26, 0.900
PnP-FFDnet | 30.50, 0.926 | 24.18, 0.828 | 32.15, 0.933 | 40.70, 0.989 | 25.42, 0.849 | 25.27, 0.829 | 29.70, 0.892
DE-RNN | 21.46, 0.697 | 19.47, 0.715 | 27.85, 0.818 | 30.16, 0.909 | 23.65, 0.832 | 24.83, 0.855 | 24.53, 0.804
DE-GAP-Unet-3D | 26.76, 0.866 | 21.42, 0.786 | 30.45, 0.894 | 33.82, 0.963 | 24.94, 0.885 | 24.83, 0.847 | 27.07, 0.878
DE-GAP-RSN-CNN | 27.33, 0.887 | 22.58, 0.829 | 30.74, 0.903 | 35.95.0.977 | 25.33, 0.899 | 25.57, 0.881 | 27.92, 0.896
DE-GAP-RSN-Unet | 28.92, 0.939 | 23.68, 0.869 | 32.37, 0.951 | 36.54, 0.972 | 25.50, 0.905 | 25.67, 0.884 | 28.80, 0.913
DE-GAP-CNN | 28.79, 0.935 | 23.55, 0.864 | 32.35, 0.950 | 38.14, 0.983 | 25.45, 0.903 | 25.84, 0.890 | 29.02, 0.921
DE-GAP-FFDnet | 29.32, 0.952 | 24.71, 0.907 | 33.06, 0.971 | 39.89, 0.992 | 25.85, 0.905 | 26.02, 0.892 | 29.81, 0.936

Figure 5: Comparison of selected reconstruction results with the spatial size
of $256\times 256\times 8$. It can be noticed in the zooming areas that GAP-TV
is severely blurry, PnP-FFDnet(GAP) and PnP-FFDnet(ADMM) is kind of over
smooth around the edges. Our model can achieve cleaner results with sharper
edges.

##  4 Experiment

###  4.1 Experiment Setting

####  4.1.1 Architecture Specifics

For our learned network, we have experimented with various network
architectures. Specifically, for the DE-RNN model, we adopt the architecture
from BIRNAT [7]. Regarding its two-stage (forward+backward) RNN as a whole, we
iteratively feed the output of the backward RNN back as the input of the
forward one. For the DE-GAP model, we employ different neural networks as
denoisers $\mathcal{D}_{\theta}$ and utilize the real spectral norm [31] for
convergence purposes. We found that some architectures can yield fairly good
performance while combining our proposed DEQ for SCI. In a summary, these
feasible network architectures are Unet [20] with real spectral norm (denoted
as RSN-Unet), Unet with 3D convolutional kernels (denoted as Unet-3D), simple
CNN networks without and with real spectral norm (denoted as CNN and RSN-CNN,
respectively), and FFDnet [36].

####  4.1.2 Training Details

Following BIRNAT [7], we choose the dataset DAVIS2017 [18] for training.
DAVIS2017 has $90$ different scenes and in total $6208$ frames. We crop its
video frames to video patch cubes with the spatial size of $256\times
256\times 8$, and obtain $26,000$ training samples with data augmentation.
Then we train the neural network for $30$ epochs, and initialize the learning
rate as $1\times 10^{-3}$, which would be reduced with a decay rate $10\%$
every $10$ epochs. During training, we utilize Anderson acceleration for both
the forward and backward pass fixed-point iterations.

Figure 6: Comparison of selected reconstruction results of real data Water
Balloon with the spatial size of $512\times 512\times 10$. Reconstruction of
the real data is more difficult than simulations due to the inevitable
measurement noise. As shown in this figure, GAP-TV, DeSCI, and PnP-FFDnet
(GAP) have more artifacts and distortions around margins. Our model can
maintain a clear and accurate image structure, thus leading to higher
performance.

###  4.2 Experiment Results

####  4.2.1 Comparisons on Datasets

For evaluations, we test our proposed DE-RNN and DE-GAP on six classical
simulation datasets including Kobe, Runner, Drop, Traffic, Vehicle, and Aerial
[34] with the spatial size of ${256\times 256}$ and compression ratio $B$=8.
The quantitative comparison results with other video SCI reconstruction
algorithms including GAP-net [16], GAP-TV [32], E2E-CNN [19] and PnP-FFDnet
[34] on Peak Signal to Noise Ratio (PSNR) and structured similarity (SSIM)
[25] are provided in Table 1. What stands out in the table is that our method
achieves around $0.1$ dB improvement in PSNR and $0.4$ in SSIM in comparison
to others. The improvement of SSIM indicates our method could reconstruct
images with relative fine structure, which is confirmed by qualitative
evaluations in Fig. 5. Specifically, we observe that: i) GAP-TV results have
obvious ghosts and fail in high-quality structure reconstruction. For
instance, the cars in the Traffic scene reconstructed by GAP-TV are all with
heavy blur. ii) in comparison to GAP-TV, our method reconstruct explicit
content. iii) PnP-FFDnet approaches often cause distortion around margins,
while our results maintain a clear and accurate structure.

To sum up, both the quantitative and qualitative comparisons demonstrate that
our method could achieve competitive performance in contrast to other
algorithms. We do notice that there are some recent work using complicated
deep networks to obtain better results than ours [6, 28, 26, 16]. However,
these handcraft designs of different network structures are not necessarily
converging to a stable point. By contrast, our paper aims to provide a stable
solution for SCI reconstruction.

Recalling Fig. 2, where we have run existing methods and our algorithm for
iterations, while RNN and PnP fail in stable recovery, our method could
converge to a fixed point and maintain at high-level results. Reconstructed
frames in Fig. 1 further verified this virtue of our proposed algorithm.

Figure 7: Comparison of selected reconstruction results of real data Chopper
Wheel with the spatial size of $256\times 256\times 14$. Reconstruction of the
real data is more difficult than simulations due to the inevitable measurement
noise. As shown in this figure, GAP-TV, DeSCI, and PnP-FFDnet (GAP) have more
ghosts in the areas with large motion. Our model can lead to higher
performance.

####  4.2.2 Real-world Data Reconstruction

We also evaluate the DE-GAP model on real-world dataset Water Balloon and
Chopper Wheel with the spatial size of ${512\times 512\times 10}$ [19] and
${256\times 256\times 14}$ [15] captured by real video SCI cameras. Note that
this is more challenging due to the unavoidable noise inside the real
measurements, which demands the high robustness of the algorithm.

We compare the results with other algorithms including GAP-TV [32], DeSCI [14]
and PnP-FFDnet [34], as shown in Figs. 7 and 6. The reconstruction results on
real-world data demonstrate the effectiveness and generalization of our
proposed method. Note that the reconstruction results of real data are
achieved by the model trained to utilize the simulation mask, which means that
our proposed model is kind of flexible and can achieve stable results by the
virtue of the fact that our model can be theoretically infinitely extended.
Specifically, we observe that: i) GAP-TV and DeSCI often generate a lot of
artifacts and show noisy texture. ii) PnP-FFDnet has artifacts and distortions
around margins. iii) In contrast to them, our method shows high-quality
results with clear content and structure.

####  4.2.3 Processing Time

Though we equivalently realize infinite optimization iterations with deep
neural networks plugged in to perform video SCI reconstruction, our designed
methods elegantly avoid long inference time. As Table 2 shows, our method only
needs a short processing time in comparison to other algorithms. The source
code of our algorithm will be made available to the public to be used for
other tasks of inverse problems.

Table 2:  Average running time per measurement in seconds by different
algorithms on classical six datasets for video SCI reconstruction. While
permitting unlimited iterative steps and infinite network depth, our method
needs shorter inference time in contrast to other algorithms.

GAP-TV | DeSCI | PnP-FFDnet | DE-RNN | DE-GAP
---|---|---|---|---
4.2 | 6180 | 3.0 | 4.68 | 1.90

##  5 Conclusion

In this paper, to solve the problems of memory requirement and unstable
recovery in existing methods, we propose deep equilibrium models for video
SCI. Fusing data-driven regularization and stable convergence in a
theoretically sound manner, we combine DEQ with existing methods and design
two novel models, i.e., DE-RNN and DE-GAP. Each equilibrium model implicitly
learns a nonexpansive operator by training the embedded neural network and
analytically computes the fixed point, thus enabling unlimited iterative steps
and infinite network depth with only a constant memory requirement in the
training and inference process. Furthermore, we derive convergence theory for
each equilibrium model to ensure the results of our models converge to
equilibrium. We evaluate our proposed models using different neural networks
as the implicit operator on a variety of simulations and real datasets. In
comprehensive comparisons to existing algorithms, both quantitative and
qualitative evaluations of our results demonstrate the effectiveness and
stability of our proposed method.

## References

  * [1] Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  Trellis networks for sequence modeling.  arXiv preprint arXiv:1810.06682, 2018.
  * [2] Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  Deep equilibrium models.  arXiv preprint arXiv:1909.01377, 2019.
  * [3] Shaojie Bai, Vladlen Koltun, and J. Zico Kolter.  Multiscale deep equilibrium models.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 5238–5250. Curran Associates, Inc., 2020.
  * [4] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends in Machine Learning, 3(1):1–122, January 2011.
  * [5] Stanley H. Chan, Xiran Wang, and Omar A. Elgendy.  Plug-and-play ADMM for image restoration: Fixed-point convergence and applications.  IEEE Transactions on Computational Imaging, 3:84–98, 2017.
  * [6] Ziheng Cheng, Bo Chen, Guanliang Liu, Hao Zhang, Ruiying Lu, Zhengjue Wang, and Xin Yuan.  Memory-efficient network for large-scale video compressive sensing.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 16246–16255, 2021.
  * [7] Ziheng Cheng, Ruiying Lu, Zhengjue Wang, Hao Zhang, Bo Chen, Ziyi Meng, and Xin Yuan.  Birnat: Bidirectional recurrent neural networks with adversarial training for video snapshot compressive imaging.  In European Conference on Computer Vision, pages 258–275. Springer, 2020.
  * [8] Raj Dabre and Atsushi Fujita.  Recurrent stacking of layers for compact neural machine translation models.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 6292–6299, 2019.
  * [9] Mostafa Dehghani, Stephan Gouws, Oriol Vinyals, Jakob Uszkoreit, and Łukasz Kaiser.  Universal transformers.  arXiv preprint arXiv:1807.03819, 2018.
  * [10] Davis Gilton, Gregory Ongie, and Rebecca Willett.  Deep equilibrium architectures for inverse problems in imaging.  arXiv preprint arXiv:2102.07944, 2021.
  * [11] Aidan N Gomez, Mengye Ren, Raquel Urtasun, and Roger B Grosse.  The reversible residual network: Backpropagation without storing activations.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 2211–2221, 2017.
  * [12] Shirin Jalali and Xin Yuan.  Snapshot compressed sensing: Performance bounds and algorithms.  IEEE Transactions on Information Theory, 65(12):8005–8024, 2019\.
  * [13] X. Liao, H. Li, and L. Carin.  Generalized alternating projection for weighted-$\ell_{2,1}$ minimization with applications to model-based compressive sensing.  SIAM Journal on Imaging Sciences, 7(2):797–823, 2014.
  * [14] Yang Liu, Xin Yuan, Jinli Suo, David J Brady, and Qionghai Dai.  Rank minimization for snapshot compressive imaging.  IEEE transactions on pattern analysis and machine intelligence, 41(12):2990–3006, 2018.
  * [15] Patrick Llull, Xuejun Liao, Xin Yuan, Jianbo Yang, David Kittle, Lawrence Carin, Guillermo Sapiro, and David J Brady.  Coded aperture compressive temporal imaging.  Optics express, 21(9):10526–10545, 2013.
  * [16] Ziyi Meng, Shirin Jalali, and Xin Yuan.  Gap-net for snapshot compressive imaging.  arXiv preprint arXiv:2012.08364, 2020.
  * [17] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32:8026–8037, 2019.
  * [18] Jordi Pont-Tuset, Federico Perazzi, Sergi Caelles, Pablo Arbeláez, Alex Sorkine-Hornung, and Luc Van Gool.  The 2017 davis challenge on video object segmentation.  arXiv preprint arXiv:1704.00675, 2017.
  * [19] Mu Qiao, Ziyi Meng, Jiawei Ma, and Xin Yuan.  Deep learning for video compressive sensing.  APL Photonics, 5(3):030801, 2020.
  * [20] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [21] Ernest Ryu, Jialin Liu, Sicheng Wang, Xiaohan Chen, Zhangyang Wang, and Wotao Yin.  Plug-and-play methods provably converge with properly trained denoisers.  In International Conference on Machine Learning, pages 5546–5557. PMLR, 2019.
  * [22] S. Sreehari, S. V. Venkatakrishnan, B. Wohlberg, G. T. Buzzard, L. F. Drummy, J. P. Simmons, and C. A. Bouman.  Plug-and-play priors for bright field electron tomography and sparse interpolation.  IEEE Transactions on Computational Imaging, 2(4):408–423, 2016\.
  * [23] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg.  Plug-and-play priors for model based reconstruction.  In 2013 IEEE Global Conference on Signal and Information Processing, pages 945–948, 2013.
  * [24] Homer F Walker and Peng Ni.  Anderson acceleration for fixed-point iterations.  SIAM Journal on Numerical Analysis, 49(4):1715–1735, 2011.
  * [25] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli.  Image quality assessment: from error visibility to structural similarity.  IEEE transactions on image processing, 13(4):600–612, 2004.
  * [26] Zhengjue Wang, Hao Zhang, Ziheng Cheng, Bo Chen, and Xin Yuan.  Metasci: Scalable and adaptive reconstruction for video compressive sensing.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2083–2092, 2021.
  * [27] Kaixuan Wei, Angelica Aviles-Rivero, Jingwei Liang, Ying Fu, Carola-Bibiane Schönlieb, and Hua Huang.  Tuning-free plug-and-play proximal algorithm for inverse imaging problems.  In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 10158–10169. PMLR, 13–18 Jul 2020.
  * [28] Zhuoyuan Wu, Jian Zhang, and Chong Mou.  Dense deep unfolding network with 3d-cnn prior for snapshot compressive imaging.  arXiv preprint arXiv:2109.06548, 2021.
  * [29] J. Yang, X. Liao, X. Yuan, P. Llull, D. J. Brady, G. Sapiro, and L. Carin.  Compressive sensing by learning a Gaussian mixture model from measurements.  IEEE Transaction on Image Processing, 24(1):106–119, January 2015\.
  * [30] J. Yang, X. Yuan, X. Liao, P. Llull, G. Sapiro, D. J. Brady, and L. Carin.  Video compressive sensing using Gaussian mixture models.  IEEE Transaction on Image Processing, 23(11):4863–4878, November 2014.
  * [31] Yuichi Yoshida and Takeru Miyato.  Spectral norm regularization for improving the generalizability of deep learning.  arXiv preprint arXiv:1705.10941, 2017.
  * [32] Xin Yuan.  Generalized alternating projection based total variation minimization for compressive sensing.  In 2016 IEEE International Conference on Image Processing (ICIP), pages 2539–2543. IEEE, 2016.
  * [33] X. Yuan, D. J. Brady, and A. K. Katsaggelos.  Snapshot compressive imaging: Theory, algorithms, and applications.  IEEE Signal Processing Magazine, 38(2):65–88, 2021.
  * [34] X. Yuan, Y. Liu, J. Suo, and Q. Dai.  Plug-and-play algorithms for large-scale snapshot compressive imaging.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [35] Xin Yuan, Yang Liu, Jinli Suo, Fredo Durand, and Qionghai Dai.  Plug-and-play algorithms for video snapshot compressive imaging.  IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 1–1, 2021.
  * [36] Kai Zhang, Wangmeng Zuo, and Lei Zhang.  Ffdnet: Toward a fast and flexible solution for cnn-based image denoising.  IEEE Transactions on Image Processing, 27(9):4608–4622, 2018.


# Multiscale approach for the network compression-friendly ordering

Ilya Safro Mathematics and Computer Science Division, Argonne National
Laboratory<EMAIL_ADDRESS>Boris Temkin Faculty of Mathematics and Computer
Science, The Weizmann Institute of Science<EMAIL_ADDRESS>
###### Abstract

We present a fast multiscale approach for the network minimum logarithmic
arrangement problem. This type of arrangement plays an important role in a
network compression and fast node/link access operations. The algorithm is of
linear complexity and exhibits good scalability which makes it practical and
attractive for using on large-scale instances. Its effectiveness is
demonstrated on a large set of real-life networks. These networks with
corresponding best-known minimization results are suggested as an open
benchmark for a research community to evaluate new methods for this problem.

Keywords: multiscale algorithms, minimum logarithmic arrangement, network
compression

##  1 Introduction

Finding a suitable compressed representation of large-scale networks is
intensively studied in both practical and theoretical branches of data mining
[12, 2, 11, 5, 26]. In particular, the success of applying some of the
recently proposed compression schemes [11, 5, 3] strongly depends on the
“compression-friendly” arrangement of network nodes. Usually, the goal of
these arrangements is to order the nodes such that the endpoints of network
links (edges) are located as close as possible. Doing so leads to a more
compact representation of links and allows a better performance of compression
schemes and network element access operations.

In [11], Chierichetti et al. propose a combinatorial optimization problem,
namely, the minimum logarithmic arrangement problem, that seeks a nearly
optimal information-theoretical compressed encoding size for all network
links. This is achieved by ordering the network nodes and assigning to them
unique integer values (ids) such that the endpoints of a link will obtain
close values. The problem has been proven to be NP-hard, and two heuristics
for ordering the social networks were given. In this paper, we present a
multiscale method for approximating a generalized link-weighted and node-
volumed version of the minimum logarithmic arrangement for general networks.
The importance of a link-weighted property (when each link is assigned by a
nonnegative weight) for this problem is twofold. First, the link weight can
measure the significance of that link. For example, in many cases we have
information regarding its access frequency, and we would prefer that
frequently accessed links would be compressed better. Second, the multiscale
algorithmic framework admits a natural aggregation of weighted links at
different scales of the problem representation (discussed later). Similary to
the link weights, all nodes are assigned by nonnegative volumes that represent
the length of their segments captured in the ordering of the network nodes.

Our multiscale algorithm is based on the algebraic multigrid (AMG) [8]
methodology for linear ordering problems [14, 23]. The main objective of the
AMG-based framework is to construct a hierarchy of problems (coarsening), each
approximating the original problem, but with fewer degrees of freedom. This is
achieved by introducing a sequence of successive projections of networks’
graph Laplacians into lower-dimensional spaces and solving the problem in
them. The multiscale framework has two key advantages that make it attractive
for applying on modern large-scale instances: it exhibits a linear complexity,
and it can be relatively easily parallelized and implemented by using standard
matrix-vector operations. Another advantage of the multiscale framework is its
heterogeneity, expressed in the ability to incorporate external appropriate
optimization algorithms (as a refinement) in the framework at different
scales. For more detailed surveys on the multiscale methods and their
parallelization for combinatorial optimization problems, we refer the reader
to [9, 27, 7].

In contrast to [11], the main goal of this work is to provide a generic solver
for the minimum logarithmic arrangement for general networks. Social networks
(including big parts of Web graph) commonly consist of a combination of
structural properties (such as degree distribution, small diameter and
expander-like topology) that enable fast greedy methods (usually based on some
preferential ordering of graph traversal) to find arrangements for further
high quality compression. These methods can be successful because most of the
links exhibit good locality [18, 28] and only a small number of them can be
considered as global edges that connect well-separated nodes/regions of a
network. Thus, these networks can hardly be considered as strongly irregular
and difficult instances. In real life, however, there occur many situations
when the structure of a stored network (or a part of it) is complex and
irregular (such as decision/detailed supply/infection spread networks and even
parts of social networks that do not exhibit power law degree distribution),
which creates many contradictions between greedy local decisions and solutions
that consider a more global picture. We address this type of problem by
introducing a multiscale solver.

The experimental part of this work shows how far are the existing state-of-
the-art ordering heuristics from being optimal. We demonstrate significant
improvement for minimization of logarithmic arrangement for various families
of networks (including social networks and other (ir)regular instances). The
comparison was performed with several recently introduced methods in [1, 5, 6,
11]. In almost all cases, our solver exhibited better numerical results than
previous best-known results.

We introduce notation and necessary definitions in Section 2. The main
algorithm is described in Section 3. The computational results and discussion
about different parameter settings are presented in Section 4. In Section
Section 5 we conclude and provide possible future research directions. For
readers interested in their own implementation, we note that some details
regarding the multiscale scheme (that are not related strongly to the
discussed problem) are omitted in this paper; they can be found in [22, 23,
20, 9].

##  2 Notation and problem definition

A network is described by a weighted directed graph $G=(V,~{}E)$, where
$V=\\{1,2,...,n\\}$ is the set of nodes (vertices) and $E$ is the set of
directed edges. If $ij\in E$, then there exists an edge $i\rightarrow j$.
Denote by $w_{ij}$ the non-negative weight of the directed edge $ij$ between
nodes $i$ and $j$; if $ij\notin E$, then $w_{ij}=0$. Let $\pi$ be a bijection

$\pi:~{}V~{}\longrightarrow~{}\\{1,2,...,n\\}~{}~{}~{}.$

The purpose of the link-weighted version of the minimum logarithmic
arrangement problem (MLogA) is to minimize

$\sum_{ij\in E}w_{ij}\lg|\pi(i)-\pi(j)|$ (1)

over all possible permutations $\pi$. The base of a logarithm ($\lg$) will be
always 2. We define the generalized form of this problem (GMLogA) that emerges
during the multilevel solver; each vertex $i$ is assigned with a $volume$ (or
$length$), denoted $v_{i}$. Given the vector of all volumes, $v$, the task now
is to minimize the cost

$\textsf{c}(G,x)=\sum_{ij\in E}w_{ij}\lg|x_{i}-x_{j}|$ (2)

over all possible $\pi$, where
$x_{i}=\frac{v_{i}}{2}+\sum_{k,\pi(k)<\pi(i)}v_{k}$; that is, each vertex is
positioned at its center of mass, capturing a segment on the real axis that
equals its length. The original form of the problem is the special case where
all the volumes are equal.

We denote the process of creating a weighted undirected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ by
$\mathcal{G}\leftarrow\textsf{un}(G)$, where $G$ is a given directed graph. By
applying $\textsf{un}(G)$, we obtain $\mathcal{V}\leftarrow V$ and

$ij\in\mathcal{E}\Leftrightarrow ij\in E\text{ or }ji\in E~{}~{}~{}.$

The edge weights are accumulated from the set of directed edges

$\forall~{}ij\in\mathcal{E}~{}w_{ij}=\sum_{ij\in E}w_{ij}+\sum_{ji\in
E}w_{ji}$

allowing use of multigraphs. Both MLogA and GMLogA are formulated for the
undirected graphs similarly to their directed versions. In this paper we will
find approximate solutions for the following problem:

###### Problem 2.1

Given a directed weighted network graph $G$, the objective is to find a
permutation $\pi$ that minimizes $\textsf{c}(G,x_{\pi})$, where $x_{\pi}$ is a
vector of the respective node coordinates restricted by $\pi$.

Note that the problem is formulated for both directed and undirected graphs.
We define111$\beta$ with no subindexes will be used where appropriate.
$\beta_{G,x_{\pi}}$ as the main unit for the empirical comparison in Section
4:

$\beta_{G,x_{\pi}}=\textsf{c}(G,x_{\pi})/\sum_{ij\in E}w_{ij}~{}.$ (3)

##  3 The algorithm

In the multiscale framework we construct a hierarchy of decreasing-size
graphs, $G_{0},G_{1},...,G_{k}$, by coarsening, starting from the given graph
$G_{0}\leftarrow\textsf{un}(G)$. In other words, we create an undirected
version of a given directed graph by assigning to the undirected edges bigger
weights if, in the original unweighted graph, the directed edges had been
reciprocal. At the coarsest level Problem 2.1 is solved exactly, and starting
from the $(k-1)$th level it is formulated and approximated by successive
prolongation of the solution from the previous coarser level. This entire
process is called a $V$-$cycle$ (see Figure 1).

Figure 1: Scheme of V-cycle. The arrows show the order in which smaller graphs
are created and revised for approximation.

###  3.1 Coarsening

In the present work, coarsening is interpreted as a modified process of
weighted aggregation reinforced by the algebraic distance couplings for
logarithmic sum minimization. For a detailed discussion about the algebraic
distance-based weighted aggregation and multiscale graph organization related
to the general graph optimization problems, we refer the reader to [20]. We
will briefly repeat its basic components for the completeness of the paper and
will concentrate on the modifications related to Problem 2.1.

####  3.1.1 Algebraic distance coupling

Algebraic distance-based coupling is a measure of connectivity strength
between two nodes connected by an edge [20, 10]. Given the Laplacian of a
graph, denoted by $L=D-W$, where $W$ is a weighted adjacency matrix of a graph
and $D$ is the diagonal matrix with entries $d_{ii}=\sum_{j}w_{ij}$, we define
an iteration matrix $H$ for Jacobi over-relaxation (JOR) as

$H=(D/\omega)^{-1}((1/\omega-1)D+W_{l}+W_{u})~{},$

where $0\leq\omega\leq 1$ (see Appendix A) and $W_{l}$ and $W_{u}$ are the
strict lower and upper triangular parts of $W$, respectively.

###### Definition 3.1

The algebraic distance coupling $\rho_{ij}$ is defined as

$\rho_{ij}=1/\big{(}\sum_{r=1}^{R}\lg|\chi_{i}^{(k,r)}-\chi_{j}^{(k,r)}|\big{)}~{},$

where $\chi^{(k,r)}=H^{k}\chi^{(0,r)}$ is a relaxed randomly initialized test
vector (i.e., $\chi^{(0,r)}$ is a random vector sampled over [-1/2, 1/2]), R
is a number of initial test vectors, and k is a number of iterations.

Note that this definition is a modified version of the original definition
from [20, 10] which makes it more suitable for Problem 2.1. Several
interesting properties (including convergence and model description) still can
be proved similarly to [10].

####  3.1.2 Coarse graph construction

Considering $\rho_{ij}$ as an edge strength measure, one can construct a
coarse graph by defining a classical AMG interpolation matrix $P$ that will
project the fine graph to the coarse graph (i.e., to the lower-dimensional
space). The projection is represented as

$L_{c}\leftarrow PL_{f}P^{T}~{},$

where $L_{f}$ and $L_{c}$ are the Laplacians of fine ($G_{f}$) and coarse
($G_{c}$) graphs, respectively.

We begin by selecting a set of seed nodes $C\subset V_{f}$ that will represent
the centers of future coarse nodes (or aggregates). In fact, $C$ is
interpreted as a dominating set of $V_{f}$ (not necessarily the minimum size)
such that other (fine) nodes $F=V_{f}\setminus C$ should be strongly coupled
to $C$. This can be done by traversing all nodes and identifying for every
visited node $i$ whether

$\sum_{j\in C}\rho_{ij}/\sum_{j\in
V_{f}}\rho_{ij}\geq\Theta_{1}~{}~{}~{}\text{ and }~{}~{}~{}\sum_{j\in
C}w_{ij}/\sum_{j\in V_{f}}w_{ij}\geq\Theta_{2}~{},$ (4)

where $\Theta_{(1,2)}$ are the parameters of coupling strength (see Appendix
A). The order in which the nodes are traversed is based on the future volume
principle [20], which is a measure of how large (representative in terms of
the current minimization problem) a coarse node can be. To keep the linear
complexity of the entire framework, the order need not be calculated exactly
but only roughly (for example using bucketing sort). In contrast to the
multiscale approach for the generalized minimum $p$-sum problem [20], we found
that automatically moving to $C$ those nodes that have exceptionally large
future volume is not necessary and even has a small negative impact in several
social networks.

After identifying $C$, we define for each $i\in F$ its coarse neighborhood
$N^{c}_{i}$ that contains a limited set of $C$-nodes to which $i$ is
connected. The criterion for choosing $C$-nodes to $N^{c}_{i}$ is also based
on $\rho_{ij}$. Let $I(j)$ be the ordinal number in the coarse graph of the
node that represents the aggregate around a seed whose ordinal number at the
fine level is $j$. The classical AMG interpolation matrix $P$ is defined by

$P_{iI(j)}~{}=~{}\left\\{\begin{array}[]{lll}w_{ij}/\sum\limits_{k\in
N^{c}_{i}}w_{ik}&\textsf{for }i\in F,~{}j\in N^{c}_{i}\\\ 1&\textsf{for }i\in
C,~{}j=i\\\ 0&\textsf{otherwise}~{}~{}~{}.\end{array}\right.$ (5)

$P_{iI(j)}$ thus represents the likelihood of $i$ belonging to the $I(j)$th
aggregate. The edge connecting two coarse aggregates $p$ and $q$ is assigned
with the weight $w_{pq}=\sum_{k\not=l}P_{kp}w_{kl}P_{lq}$. The volume of the
$p$th coarse aggregate is $\sum_{j}v_{j}P_{jp}$.

###  3.2 Uncoarsening

The uncoarsening process starts by solving Problem 2.1 at the coarsest level.
Since the number of nodes at the coarsest level is very small (in our tests it
was 9), the problem can be solved exactly by exhaustive search.

####  3.2.1 Minimizing the contribution of one node

Before proceeding to the stages of coarse-to-fine projection of a coarse-level
solution, we describe how to approximate a solution of Problem 2.1 for a
single node only. This will be a basic step in the initialization and
relaxation described in Sections 3.2.2 and 3.2.3.

Denote by $N_{i}$ the set of $i$th neighbors with already assigned coordinates
$\tilde{x}_{j}$. To minimize the local contribution of $i$ to the total energy
(2), we assign to it a coordinate $x_{i}$ that minimizes

$\sum_{j\in N_{i}}w_{ij}\lg|x_{i}-\tilde{x}_{j}|~{}.$ (6)

Since for every $j\in N_{i}$, $x_{i}=\tilde{x}_{j}$ implies that the sum (6)
is minus infinity, we resolve this ambiguity by setting

$x_{i}=\tilde{x}_{t}\iff t=\arg\min_{k\in N_{i}}\sum_{k\neq j\in
N_{i}}w_{kj}\lg|\tilde{x}_{k}-\tilde{x}_{j}|~{}.$ (7)

The trivial exact solution has a running time $O(|N_{i}|^{2})$, as it requires
to compute $|N_{i}|$ sums, each one with $|N_{i}|$ terms. Thus, to preserve
the linear complexity of the entire algorithm, one can use the trivial
solution for nodes with small $|N_{i}|$ only. We will approximate (7) using
the heuristic that seeks the nearly minimum sum in the point of maximal
density.

Consider set $\\{\tilde{x}_{j}:~{}j\in N_{i}\\}$ as independent and
identically distributed samples of a random variable with unknown distribution
and $w_{ij}$ as the posteriori probability of a sample $\tilde{x}_{j}$.
Assuming that, we have to choose a point where the estimated probability
density is maximized. Various approaches exist for density estimation [25].
One of the most popular is called the “Parzen window” (or kernel density
estimation) method. In this method, the density at point $x$ is estimated as

$d(x)=\frac{1}{|N_{i}|h}\sum_{j\in
N_{i}}K\left(\frac{|x-\tilde{x}_{j}|}{h}\right)~{},$ (8)

where $K$ is a kernel and $h$ is a smoothing parameter called the bandwidth.
To simplify calculations, we choose a kernel similar to the Gaussian kernel.
Thus, the estimated density at point $x$ is

$\hat{d}(x)=\frac{1}{|N_{i}|h}\sum_{j\in
N_{i}}w_{ij}2^{|x-\tilde{x}_{j}|/h}~{}.$

Note that the factor $\frac{1}{|N_{i}|h}$ may be omitted without changing the
maximal’s location. Moreover, it is sufficient to calculate
$\hat{d}(\tilde{x}_{j})$ for every $j\in N_{i}$. This approach allows us to
compute the density in all $|N_{i}|$ points in linear time using following
algorithm.

Input: $x_{1},\ldots,x_{k}$ \- samples in increasing order

Input: $h$ \- smoothing parameter (window)

Input: $p_{1},\ldots,p_{k}$ \- posteriory probabilities

Output: $\hat{d}_{1},\ldots,\hat{d}_{k}$ \- estimated density

/* calculate incrementally influence of left neighbors */

$s_{1}\leftarrow p_{1}$

for _$t=2$ to $k$_ do $s_{t}=p_{t}+s_{t-1}2^{(x_{t-1}-x_{t})/h}$ end /*
calculate incrementally influence of right neighbors */

$r_{k}\leftarrow p_{k}$

for _$t=k-1$ to $1$_ do $r_{t}=p_{t}+s_{t+1}2^{(x_{t}-x_{t+1})/h}$ end

/* aggregate results */

for _$t=1$ to $k$_ do $\hat{d}_{t}=s_{t}+r_{t}-p_{t}$ end

Algorithm 1 Linear algorithm for density estimation

Likewise, for all density estimations techniques the crucial step is to choose
the bandwidth $h$. On the one hand, it should be big enough to smooth the
peaks and highly-oscillatory components. On the other hand, if $h$ is too big,
then the maximum density will always be located at the middle. Based on
experiments, we choose $h=N/2\lg(N)$, where
$N=\max(\tilde{x}_{i})-\min(\tilde{x}_{i})$.

To justify the chosen heuristic, we present our results from experiments on
real data. In Figure 2 we show the error distribution for our approximation.
The error is calculated as

$\text{error}=\hat{\theta}-\theta^{*}~{},$

where $\theta^{*}$ is the energy (excluding $-\infty$ terms) at the exact
solution of (7) and $\hat{\theta}$ is an approximation. We can easily see that
most of the time the approximation error is either zero or very small.

Figure 2: Error distribution for density approximation on real-life networks.
Each curve corresponds to one network. Each point on a curve corresponds to
one difference between the exact (quadratic) solution of 7 and its
approximation. The comparisons (points) are ordered by their errors in an
ascending order.

####  3.2.2 Initialization

Given is the arrangement of the coarse-level aggregates in its generalized
form, where the center of mass of each aggregate $j\in C$ is positioned at
$x_{I(j)}$ along the real axis. We initialize the fine-level arrangement by
letting each seed $j\in C$ inherit the position of its respective aggregate:
$x_{j}\leftarrow x_{I(j)}$. However, in contrast to the initialization by
stages (that was done for the minimum $p$-sum and the minimum workbound
problems [23]) in which the $F$-nodes’ positions are calculated according to
all already initialized fine-level neighbors, we found that the initialization
that is inspired by the principles of the classical AMG coarse solution
projection $x^{f}~{}\leftarrow~{}P^{T}x^{c}$ can produce significantly better
approximations. In other words, after initializing the fine-level $C$-nodes,
$F$-nodes will be placed one by one using Algorithm 1, which takes into
account only the $C$-neighbors of a current $F$-point. The entire
initialization scheme is presented in Algorithm 2.

Input: arrangement of a coarse graph

Output: initial ordering of fine level nodes

/* Initialize $C$-points */

for _$\text{all }j\in C$_ do $x_{j}\leftarrow x_{I(j)}$ end /* Initialize
$F$-points */

for _$\text{all }j\in F$_ do $x_{j}\leftarrow$ find best coordinate for $j$
using Algoriothm 1 with $N_{j}$ end /* Legalization of the coordinates */

for _$\text{all }j\in V$_ do $x_{i}=v_{i}/2+\sum_{x_{k}<x_{i}}v_{k}$ end

Algorithm 2 Initialization of the fine level

####  3.2.3 Relaxations and strict minimization

The initialization produces a feasible (that satisfies the constraint of (2))
solution inherited from the coarse level. This solution is further improved by
employing a small number of relaxation sweeps of two types: compatible and
Gauss-Seidel (GS). The goal of the compatible relaxation is to improve the
positions of $F$-nodes while keeping the $C$-nodes invariant. One sweep of the
compatible relaxation can be described by the two last “for” cycles of
Algorithm 2 when in the first cycle the best coordinate is calculated over all
neighbors.

Having the improved ordering of $F$-nodes relative to the $C$-nodes, we apply
GS relaxation. This relaxation traverses all nodes one by one and tries to
find a best position for every node such that its own contribution to the
total sum of logarithms will be minimized. In all our experiments we observed
that, after improving the order with the compatible relaxation, the GS
relaxation has almost no influence on the order. This is in contrast to other
multiscale linear ordering methods in which the GS relaxation was one of the
most powerful and crucial components.

####  3.2.4 Strict minimization

We applied two methods of strict minimization (SM, also known as “local
refinement”): node-by-node minimization (N-N) and window minimization (WM).
The principal difference between relaxations and SM is that in SM each change
is accepted if and only if it reduces the entire sum (2) and not the local
node contribution only (6).

In the N-N minimization all nodes are scanned according to their current
order, and each vertex $i$, in its turn, is checked for the best position over
some small enough segment that $i$ belongs to. The $k$ left and $k$ right
candidate positions are scanned, and the one with the minimal cost is chosen.

In the WM (see [23]), the total cost of the arrangement is reduced by
numerical minimization of a collective contribution of a small group of
consecutive nodes. Given a current feasible solution $\tilde{x}$ of the graph
logarithmic arrangement, denote by $\delta_{i}$ a small correction to
$\tilde{x}_{i}$. Denote by
$\mathfrak{W}=\\{i_{1}=\pi^{-1}(s+1),...,i_{q}=\pi^{-1}(s+q)\\}$ a subset of
successive $q$ vertices in the current arrangement. The goal of the local
minimizaiton problem is then to find $\delta$ such that

$\sum_{\begin{subarray}{c}i,j\in\mathfrak{W}\end{subarray}}w_{ij}\lg|\tilde{x_{i}}+\delta_{i}-\tilde{x_{j}}-\delta_{j}|+\sum_{\begin{subarray}{c}i\in\mathfrak{W}\\\
j\not\in\mathfrak{W}\end{subarray}}w_{ij}\lg|\tilde{x_{i}}+\delta_{i}-\tilde{x_{j}}|$
(9)

is minimized. The solution to this optimization problem can be approximated
numerically by using techniques similar to 3.2.1.

The situation with strict minimization is similar to that with GS. In all our
experiments, we found no need to use the more complicated WM. The simple N-N
minimization with small distance parameter $k<5$ was enough to obtain the best
results. This is a major advantage because in most multiscale schemes the
significant part of the running time is spent on the uncoarsening. The lack of
need (or, more correctly, optionality) for GS and WM (or other collective
optimization) also strongly advocates the choice of the C-nodes because the
further compatible relaxation solves the problem well.

After having defined all the components of the V-cycle, we present the full
scheme in Algorithm 3.

Input: $L_{f}$ is the Laplacian of undirected graph $\mathcal{G}$ (initially
obtained from $\textsf{un}(G)$)

Output: approximated graph minimum logarithmic arrangement

if _$\mathcal{G}$ is small enough_ then solve the problem exactly

return the arrangement end /* Coarsening */

calculate algebraic distances $\rho_{ij},~{}\forall ij\in\mathcal{E}$

identify the set of $C$-points (see (4))

construct the fine-to-coarse projection $P$ (see (5))

$L_{c}\leftarrow PL_{f}P^{T}$

/* Recursive call */

$\Pi\leftarrow~{}$V-cycle($L_{c}$)

/* Uncoarsening */

$\pi_{1}\leftarrow~{}$initialization($\Pi$)

$\pi_{2}\leftarrow~{}$compatible-relaxation($\pi_{1}$)

$\pi_{3}\leftarrow~{}$GS-relaxation($\pi_{2}$)

$\pi_{4}\leftarrow~{}$N-N refinement($\pi_{3}$)

return $\pi_{4}$

Algorithm 3 ms-GMLogA: Full scheme of one V-cycle for GMLogA.

##  4 Computational results

The implementation of Algorithm 3 is based on the multiscale framework used in
[23]. The implementation is nonparallel and has not been fully optimized. The
results (arrangements and running times) should be considered only
qualitatively and can certainly be further improved by more advanced
implementation and multiscale techniques.

Because of the practical significance of MLogA, we designed an open site [24]
with the benchmark graphs and a set of numerical results. In this site we
present the results listed below and invite the scientists who work on this
problem to submit interesting new networks, their solutions, and improved
arrangements for the existing networks.

Our benchmark consists of 100 graphs of different nature and size (most of
them are taken from [13] and [17]). For these graphs we created their
undirected versions (by applying $\textsf{un}(\cdot)$) and evaluated our
algorithms on both the directed and undirected versions. We evaluated two
baseline algorithms:

  * •

MinLA+N-N: the multiscale solver that finds an approximation of MinLA [22]
followed by GMLogA-oriented N-N fast postprocessing (see Section 3.2.4 with
$k=10$);

  * •

ms-GMLogA: GMLogA solver described in Algorithm 3.

Clearly, MinLA and GMLogA can have different optimal orderings. However, we
observed that qualitative MinLA orderings reinforced by appropriate
postprocessing can also lead to good solutions of GMLogA. In spite of the fact
that most of the best solutions were obtained with ms-GMLogA and because the
observation was done on a benchmark of large-scale real-life graphs (and there
exists intensive research on MinLA [14, 16]), we decided to present these
results too. We also mention that the Fiedler vector-based solutions (even
those that try to achieve a good local minimum for MinLA) usually produce
significantly fewer qualitative solutions even with appropriate
postprocessing.

In our experiments we use a WebGraph framework [5], which provides a simple
way to manage large graphs. In particular, the last version of WebGraph
produces random, lexicographical [5], Gray [6], double shingle [11], and
LayeredLPA [1] orderings for the graphs. The important feature of these
orderings is that they are completely endogenous (i.e., determined by the
graph itself), contrary to natural ordering, that is, the ordering in which
the graph was created. It would seem that the random ordering always has to be
significantly worser, than the lexicographic and natural orderings. However,
this is not always true, mostly because of the ways the latter have been
produced. For example, in highly parallel systems with shared memory, the
parallelized graph traversal algorithms can produce arrangements that have
been dumped in parallel from different processors. Each processor produces an
arrangement of good locality for a small subset of nodes, but the entire
arrangement does not possess this property. In particular, in our benchmark we
found two graphs whose random ordering was better than the natural and
lexicographic ones.

We use the average number of bits per link $\beta$ (3) as a unit of comparison
for the main results presented in Figures 3(a-f) and 4(a-b). In both figures
the left and right columns correspond to the comparison results of directed
and undirected graphs, respectively. In Figures 3(a-f), each plot area
consists of two curves. Each point in bold curves represents a ratio between
$\textsf{c}(G)$ obtained by ms-GMLogA and
$\min(\textsf{c}(G,\pi_{nat}),\textsf{c}(G,\pi_{lex}),\textsf{c}(G,\pi_{rnd}))$
for a particular graph from the benchmark, where $\pi_{nat}$, $\pi_{lex}$, and
$\pi_{rnd}$ are the best orderings obtained by natural, lexicographic, and
randomized arrangements, respectively. Similarly, each point in regular curves
corresponds to the ratio with MinLA+N-N in the nominator. To demonstrate the
robustness of the algorithm, we show several sets of its parameters. Each row
pair of subfigures corresponds to the same set of parameters for directed and
undirected graphs, respectively.

In the multiscale algorithms, there are several possible sources (see [9]) of
potentially higher than linear complexity (or linear with coefficients that
are too big). Here we address the four most important: (a) the number of
iterations in the relaxations/refinements, (b) the complexity of one iteration
of relaxation/refinement, (c) the order of interpolation that can increase the
complexity of the coarse graphs, and (d) too many C-nodes (addressed in
Appendix A). The first set of parameters (Figures 3(a,b)) (which is a
suggested default set) contains a mild configuration. The order of
interpolation (number of nonzero entries in one row of $P$) is only 1; $k$ in
N-N refinement is 5; and the number of relaxations of any kind is at most 20.
The calculation of the algebraic distance couplings in this set is based on 5
randomly initialized vectors and only 20 iterations of JOR. For almost all
graphs, we observed that ms-MLogA produces better orderings than does
MinLA+N-N. On average, ms-MLogA improves the graphs by more than 40% of their
ordering cost in comparison to natural, lexicographic, and random
arrangements.

We note that the most beneficial graphs (those that have ratios from 0.2 to
0.4) in both directed and undirected versions come from different collections
such as VLSI design networks, part of Web network, road maps, and Amazon
links. In general, these graphs differ significantly in their structural
properties. We observed that no particular part of the benchmark was less
beneficial than another, thus attesting to the generality of the proposed
method.

The difference between improvements obtained with different parameter sets is
not significant; however, we provide them to demonstrate the robustness of the
method. In the second set of parameters (Figure 3(c,d)) the number of random
initial vectors in the algebraic distance coupling is reduced to 1 only. It
leads to the faster coarsening (while keeping all convergence and model
properties of the algebraic distance [10]) and potentially weaker decisions
regarding edge strengths based on the algebraic distance. In the third set of
parameters (Figure 3(e,f)) we increase $k$ in N-N refinement from 5 to 10 and
the interpolation order from 1 to 2 at the finest level of hierarchy. At all
coarse levels these parameters are increased too, according to the logarithmic
increase scale described in [22].

In Figure 4(a,b) we show the difference between the very fast version of ms-
GMLogA with no refinement and fast relaxation (5 sweeps) and the slowest
version with very aggressive relaxation/refinement parameters ($k=25$ and 40
sweeps). Each subfigure contains two curves: regular and bold. The regular
corresponds to the fast version and the bold to the slowest one. It is easy to
see that the difference between these two versions is not significant and one
can certainly use a fastest version to obtain qualitative results.

The Gray [6] and double shingle [11] orderings have proven themselves as more
successful heuristics than lexicographic and natural orderings [24] and, thus,
we present their comparison with ms-GMLogA in separate Figure 5(a-b). Two
curves depicted in Figure 5(a) correspond to the comparison of Gray ordering
on directed (regular curve) and undirected (bold curve) graphs, respectively.
Each point on the curves corresponds to the ratio between $\beta$ values of
ms-GMLogA with default parameters and Gray ordering, respectively. Similarly,
the comparison of double shingle ordering [11] is shown in Figure 5(b). The
(double) shingle heuristic is based on the node similarity ordering derived
from estimation of Jaccard coefficients. In [11], the authors prove that
(double) shingle heuristic will be beneficial for networks that can be
described by the preferential attachment model.

Recently, the layered label propagation algorithm (LayeredLPA) has been
introduced in [1]. This propagation method is based on the Potts model [19].
This algorithm is significantly more successful than natural, random,
lexicographic, Gray and (double) shingle orderings. The success of the label
propagation methods has a similar nature to the success of the algebraic
distance coupling (see Section 3.1.1 and [10]) in which the propagation and
averaging of random values over the node neighborhoods is employed. However,
our multiscale method shows better results (see Figure 6). Note that the
LayeredLPA is introduced for undirected graphs only. We believe, that
introducing the AMG-based framework to the label propagation model can
significantly improve its quality.

(a) directed

(c) directed

(e) directed

(b) undirected

(d) undirected

(f) undirected

Figure 3: Main comparison results. Notations ’lg’ and ’la’ correspond to the
reults of ms-GMLogA and MinLA+N-N, respectively.

(a) directed

(b) undirected

Figure 4: Comparison of the fastest and the slowest versions of ms-GMLogA.
Notations ’no-N-N-lg’ and ’k25-N-N’lg’ correspond to the fast and slow
versions of ms-GMLogA, respectively.

(a) Gray ordering versus ms-GMLogA

(b) Double shingle versus ms-GMLogA

Figure 5: Comparison of Gray and double shingle orderings versus ms-GMLogA.
Figure 6: Comparison of LayeredLPA versus ms-GMLogA.

###  4.1 Scalability of ms-GMLogA

Good scalability is one of the most important advantages of the multiscale
algorithms. All components of our scheme are of linear complexity, and the
total complexity is $O(|V|+|E|)$. The dependence of the running time on the
graph size is depicted in Figure 7. Each small circle in the figure
corresponds to a particular graph. We added to this figure a regression line
whose slope is close to 1. The dependence is presented in logarithmic scale.
We do not expect a more precise relationship between $|V|+|E|$ and the running
time since the implementation is far from being optimized and since it may
depend also on structural factors such as the degree distribution of the graph
and its diameter.

Figure 7: ms-GMLogA: running time vs graph size.

###  4.2 Spectral approach

Arrangement of the graph vertices according to the eigenvector (called the
Fiedler vector) corresponding to the second smallest eigenvalue of the graph
Laplacian is a well-known and successful heuristic used for many ordering
problems such as the minimum $p$-sum and the minimum envelope reduction (see
[4, 15]). Usually, for large graphs the computation of the Fiedler vector is
too expensive, and some approximation is employed [21]. Most of these
approximations can be viewed as a global averaging process that works until a
particular convergence. As a result, during this process the vertex tends to
be located at the weighted average of its neighbors, which makes it suitable
for solving the quadratic functionals such as the minimum $2$-sum problem [15]
but creates a poor solution for the sum of logarithms. Several examples (from
[17]) of a comparison between spectral approach and ms-GMLogA are presented in
Table 1.

Table 1: Comparison of spectral method and ms-GMLogA Network | Spectral | ms-GMLogA
---|---|---
as-caida20071105 | 12.4471 | 6.32
email-Enron | 10.3117 | 7.06
oregon1_010407 | 11.0662 | 6.18
p2p-Gnutella06 | 10.6513 | 7.13
wiki-Vote | 9.2152 | 7.41

###  4.3 Compressing the ordered networks

In [11], Chierichetti et al. conjectured that minimizing the unweighted
version of MLogA will automaticaly improve the graph compression. The
conjecture has been supported by experimental observations on few web graphs
and social networks. However, for general networks, a better ordering does not
always imply a better compression ratio produced by the kind of compressions
described in [5]. To illustrate that we present in Figure 8 a dependence of
order improvement on the compression produced in [5].

Given an initial ordering (natural or Gray), we apply our algorithm to obtain
a new ordering with smaller $\beta$. The x-axis corresponds to the ratio
between $\beta$ values of initial ordering and ms-GMLogA. The y-axis
corresponds to the ratio obtained by compressing these orders. On most of the
graphs, the conjecture is confirmed. However, on some graphs we can see the
degradation of the compression results. It is observable, in particular, in
Figure 8(b), where the initial arrangement was obtained by Gray ordering that
is significantly better than the natural ordering.

a) initial ordering is natural

b) initial ordering is Gray

Figure 8: Correlation between ordering and compression improvements

##  5 Conclusions and future work

We have proposed a fast linear algorithm (ms-GMLogA), for compression-friendly
graph reordering. The algorithm belongs to the family of multiscale methods
[9]. It uses a novel AMG-based coarsening scheme that is reinforced by a
modification of algebraic distance couplings [20]. The empirical results on a
large benchmark demonstrate its quality and scalability. The model takes into
account the edge and node weights that can express different properties of a
network model such as known link importance or its access frequency. Overall,
we recommend this multiscale framework for practical, compression-friendly
network orderings.

As future work we identify the following research directions that have the
potential to improve ms-GMLogA:

  * •

development of a learning algorithm to improve the bandwidth $h$ in (8);

  * •

development of a more sophisticated collective refinement of nodes;

  * •

development of a multiscale method for MLinGapA (defined below); and

  * •

more sophisticated multiscale organization of directed graphs.

A more sophisticated collective refinement of nodes can allow interruption of
the coarsening earlier. Thus, if some more qualitative approximation can be
achieved at the very coarse levels instead of solving exactly the coarsest and
applying fast simple refinement at the next few finer level, it can
potentially lead to better solutions at the fine levels.

The next natural step is to design first a refinement and then a full
multiscale approach for the minimum logarithmic gap arrangement problem
(MLinGapA) defined in [11] as an alternative model for compression. In this
problem the goal is to find a permutation of nodes that minimizes

$\sum_{i\in V}f_{\pi}(i,out(i)),$ (10)

where $f_{\pi}(i,out(i))$ is a sum of logarithms of gaps between consecutive
neighbors of $i$ ordered by $\pi$. This problem is also NP-hard [11]; however,
we believe that having a suitable refinement and relaxation algorithms would
make it possible to adopt the multiscale framework for this problem as well.

The last research direction is related to finding a better multiscale scheme
for directed graphs. The presented model constructs an undirected edge $ij$
from a given set of directed edges between $i$ and $j$ by reflecting both
directions as an edge weight. In fact, this way of representation is not far
from the attempts to solve the problems on nonsymmetric matrix $A$ by applying
the known techniques on $A+A^{\prime}$ or $A\cdot A^{\prime}$. It is known
that these methods suffer from several drawbacks, however; thus, we identify a
finding of a more advanced multiscale representation for directed graphs as
one of the major future research directions.

## Acknowledgments

This work was partially funded by the CSCAPES institute, a DOE project, and in
part by DOE Contract DE-AC02-06CH11357. We express our gratitude to Marco Rosa
for helping to design experiments with Webgraph package [1].

## Appendix A: Parameters

$\Theta_{1,2}$. The threshold parameters $\Theta_{1}$ and $\Theta_{2}$ from
(4) are responsible for a control of complexity of the coarse level problem
and the quality of coarsening. It is important to check their robustness when
designing a multiscale framework. The presented computational experiments have
been executed with $\Theta_{1}=\Theta_{2}=1/2$, however, other values
$.2\leq\Theta_{1},\Theta_{2}\leq.8$ exhibit a good behavior of numerical
results as well. Of course, the larger values increase the running time as the
coarse graphs become bigger.

$\omega$. Detailed discussion about choosing parameter $\omega$ for Jacobi
over-relaxation on general graphs is discussed in [10]. In all our
computational experiments this parameter was 0.5.

## References

  * [1] Laboratory for Web Algorithmics.  http://law.dsi.unimi.it/.
  * [2] Micah Adler and Michael Mitzenmacher.  Towards compressing web graphs.  In DCC ’01: Proceedings of the Data Compression Conference, page 203, Washington, DC, USA, 2001. IEEE Computer Society.
  * [3] Alberto Apostolico and Guido Drovandi.  Graph compression by BFS.  Algorithms, 2(3):1031–1044, 2009.
  * [4] S.T. Barnard, A. Pothen, and H. Simon.  A spectral algorithm for envelope reduction of sparse matrices.  Numerical Linear Algebra with Applications, 2(4):317–334, 1995\.
  * [5] P. Boldi and S. Vigna.  The webgraph framework I: compression techniques.  In WWW ’04: Proceedings of the 13th International Conference on World Wide Web, pages 595–602, New York, NY, USA, 2004. ACM.
  * [6] Paolo Boldi, Massimo Santini, and Sebastiano Vigna.  Permuting web graphs.  In Konstantin Avrachenkov, Debora Donato, and Nelly Litvak, editors, WAW, volume 5427 of Lecture Notes in Computer Science, pages 116–126. Springer, 2009.
  * [7] Erik G Boman, Umit V Catalyurek, CÃ©edric Chevalier, Karen D Devine, Ilya Safro, and Michael M Wolf.  Advances in parallel partitioning, load balancing and matrix ordering for scientific computing.  Journal of Physics: Conference Series, 180(1), 2009.
  * [8] A. Brandt, S. McCormick, and J. Ruge.  Algebraic multigrid (AMG) for automatic multigrid solution with application to geodetic computations.  Technical report, Institute for Computational Studies, Fort Collins, CO, POB 1852, 1982.
  * [9] A. Brandt and D. Ron.  Chapter 1 : Multigrid solvers and multilevel optimization strategies.  In J. Cong and J. R. Shinnerl, editors, Multilevel Optimization and VLSICAD. Kluwer, 2003.
  * [10] Jie Chen and Ilya Safro.  Algebraic distance on graphs.  Technical Report ANL/MCS-P1696-1009 (Preprint), Argonne National Laboratory.
  * [11] Flavio Chierichetti, Ravi Kumar, Silvio Lattanzi, Michael Mitzenmacher, Alessandro Panconesi, and Prabhakar Raghavan.  On compressing social networks.  In KDD ’09: Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 219–228, New York, NY, USA, 2009. ACM.
  * [12] Yongwook Choi and Wojciech Szpankowski.  Compression of graphical structures.  In ISIT’09: Proceedings of the 2009 IEEE international conference on Symposium on Information Theory, pages 364–368, Piscataway, NJ, USA, 2009. IEEE Press.
  * [13] T. Davis.  University of Florida Sparse Matrix Collection.  NA Digest, 97(23), 1997.
  * [14] J. Díaz, J. Petit, and M. Serna.  A survey of graph layout problems.  ACM Comput. Surv., 34(3):313–356, 2002.
  * [15] Martin Juvan and Bojan Mohar.  Optimal linear labelings and eigenvalues of graphs.  Discrete Appl. Math., 36(2):153–168, 1992.
  * [16] Y. Lai and K. Williams.  A survey of solved problems and applications on bandwidth, edgesum, and profile of graphs.  J. Graph Theory, 31:75–94, 1999.
  * [17] J. Leskovec.  Stanford Network Analysis Package (SNAP).  http://snap.stanford.edu/index.html.
  * [18] Jure Leskovec, Lars Backstrom, Ravi Kumar, and Andrew Tomkins.  Microscopic evolution of social networks.  In KDD ’08: Proceeding of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 462–470, New York, NY, USA, 2008. ACM.
  * [19] Renfrey B. Potts.  Some generalized order-disorder transformations.  Proceedings of the Cambridge Philosophical Society, 48:106–109, 1952.
  * [20] Dorit Ron, Ilya Safro, and Achi Brandt.  Relaxation-based coarsening and multiscale graph organization.  Technical Report ANL/MCS-P1741-0410 (Preprint), Argonne National Laboratory, 2010.
  * [21] Y. Saad.  Numerical Methods for Large Eigenvalue Problems.  Halsted Press, 1992.
  * [22] I. Safro, D. Ron, and A. Brandt.  Graph minimum linear arrangement by multilevel weighted edge contractions.  Journal of Algorithms, 60(1):24–41, 2006.
  * [23] I. Safro, D. Ron, and A. Brandt.  Multilevel algorithms for linear ordering problems.  Journal of Experimental Algorithmics, 13:1.4–1.20, 2008.
  * [24] I. Safro and B. Temkin.  Benchmark for the Minimum Logarithmic Arrangement Problem.  http://www.mcs.anl.gov/ safro/mloga.html.
  * [25] B. W. Silverman.  Density Estimation for Statistics and Data Analysis.  Chapman & Hall, London, 1986.
  * [26] Jie Sun, Erik M Bollt, and Daniel Ben-Avraham.  Graph compression-save information by exploiting redundancy.  Journal of Statistical Mechanics: Theory and Experiment, 2008(06):P06001, 2008.
  * [27] Chris Walshaw.  Multilevel Refinement for Combinatorial Optimisation: Boosting Metaheuristic Performance, pages 4261–289.  Springer, 2008.
  * [28] H. Zhang, B. Qiu, K. Ivanova, C. L. Giles, H. C. Foley, and J. Yen.  Locality and attachedness-based temporal social network growth dynamics analysis.  Journal of American Society of Information Science and Technology, 61:964–977, 2009.

The submitted manuscript has been created in part by UChicago Argonne, LLC,
Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S.
Department of Energy Office of Science laboratory, is operated under Contract
No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others
acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in
said article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government.