# Mode-Assisted Joint Training of Deep Boltzmann Machines

Haik Manukian  email<EMAIL_ADDRESS>Department of Physics, University of
California, San Diego, La Jolla, CA 92093  Massimiliano Di Ventra  email:
<EMAIL_ADDRESS>Department of Physics, University of California, San
Diego, La Jolla, CA 92093

###### Abstract

The deep extension of the restricted Boltzmann machine (RBM), known as the
deep Boltzmann machine (DBM), is an expressive family of machine learning
models which can serve as compact representations of complex probability
distributions. However, jointly training DBMs in the unsupervised setting has
proven to be a formidable task. A recent technique we have proposed, called
mode-assisted training, has shown great success in improving the unsupervised
training of RBMs. Here, we show that the performance gains of the mode-
assisted training are even more dramatic for DBMs. In fact, DBMs jointly
trained with the mode-assisted algorithm can represent the same data set with
orders of magnitude lower number of total parameters compared to state-of-the-
art training procedures and even with respect to RBMs, provided a fan-in
network topology is also introduced. This substantial saving in number of
parameters makes this training method very appealing also for hardware
implementations.

One of the most highly influential models in Artificial Intelligence (AI) and
deep learning in particular is the Boltzmann machine (BM). They were
constructedackley1985learning ; smolensky1986information  as a powerful
stochastic generalization of Hopfield networks hopfield1982neural  that
possess a simple expression for their log-likelihood gradient. However, they
remained impractical to train due to their reliance on high-dimensional
sampling to calculate that gradient. A relatively efficient learning
algorithm, called contrastive divergence (CD), was discovered for BMs with a
simplified topology called restricted Boltzmann machines
(RBMs)hinton2002training , which have since gone on to see success in various
domains goodfellow2016deep . However, the extension of RBMs to deep Boltzmann
machines (DBMs) has been difficult goodfellow2012scaling , and as such, DBMs
are now mostly overshadowed by their deep feedforward cousins
goodfellow2014generative  in generative applications.

This is not due to DBM’s lack of ability, but rather the absence of effective
means to train these models. There remain quite a few reasons to search for
better learning algorithms for DBMs, as they are a versatile computational
medium. A principle use is as compact generative models for complex
probability distributions in unsupervised settings, considered a critical
component of the forthcoming “third-wave” of AI launchbury2017darpa . Trained
DBMs can also serve as an informed prior for feedforward networks, leading to
better generalization in supervised tasks erhan2010does . In the physical
sciences, DBMs serve as powerful variational representations of many body
wavefunctions, more efficiently gao2017efficient  than RBMs carleo2017solving
; melko2019restricted , and have potential applications in condensed matter
physics and quantum computing carleo2019machine .

Figure 1: Schematic of a deep Boltzmann Machine with a visible layer, $v$, and
hidden layers, $h^{(j)}$, with $j=1,2,\dots$. Connections between nodes are
symmetric and undirected, in contrast to typical directed feedforward
networks.

Various attempts have been made to climb the summit of DBM training
salakhutdinov2009deep ; salakhutdinov2010efficient ; hinton2012better ;
goodfellow2013joint ; melchior2016center . Most approaches rely on pre-
training by breaking up layers into RBMs and training them sequentially, after
which the DBM is fine-tuned jointly. Correlations between layers are ignored
during pre-training, minimizing the potential advantages of the deep
architecture and can be disruptive to joint training melchior2016center .
Recently, the authors introduced mode-assisted training AcceleratingDL ;
manukian2020mode , which combines CD with samples of the model distribution
mode. This stabilizes training, allows the learning of very accurate
densities, and strikes a better tradeoff between accuracy and computational
cost compared to CD. As the method is agnostic to the connectivity of the
network, we can apply mode-assistance to DBMs with the hope of capturing the
model capacity missed by other approaches.

In this work, we find that the benefits of mode-assisted training are even
more dramatic in the case of DBMs. In fact, it produces more accurate models
without requiring pre-training while also utilizing orders of magnitude less
parameters, compared to pre-trained salakhutdinov2009deep  or centered DBMs
melchior2016center . The role of the network topology is also discussed, where
we discover that DBMs are easier to train if the size of the layers decreases
with depth. We evaluate the density modeling performance of mode-assisted DBMs
by computing exact log-likelihoods achieved on small data sets and
approximating the likelihood on the MNIST data set. The approach we propose
can be extended to other types of neural networks, and is relevant also for
the hardware implementation of these models, where a much smaller number of
parameters directly translates into components and energy savings.

To see how mode-assisted training can be extended to deep architectures, we
give a quick overview of DBMs and the basic approach to their training. DBMs
are undirected weighted graphs that differentiate between $n_{v}$ visible
nodes, and $\ell$ layers of $n_{\ell}$ latent, or ‘hidden’, nodes, not
directly constrained by the data salakhutdinov2009deep . We assume that
$\ell>1$ as $\ell=1$ recovers the RBM, and like an RBM, there are no
connections within a layer. Each state of the machine corresponds to an energy
of the form

$\displaystyle E({\bf v},{\bf h}^{(1)},$ $\displaystyle\cdots,{\bf
h}^{(\ell)})=-{\bf a}^{T}{\bf v}-{\bf v}^{T}{\bf W}^{(1)}{\bf h}^{(1)}$ (1)
$\displaystyle-\sum_{i=1}^{\ell}{\bf b}^{(i)T}{\bf
h}^{(i)}-\sum_{i=2}^{\ell}{\bf h}^{(i-1)T}{\bf W}^{(i-1)}{\bf h}^{(i)},$

where the biases ${\bf a}\in\mathbb{R}^{n_{0}}$, ${\bf
b}^{\ell}\in\mathbb{R}^{n_{\ell}}$, and weights ${\bf
W}^{\ell}\in\mathbb{R}^{n_{\ell-1}\times n_{\ell}}$ are the learnable
parameters. The energy function in Eq. (1) induces a Boltzmann-Gibbs
distribution over states,

$p({\bf x})=\frac{e^{-E({\bf x})}}{\mathcal{Z}},$ (2)

where ${\bf x}=({\bf v},{\bf h}^{(1)},\dots,{\bf h}^{(\ell)})$. The partition
function, $\mathcal{Z}=\sum_{\\{{\bf x}\\}}e^{-E({\bf x})}$, involves the sum
of an exponentially growing number of states, making the exact computation of
its value infeasible for most data sets.

During learning, a DBM is tasked to match its marginal distribution over the
visible layer, $p({\bf v})=\sum_{\\{{\bf h}\\}}p({\bf v},{\bf h})$, to an
unknown data distribution, $q({\bf v})$, represented by a data set,
$\mathcal{D}$. Training a DBM amounts to a search for the appropriate weights
and biases that will minimize the quantity known as the Kullback-Leibler (KL)
divergence between the two distributions,

$\text{KL}(q||p)=\sum_{\\{{\bf v}\\}}q({\bf v})\log\frac{q({\bf v})}{p({\bf
v})},$ (3)

or, equivalently, maximizing the log-likelihood of the dataset,
$\textrm{LL}(p)=\sum_{{\bf v}\in\mathcal{D}}\log p({\bf v})$. The optimization
of the non-linear, and typically high dimensional Eq. (3) (or log-likelihood),
is often done via stochastic gradient descent with respect to the DBM
parameters, which leads to weight updates of the form fischer2012introduction
,

$\Delta w_{ij}\propto\langle x_{i}x_{j}\rangle_{D}-\langle
x_{i}x_{j}\rangle_{M}.$ (4)

For every gradient update in Eq. (4), nodal statistics must be computed under
two different distributions. The first one on the RHS of Eq. (4) is called the
“data term”, and is an expectation over the data induced distribution, $q({\bf
v})p({\bf h}|{\bf v})$, with the network’s visible layer fixed to the data.
The second term on the RHS of Eq. (4) is called the “model term” which is an
expectation over the entire model distribution in Eq. (2). In the case of
RBMs, the data term can be sampled from exactly, but the model term must be
approximated. With DBMs, the data term must also be approximated, most
popularly with an iterative mean field procedure (see Methods).

In both cases, model statistics are collected via a Markov Chain Monte Carlo
(MCMC) procedure dubbed ‘contrastive divergence’ (CD) hinton2002training .
CD-$k$ is a form of Gibbs sampling that initializes chains of length $k$ from
elements of the dataset. Trouble arises when the model distribution contains
’spurious’ modes where the data distribution has negligible probability. In
these cases, ergodicity breaks down, and mixing times become prohibitively
long, frequently resulting in CD becoming biased enough to cause training to
diverge fischer2010empirical . Training a DBM jointly with CD has proven to be
a formidable task. Even a two-layer DBM on MNIST has not seen success without
some kind of modification. salakhutdinov2009deep ; melchior2016center ;
goodfellow2012scaling ; desjardins2012training

Here, instead, we use mode-assisted training in the joint and unsupervised
learning of DBMs. The essence of mode-assisted training is the replacement of
some gradient updates in Eq. 4 with ones of the form,

$\Delta w_{ij}\propto[x_{i}x_{j}]_{D}-[x_{i}x_{j}]_{M}$ (5)

Where the notation $[f({\bf x})]_{q}$ represents $f({\bf x}_{\text{mode}})$,
evaluated at ${\bf x}_{\text{mode}}$, the mode of some distribution, $q({\bf
x})$. One may employ any optimization solver to sample the mode of the above
distributions. Due to its proven efficiency, here we employ a memcomputing one
as reported in our previous work manukian2020mode .

The mode-assisted update can be thought of as a saddle-point approximation of
an expectation bender2013advanced . This technique, also known as Laplace’s
method, is commonly employed to approximate integrals (expectations) of the
form,

$\langle f({\bf x})\rangle=\frac{\int e^{-E({\bf x})}f({\bf x})d{\bf x}}{\int
e^{-E({\bf x})}d{\bf x}}\approx f({\bf x}_{\text{mode}}).$ (6)

Figure 2: Average converged log-likelihood performance (lower bound) between
RBMs trained with CD-1, mode-assisted DBMs (MA), and unassisted DBMs with
CD-1. The DBMs have two hidden layers. The networks were trained on the
shifting bar data set with $n_{v}=24$ (left plot) and $n_{v}=12$ (right plot)
for 200,000 and 100,000 gradient updates respectively, following a linearly
decaying learning rate schedule from $\epsilon=1\to 0.001$. For the DBMs the
hidden layer ratio was fixed at $\alpha=n_{h^{(2)}}/n_{h^{(1)}}=0.2$.
Performance is shown as a function of total number of hidden nodes,
$n_{h}=n_{h^{(1)}}+n_{h^{(2)}}$. The solid lines are the median obtained
across an ensemble of 50 networks, and the shaded regions enclose the 95th and
5th percentiles.

The weight updates driven by the mode are incorporated in a probabilistic way,
with the probability of a mode driven update following a sigmoid, starting low
in the initial phases of the training and reaching a maximal value at the end
of training:

$P_{\text{mode}}(n)=P_{\text{max}}\sigma(\alpha n+\beta).$ (7)

Here, $n$ is the current epoch, and $\alpha$, $\beta$, $P_{\text{max}}$
control the shape of the sigmoid. Throughout the work we set, $\alpha=20/N$
($N$ is the total number of epochs), $\beta=-6$ and $P_{\text{max}}=0.1$. This
schedule was introduced in Ref. manukian2020mode, , based on the empirical
observation that the mode samples work best when the support has been
‘discovered’ by CD.

The key insight of mode assisted training is not to approximate the
expectation, but rather directly prevent spurious modes from appearing in the
model distribution. This allows the inexact but efficient approximation of
CD-$k$ to minimize the KL divergence (or maximize the log-likelihood) without
diverging. Searching for the mode with a highly specialized optimizer results
in a better trade-off between log-likelihood performance and computation time.

## Results

To illuminate the effectiveness of mode assistance on DBMs, we first evaluate
performance on two differently sized synthetic shifting bar data sets. In Fig.
2, the mode-assisted algorithm on a two-layer DBM is compared to two
baselines, CD on the same DBM and to CD on an RBM with the same number of
hidden nodes. The converged log-likelihoods are reported, and results are
shown as a function of the total number of hidden nodes. Overall, mode-
assisted training almost always converges to more accurate densities than CD
alone on a DBM, and in most cases also does better than the corresponding RBM
with CD. Mode assistance is also seen to prevent the divergence sometimes seen
with Gibbs sampling, as well as reducing the variance of the converged models,
resulting in smaller error bars. Past a certain point, all methods expectedly
incur a loss in performance as the number of hidden nodes increases for a
fixed number of training iterations. Mode-assisted training suffers the least
in this regard.

Although the DBMs in Fig. 2 possess the same total number of hidden nodes as
the RBMs, they contain fewer total parameters. This means that a better
trained DBM takes advantage of abstract features composition afforded by the
depth of the network, mirroring similar gains found in deep feedforward neural
networks compared to single layer perceptrons goodfellow2016deep .

Note, however, that this parameter efficiency in DBMs is not present when
training jointly with CD. In fact, they perform systematically worse than an
equivalently sized RBM in terms of log-likelihood. Mode assisted DBMs on the
other hand perform as well as RBMs or better, all the while maintaining
parameter efficiency.

To quantify this efficiency gain, let us consider the case of a DBM with two
hidden layers. If the total number of hidden nodes is fixed,
$n_{h}=n_{h^{(1)}}+n_{h^{(2)}}$, then a measure of parameter efficiency
compared to an RBM with the same number of nodes is captured by the parameter
$\alpha=n_{h^{(2)}}/n_{h^{(1)}}$. The total number of parameters (weights and
biases) in an RBM with $n_{v}$ visible nodes and $n_{h}$ hidden nodes is,

$n_{\text{RBM}}=n_{v}n_{h}+n_{v}+n_{h}.$

A DBM with the same total number of hidden nodes would have the following
number of parameters

$n_{\text{DBM}}=n_{v}n_{h^{(1)}}+n_{h^{(1)}}n_{h^{(2)}}+n_{v}+n_{h^{(1)}}+n_{h^{(2)}}.$

Considering typical training scenarios of a large visible layer, $n_{v}\gg
n_{h}\gg 1$, the fraction of total parameters used in the DBM compared to the
RBM (the parameter efficiency) can be simplified as,

$e=\frac{n_{\text{DBM}}}{n_{\text{RBM}}}\approx\frac{1}{1+\alpha}.$ (8)

Figure 3: Average Log-Likelihood achieved on an $n_{v}=24$ dimensional
shifting bar data set as a function of DBM shape, $\alpha$. The total number
of hidden nodes are kept fixed at $n_{h^{(1)}}+n_{h^{(2)}}=22$. An ensemble of
50 DBMs were trained with CD-1 and MA, where the solid line shows the median
average log-likelihood achieved after $10^{5}$ gradient iterations with a
linearly decaying learning rate $\epsilon=1\to 0.01$, and shaded regions
delineate the 95th and 5th percentiles.

The smaller $e$, the fewer parameters a DBM possesses with respect to an RBM
with the same number of nodes. We then see that a naïve interpretation of Eq.
(8) would suggest that the DBM parameter efficiency becomes more significant
with increasing $\alpha$, meaning the resulting DBM has more neurons in the
deeper layers (fans out). However, this assumes equal log-likelihood
(performance) of the two models, which is not obvious. In fact, for DBMs we
expect that a redistribution of the relative number of nodes in the different
hidden layers plays a significant role in their performance. This tradeoff
between parameter efficiency and performance merits further investigation.

This is shown in Fig. 3, where DBMs trained with and without mode assistance
are compared as a function of $\alpha$. For all cases, the mode-assisted
algorithm performs better than its unassisted counterpart, which is nothing
other than the confirmation of the results of Fig. 2. Importantly, however, we
also observe that the log-likelihood performance increases as the network
becomes less parameter efficient compared to an RBM $(\alpha\to 0)$, namely
when it has a fan-in topology. We find that a balance is reached around
$\alpha\sim 0.15$, namely the second hidden layer has only about 15% of nodes
than the first hidden layer. The message here is that for a fixed number of
hidden nodes, it is best to organize the network as a fan-in DBM rather than
an RBM, if one has a method to train the DBM close to capacity in the first
place.

Figure 4: Average log-likelihood on the MNIST data set achieved after 100
epochs with mode-assisted training (MA) compared to the recent best achieved
results on DBMs using CD as well as pre-training. DBMs trained with MA were
all with fixed $\alpha=0.15$ and an increasing number of hidden nodes,
$n_{h}$. The learning rate was chosen to follow a linear decay,
$\epsilon=0.05\to 0.0005$ and no pre-training was employed. The best log-
likelihood was reported out of 10 randomly initialized runs.

Finally, to show the substantial improvements provided by the mode-assisted
training of DBMs, we compare it to the training of centered
DBMsmelchior2016center  and pre-trained DBMssalakhutdinov2009deep  using CD on
the MNIST data set. In Fig. 4, we report similar trends we observed in the
smaller data sets, but scaled to more dramatic results. As the number of
hidden nodes increases, small mode-assisted DBMs surpass the performance of
significantly larger standard DBMs, centered DBMs and eclipse the performance
of much larger pre-trained DBMs. In fact, with a DBM of only $120\times 18$
hidden nodes (and no pre-training), trained with our mode-assisted approach,
we reach the same log-likelihood of a DBM with $500\times 1000$ hidden nodes
and pre-training – a parameter savings of two orders of magnitude.

For a fair comparison, the total training iterations were set to 100 epochs in
each case, and the log-likelihood results with our approach are shown as a
function of total number of hidden nodes at a fixed ratio, $\alpha=0.15$. The
partition function was approximated using Annealed Importance Sampling (AIS)
(see Methods), identical to the other referenced results melchior2016center .
It’s worth noting that AIS tends to under-approximate the partition function,
leading to an over-estimate of the log-likelihood. This effect is exaggerated
in larger models. Even with this negative impact, a mode-assisted DBM with a
total of 138 nodes without pre-training achieves about the same performance as
one with 1,500 that was pre-trained.

## Discussion

We conclude by providing an intuitive understanding of the dramatic
performance improvement of the mode-assisted training over Gibbs sampling
(CD). The standard method of pretraining DBMs initializes weights at a greedy
starting point but critically ignores correlations within the system. The end
result is a dramatic loss of parameter efficiency compared to mode-assisted
joint training of DBMs. Even training methods like centering, with more
advanced sampling techniques like parallel tempering swendsen1986replica ,
fail to train the DBM near its capacity. The issue is the reliance on Gibbs-
like sampling to explore the states of the DBM, and using only these
statistics to compute the likelihood gradient. The breakdown occurs when Gibbs
chains fail to explore adequately the states of the system.

Because of its random-walk nature, the local dynamics of Gibbs sampling
suffers long auto-correlation times when equally likely states are separated
by an extensive number of variables flips. As a consequence of long mixing
times between such states, statistics become heavily biased. In the initial
phase of DBM training, when weights are small and randomly distributed, these
issues are minimized. Interactions between nodes in this ‘high temperature’
regime are weak, and are well approximated by a mean-field theory
sherrington1975solvable . In this phase, CD works well without much bias, and
is the reason mode-assistance is really not necessary early on in the
training.

However, as training continues, the DBM is effectively ‘cooled’ as weights
learn from the data and grow larger in magnitude. Correlations between nodes
become significant, and can no longer be ignored or averaged over: mean-field
theory is inadequate to describe this phase. Gibbs sampling can then
dramatically fail in these conditions, that is why mode-assistance is
primarily applied in this phase. During mode updates, information is
propagated from the ground state, the most ‘collective’ or ‘coordinated’ state
of the DBM, to all the weights. This prevents probability mass density from
accumulating far away from the current state of the chain, keeping Gibbs
sampling in its effective regime.

Aside from mode-assistance, we have discovered that the topology of a DBM
plays a major role in its performance as well as its parameter efficiency. Too
many hidden nodes act as a noise source during training, slowing down
learning. Too few hidden nodes reduce the capacity of the network. For the
network sizes considered, a ‘sweet spot’ was found where these two effects are
in balance. At this stage, it is not clear if this is particular to DBMs, or
there are deeper reasons why this is the case, and further work in this
direction would be interesting.

Finally, mode training relies on an efficient search for the mode of a DBM,
which is an NP-hard computational problem on its own. On this front, we employ
a novel dynamical systems based optimization technique called memcomputing,
which has shown great promise in efficient optimization of non-convex energy
landscapes like the DBM sheldon2019taming ; manukian2020mode . Since the DBM
energy landscape can be represented within the optimizing dynamics, a path
exists to extend these systems to (learning) samplers. In doing so, the need
for CD would be entirely eliminated. We leave these avenues of research to be
explored in future work.

## Methods

DBMs are powerful models in machine learning, but bring with them considerable
computational burdens during evaluation. To deal with them, we follow the
procedure outlined in Ref. salakhutdinov2009deep, . First, due to additional
complexity, the log-likelihood is not directly maximized as in the case of
RBMs. Instead, a variational lower bound to the log-likelihood is optimized.
Second, the most common variational form chosen leads to mean-field updates
for the data term. Finally, for evaluation of log-likelihood lower bounds, an
approximation of the partition function is necessary. We provide a short
description of all these approximations which are frequently encountered in
DBM training scenarios.

Likelihood Lower Bound \- Since the data expectation in Eq. (4) is no longer
in closed form for the DBM, data-dependent statistics must be approximated
with a sampling technique over the conditional distribution, $p({\bf h}|{\bf
v})$, where ${\bf h}=\\{{\bf h}^{(1)},\cdots,{\bf h}^{(\ell)}\\}$. In
practice, a variational lower bound to the log-likelihood is maximized
instead, which is tractable and is found to work well (as in the model term)
neal1998view ; salakhutdinov2009deep ; salakhutdinov2012efficient .

The variational approximation replaces the original posterior distribution,
$p({\bf h}|{\bf v})$ with an approximate distribution, $r({\bf h}|{\bf v})$,
where the parameters of $r$ follow the gradient of the resulting lower bound
on the original log-likelihood neal1998view ,

$\displaystyle\log p({\bf v})$ $\displaystyle\geq\sum_{{\bf h}}r({\bf h}|{\bf
v})\log p({\bf v},{\bf h})+H_{r}({\bf v})$ (9) $\displaystyle=\log p({\bf
v})-\textrm{KL}(r({\bf h}|{\bf v})||p({\bf h}|{\bf v})),$

where $H_{r}({\bf v})=-\sum_{\bf h}r({\bf h}|{\bf v})\log r({\bf h}|{\bf v})$
is the Shannon entropy of $r$. This variational loss simultaneously attempts
to maximize the log-likelihood of the data set and minimize the KL divergence
between the true conditional distribution, $p$ and its approximation, $r$.

Data Term \- A fully factorial mean field ansatz is often used in the
variational approach, $r({\bf h}|{\bf v})=\prod_{i}r(h_{i}|{\bf v})$ with
$r(h_{i}=1|{\bf v})=\mu_{i}$, with $\mu_{i}\in[0,1]$ randomly initialized from
a uniform distribution. Maximizing the lower bound in Eq. (9) results in
updates of the form,

$\mu_{i}\leftarrow\sigma\left(\sum_{j}W_{ij}^{(1)}v_{j}+\sum_{i\neq
j}J_{ij}\mu_{i}+b_{i}\right).$ (10)

Here $J$ is a block matrix containing the weights of the hidden nodes and
$b_{i}$ is the bias of the $i$-th hidden node. Convergence is typically fast,
(in our experiments less than 30 iterations are enough) and these states are
then used to calculate the data expectation in Eq. (4) during the Gibbs phase
of the training.

Partition Function \- Computing the partition function in Eq. (2) exactly is
infeasible for large DBMs. Its value is only required to evaluate the
performance of the networks and appears in the log-likelihood as a normalizing
constant. Annealed Importance Sampling (AIS) salakhutdinov2008learning  is the
procedure often used to approximate the partition function of large RBMs and
DBMs. AIS estimates the ratio of partition functions, $Z_{N}/Z_{0}$, using a
sequence of probability distributions between a chosen initial distribution
and the desired one. The initial distribution, $p_{0}$, is chosen to have an
exactly known partition function and to be simple to sample from (e.g.
uniform) and $p_{N}$ is the desired distribution whose partition function one
wants to compute.

The sequence in the case of DBMs is parameterized by $\beta_{k}$ (inverse
temperatures), giving $p_{k}({\bf x})=e^{-\beta_{k}E({\bf x})}/Z_{k}$. Markov
chains are initialized uniformly according to $p_{0}$ and $0\leq\beta_{k}\leq
1$ is slowly annealed to unity according to a desired schedule, all the while
allowing the chains to run. The ratio is then approximated by the product of
ratios of the unnormalized intermediate distributions,

$\frac{Z_{N}}{Z_{0}}=\frac{p_{1}^{*}({\bf x})}{p_{0}^{*}({\bf
x})}\cdots\frac{p_{k}^{*}({\bf x})}{p_{k-1}^{*}({\bf x})}.$ (11)

When dealing with bipartite graphs like DBMs, either all the even or all the
odd layers can be analytically traced out, resulting in a tighter
approximation. For a direct comparison with previous work
salakhutdinov2009deep ; melchior2016center , we average over an identical
number of 100 AIS runs with 29,000 linearly spaced intermediate distributions.

Acknowledgments
Work supported by DARPA under grant No. HR00111990069. H.M. acknowledges
support from a DoD-SMART fellowship. All memcomputing simulations reported in
this paper have been done on a single core of an AMD EPYC server.

## References

  * (1) Ackley, D. H., Hinton, G. E. & Sejnowski, T. J.  A learning algorithm for boltzmann machines.  _Cognitive Science_ 9, 147–169 (1985).
  * (2) Smolensky, P.  Information processing in dynamical systems: Foundations of harmony theory (1986).
  * (3) Hopfield, J. J.  Neural networks and physical systems with emergent collective computational abilities.  _Proceedings of the national academy of sciences_ 79, 2554–2558 (1982).
  * (4) Hinton, G. E.  Training products of experts by minimizing contrastive divergence.  _Neural Computation_ 14, 1771–1800 (2002).
  * (5) Goodfellow, I., Bengio, Y., Courville, A. & Bengio, Y.  _Deep learning_ , vol. 1 (MIT Press, Cambridge, 2016).
  * (6) Goodfellow, I. J., Courville, A. & Bengio, Y.  Scaling up spike-and-slab models for unsupervised feature learning.  _IEEE transactions on pattern analysis and machine intelligence_ 35, 1902–1914 (2012).
  * (7) Goodfellow, I. J. _et al._ Generative adversarial networks.  _arXiv preprint arXiv:1406.2661_ (2014).
  * (8) Launchbury, J.  A darpa perspective on artificial intelligence.  _Retrieved November_ 11, 2019 (2017).
  * (9) Erhan, D. _et al._ Why does unsupervised pre-training help deep learning?  _Journal of Machine Learning Research_ 11, 625–660 (2010).
  * (10) Gao, X. & Duan, L.-M.  Efficient representation of quantum many-body states with deep neural networks.  _Nature communications_ 8, 1–6 (2017).
  * (11) Carleo, G. & Troyer, M.  Solving the quantum many-body problem with artificial neural networks.  _Science_ 355, 602–606 (2017).
  * (12) Melko, R. G., Carleo, G., Carrasquilla, J. & Cirac, J. I.  Restricted boltzmann machines in quantum physics.  _Nature Physics_ 15, 887–892 (2019).
  * (13) Carleo, G. _et al._ Machine learning and the physical sciences.  _Reviews of Modern Physics_ 91, 045002 (2019).
  * (14) Salakhutdinov, R. & Hinton, G.  Deep boltzmann machines.  In _Artificial intelligence and statistics_ , 448–455 (2009).
  * (15) Salakhutdinov, R. & Larochelle, H.  Efficient learning of deep boltzmann machines.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , 693–700 (2010).
  * (16) Hinton, G. E. & Salakhutdinov, R. R.  A better way to pretrain deep boltzmann machines.  In _Advances in Neural Information Processing Systems_ , 2447–2455 (2012).
  * (17) Goodfellow, I. J., Courville, A. & Bengio, Y.  Joint training deep boltzmann machines for classification.  _arXiv preprint arXiv:1301.3568_ (2013).
  * (18) Melchior, J., Fischer, A. & Wiskott, L.  How to center deep boltzmann machines.  _The Journal of Machine Learning Research_ 17, 3387–3447 (2016).
  * (19) Manukian, H., Traversa, F. L. & Di Ventra, M.  Accelerating deep learning with memcomputing.  _Neural Networks_ 110, 1–7 (2019).
  * (20) Manukian, H., Pei, Y. R., Bearden, S. R. & Di Ventra, M.  Mode-assisted unsupervised learning of restricted boltzmann machines.  _Communications Physics_ 3, 1–8 (2020).
  * (21) Fischer, A. & Igel, C.  An introduction to restricted boltzmann machines.  In Alvarez, L., Mejail, M., Gomez, L. & Jacobo, J. (eds.) _Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications_ , 14–36 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).
  * (22) Fischer, A. & Igel, C.  Empirical analysis of the divergence of gibbs sampling based learning algorithms for restricted boltzmann machines.  In _International conference on artificial neural networks_ , 208–217 (Springer, 2010).
  * (23) Desjardins, G., Courville, A. & Bengio, Y.  On training deep boltzmann machines.  _arXiv preprint arXiv:1203.4416_ (2012).
  * (24) Bender, C. M. & Orszag, S. A.  _Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory_ (Springer Science & Business Media, 2013).
  * (25) Swendsen, R. H. & Wang, J.-S.  Replica monte carlo simulation of spin-glasses.  _Physical review letters_ 57, 2607 (1986).
  * (26) Sherrington, D. & Kirkpatrick, S.  Solvable model of a spin-glass.  _Physical review letters_ 35, 1792 (1975).
  * (27) Sheldon, F., Traversa, F. L. & Di Ventra, M.  Taming a nonconvex landscape with dynamical long-range order: Memcomputing ising benchmarks.  _Physical Review E_ 100, 053311 (2019).
  * (28) Neal, R. M. & Hinton, G. E.  A view of the em algorithm that justifies incremental, sparse, and other variants.  In _Learning in graphical models_ , 355–368 (Springer, 1998).
  * (29) Salakhutdinov, R. & Hinton, G.  An efficient learning procedure for deep boltzmann machines.  _Neural computation_ 24, 1967–2006 (2012).
  * (30) Salakhutdinov, R.  Learning and evaluating boltzmann machines (2008).


# Social Network Intelligence Analysis to Combat
Street Gang Violence

Damon Paulo, Bradley Fischl, Tanya Markow, Michael Martin, Paulo Shakarian
Network Science Center and
Dept. of Electrical Engineering
and Computer Science
U.S. Military Academy
West Point, NY 10996
$\\{$damon.paulo, bradley.fischl, tanya.markow, michael.martin1$\\}$
@usma.edu<EMAIL_ADDRESS>
###### Abstract

In this paper we introduce the Organization, Relationship, and Contact
Analyzer (ORCA) that is designed to aide intelligence analysis for law
enforcement operations against violent street gangs. ORCA is designed to
address several police analytical needs concerning street gangs using new
techniques in social network analysis. Specifically, it can determine “degree
of membership” for individuals who do not admit to membership in a street
gang, quickly identify sets of influential individuals (under the tipping
model), and identify criminal ecosystems by decomposing gangs into sub-groups.
We describe this software and the design decisions considered in building an
intelligence analysis tool created specifically for countering violent street
gangs as well as provide results based on conducting analysis on real-world
police data provided by a major American metropolitan police department who is
partnering with us and currently deploying this system for real-world use.

###### Index Terms:

complex networks; social networks; criminology

##  I Introduction

Violent street gangs are a major cause of criminal activity in the United
States [1]. In this paper, we present a new piece of software, Organizational,
Relationship, and Contact Analyzer (ORCA) that is designed from the ground-up
to apply new techniques in social network analysis and mining to support law
enforcement. In particular, we look to enable improved intelligence analysis
on criminal street gangs. The software combines techniques from logic
programming, [2] viral marketing, [3, 4] and community detection [5, 6] in a
usable application custom-tailored for law enforcement intelligence support.
This work is inspired by recent work in law enforcement that recognizes
similarities between gang members and insurgents and identifies adaptations
that can be made from current counter-insurgency (COIN) strategy to counter
gang violence. [1, 7] Due to the striking similarities between gang-violence
and COIN, the authors from the U.S. Military Academy (West Point) responded to
requests from a major metropolitan police department to transition recent
social network mining software. As a result, several West Point cadets were
able to not only conduct research, but also gain a better understanding of a
COIN environment.

The main contribution of this paper is the ORCA software which is the first
software, to our knowledge, that combine the aforementioned techniques into a
single piece of software designed for law-enforcement intelligence analysis.
The paper is organized as follows. Section II describes ORCA and its various
components while Section III provides the results of our preliminary
evaluation. Finally, related work is discussed in Section IV.

##  II System Design and Implementation

The police department we worked with to develop ORCA described several issues
concerning the intelligence analysis of street gangs. They desired a software
system that could accomplish the following tasks.

  1. 1.

Ability to ingest police arrest data and visualize network representations of
such data. The police data in question primarily consists of arrest reports
which include the individual’s personal information as well as claimed gang
membership (if disclosed). This data also includes relationships among
individuals arrested together.

  2. 2.

Ability to determine degree of group membership. While many gang members will
disclose their gang affiliation, some will not - likely fearing legal
consequences. Hence, to better allocate police efforts and intelligence
gathering resources, it is important to assign these unaffiliated members to a
gang (with a degree of confidence).

  3. 3.

Ability to identify sets of influential members. Though criminal street gangs
are decentralized, it is suspected that there are groups of individuals that
can exert influence throughout a given gang - encouraging criminal activity
that is more violent and risky than the norm. Identifying groups of
individuals in this position of influence would provide law enforcement
professionals insight into the ability of these organizations to easily adopt
such behavior.

  4. 4.

Ability to map the “ecosystem” of a given gang. Criminal street gangs tend to
be highly modular organizations, with identifiable sub-groups. In particular,
“corner crews” - groups of individuals in a gang who conduct illicit drug
transactions on the same street corner - will form highly connected clusters
in a social network representation of a street gang. Further, understanding
the relationships among these sub-organizations - both within a given gang and
between different gangs provides substantial insight into inter/intra gang
dynamics as well as in identifying individuals that connect different
organizations.

Note that all of the above described functions were included in ORCA in a
direct response to the needs of the law enforcement professionals that we have
met. These individuals are directly from our target user population to whom we
are transitioning the software to in mid-2013.

ORCA was written in Python 2.7.3 on top of the NetworkX 1.7 library
111http://networkx.github.io/documentation/latest/index.html. In Figure 1 we
show our overall scheme for the design of ORCA to address the above
challenges. The key components we leveraged in building the system included
the MANCaLog framework which we used to determine degree of membership,
TIP_DECOMP was leveraged to find sets of influential members, and the Louvain
algorithm which we built on to map the ecosystems of the various gangs. We
describe how we designed each major part throughout the remainder of this
section.

Figure 1: Overview of components and functions of ORCA.

Police data was ingested into the system in the form of spreadsheet data that
was derived from SQL queries from a police database. For purposes of
development and experimental evaluation, all personally-identifying
information was anonymized. From this data, we utilized the “individual
record” number of each person arrested and created a social network between
two individuals that were arrested together (co-arrestees). To handle this
network data structure, we utilized the Python library NetworkX mentioned
earlier. The police also required visualizations of these social networks
(which are shown later). These visualizations were also included in the
context of intelligence reports automatically generated by the software (which
are based on the subsequently described components).

###  II-A Determining Degree of Membership

Analytically associating arrested individuals with a criminal gang is an
important piece of intelligence to law enforcement officials. Even though
analysis itself is not direct evidence against an individual, it may be used
as a starting point to gather further information. Individuals who do not
admit to being in a gang may still have associates that do admit to membership
of a certain gang. Hence, we looked to contribute to the solution to this
problem by assigning such individuals a degree of membership \- a real number
in the interval $[0,1]$ that represents the confidence that the individuals is
in a given gang. A value of $0$ may be interpreted as having no information on
the individual’s affiliation with a given gang while a value of $1$ would be
interpreted as the individual as having admitted to being in the gang. To
assign these degree of membership values, we utilize MANCaLog, a logic-
programming based framework [2]. This framework considers a network
represented as a graph $G=(V,E)$ where $V$ is a set of vertices and $E$ is a
set of (undirected) relationships. For a given $i\in V$,
$\eta_{i}=\\{j|(i,j)\in E\\}$ and $d_{i}=|\eta_{i}|$. The vertices can be
assigned with a label from the set $\mathcal{L}$. Each label assigned to a
vertex is given a value (in the interval $[0,1]$ which associates that label
with the vertex with a degree of confidence. Hence, we use of network
representation of the co-arrestees as the graph and the set of different gangs
are the labels. Individuals who are known to be in a certain gang are assigned
a confidence value of $1$ for that gang and $0$ for the others. Through
MANCaLog, we can then assign a degree of confidence to the remaining nodes
derived from rules of the following form:

$\mathit{grp}_{1}\leftarrow\bigvee_{i}\neg\langle\mathit{grp}_{i},1\rangle,(\mathit{grp}_{1})_{\mathit{ifl}}$
(1)

Intuitively, the above rule says that a node will be assigned to group 1 if it
has not already been assigned to another group (with a confidence of $1$) and
has a certain number of neighbors who are also in group 1. A set of such rules
is referred to as a MANCaLog program and denoted $P$. A program may also have
facts of the form $(\langle\mathit{grp}_{i},x\rangle,v)$ which means that
vertex $v$ has a degree of membership in group $i$ with a confidence of $x$.

In the above rule, the confidence value assigned to a node being in group 1 is
based on the influence function ($\mathit{ifl}$) which maps natural numbers to
reals in the interval $[0,1]$. For instance, an influence function may look
something like the following:

$\textit{ifl}(x)=\begin{cases}0.0&\textit{if }x=0\\\ 0.1&\textit{if }1\leq
x\leq 3\\\ 0.5&\textit{if }x\geq 4\end{cases}$ (2)

In [2], the authors showed that such confidence values can be assigned to all
vertices in the network for all labels in polynomial time. However, a key
issue is to derive the influence function for the rules. We devised a simple
strategy for learning the influence function for each label (in this case
group) and have included it as algorithm IFL_LEARN.

Algorithm 1 IFL_LEARN

0: Group label $g$, program $P$

0:

1: Let $R$ be an array indexed from $1$ to $d_{*}$ (the maximum degree of the
network).

2: for $i=1,\ldots d_{*}$ do

3: Set $pos=0$, $tot=0$

4: for $v\in V$ do

5: Set $X_{v}=\\{v^{\prime}|(v^{\prime},v)\in E\wedge(\langle
g,1\rangle,v^{\prime})\in P\\}$

6: if $|X_{v}|\geq i$ then

7: $tot=tot+1$

8: if $(\langle g,1\rangle,v)\in P$ then

9: $pos=pos+1$

10: end if

11: end if

12: end for

13: If $i=0$, $R[i]=pos/tot-1.96*SER(pos,tot)$, else
$R[i]=\max(R[i-1],pos/(pos+neg)-1.96*SER(pos,tot))$ where $SER$ is the
standard error on the fraction of positive neighbors.

14: end for

15: return $R$

Algorithm IFL_LEARN considers one group label ($g$) and finds the fraction of
neighbors. We view each node as getting “signals” from its neighbors who are
in that group. For each potential number of signals (which is between one and
the maximum in-degree of the network) Algorithm 1 below computes the fraction
of nodes with at least that number of signals who are in the group. It then
computes the $95\%$ lower bound confidence interval based on standard error.
This is used as the lower bound of the degree of membership as returned by the
function. Then, based on an assumption of monotonicity (based on the results
of [8]), the algorithm sets the lower bound to be the maximum between the
previous lower bound and the current.

The algorithm returns a matrix $R$, and the resulting influence function is
$\mathit{ifl}(x)=R[x]$. This algorithm is a simple approach to learning rules
from available data and no novelty is claimed; more complex and general
approaches to rule learning are outside the scope of this paper and are the
topic of ongoing work.

###  II-B Identifying Seed Sets

According to the law enforcement professionals we met, the idea of influence
is critical to understanding the behavior of violent street gangs. Of
particular concern is the influence of radicalizing gang members - charismatic
individuals that not only participate in abnormally risky and violent
behavior, but have the ability to encourage others to do the same. In order to
identify sets of individuals who conduct such behavior, we consider the idea
of influence spread through the “tipping model” [9]. In such a model, we again
consider a population structured as a social network ($G=(V,E)$) where each
vertex $i\in V$ is adjacent to $d_{i}$ edges. We model a social contagion that
spreads through the network as follows: consider unaffected node $i$. If at
least $\lceil d_{i}/2\rceil$ of $i$’s neighbors have the contagion, then $i$
has the contagion – otherwise it does not. A key question regarding this model
is to identify a seed set \- a set of vertices in the network that, if
initially infected, will lead to the entire population receiving the
contagion. Identifying a seed set provides the analyst a set of individuals
that, when taken together, are very influential to the overall network.
Further, if the seed set is of minimal size, then the size of the set can be
used as a proxy to measure how easily the network can be influenced.
Unfortunately, a simple reduction from set cover shows that finding a seed set
of minimal size is an NP-hard problem. [3] However, we have designed a fast
heuristic algorithm that finds seed sets of very small size in practice. [4]
This algorithm is based on the idea of shell decomposition often cited in
physics literature [10, 11, 12, 13] but modified to ensure that the resulting
set will lead to all nodes being infected. The algorithm, TIP_DECOMP is
presented in this section.

Algorithm 2 TIP_DECOMP

0: Threshold function, $\theta$ and directed social network $G=(V,E)$

0: $V^{\prime}$

1: For each vertex $i,\ dist_{i}=\lfloor d_{i}/2\rfloor$.

2: FLAG = TRUE.

3: while FLAG do

4: Let $i$ be the element of $V$ where $dist_{i}$ is minimal.

5: if $dist_{i}=\infty$ then

6: FLAG = FALSE.

7: else

8: Remove $i$ from $G$ and for each $j$ in $\eta_{i}$, if $dist_{j}>0$, set
$dist_{j}=dist_{j}-1$. Otherwise set $dist_{j}=\infty$.

9: end if

10: end while

11: return All nodes left in $G$.

Intuitively, the algorithm proceeds as follows (Figure 1). Given network
$G=(V,E)$,at each iteration, pick the node for which $\lfloor d_{i}/2\rfloor$
is the least but positive (or $0$) and remove it. Once there are no nodes for
which this quantity is positive (or $0$), the algorithm outputs the remaining
nodes in the network.

Figure 2: Example of TIP_DECOMP for a simple network depicted in box A. Next
to each node label (lower-case letter) is the value for $\lfloor
d_{i}/2\rfloor$. In the first four iterations, nodes e, f, h, and i are
removed resulting in the network in box B. This is followed by the removal of
node j resulting in the network in box C. In the next two iterations, nodes a
and b are removed (boxes D-E respectively). Finally, node c is removed (box
F). The nodes of the final network, consisting of d and g, have negative
values for $\lfloor d_{i}/2\rfloor$ and become the output of the algorithm.

In addition to providing the seed set for each street gang, ORCA also provides
information about influential individuals. The shell number for each vertex
based on the process of shell decomposition [10] has been previously shown to
correlate with the vertex’s influence based on the networked version of
various epidemic models. [12]

###  II-C Identifying Ecosystems

Another important capability for ORCA was to decompose street gangs into
component organizations. In some cases, street gangs will not be monolithic
organizations, but rather confederations of various factions - many of which
may be ill identified by authorities. Further, to maintain a gang’s lines of
funding, sub-organizations known as “corner crews” operate as an informal unit
to conduct narcotics sales within certain parts of a given gang’s territory.

A common method to identify communities in a social network is to partition it
in a way to maximize a quantity known as modularity. [5] We shall use the
notation $C=\\{c_{1},\ldots,c_{q}\\}$ to denote a partition over set $V$ where
each $c_{i}\in C$ is a subset of $V$, for any $c_{i},c_{j}\in C$, $c_{i}\cap
c_{j}=\emptyset$ and $\bigcup_{i}c_{i}=V$. For a given partition, $C$, the
modularity $M(C)$ is a number in $[-1,1]$ . The modularity of a network
partition can be used to measure the quality of its community structure.
Originally introduced by Newman and Girvan [5], this metric measures the
density of edges within partitions compared to the density of edges between
partitions. A formal definition of this modularity for an undirected network
is defined as follows.

Given partition $C=\\{c_{1},\ldots,c_{q}\\}$, modularity,

$M(C)=\frac{1}{2m}\sum_{c\in C}\sum_{i,j\in c}w_{ij}-P_{ij}$

where $P_{ij}=\frac{k_{i}k_{j}}{2m}$.

The modularity of an optimal network partition can be used to measure the
quality of its community structure. Though modularity-maximization is NP-hard,
the approximation algorithm of Blondel et al. [6] (a.k.a. the “Louvain
algorithm”) has been shown to produce near-optimal partitions.222Louvain
modularity was computed using the implementation available from CRANS at
http://perso.crans.org/aynaud/communities/. The modularity associated with
this algorithm is often called the “Louvain modularity” and the associated
partition is a “Louvain partition.”

ORCA not only finds the Louvain partition, but also explores the relationships
among the sub-groups within and between gangs. For a given gang, ORCA
generates a graph showing the relationship among all sub-groups in that gang
and the neighboring sub-groups from different gangs. We call this the
“ecosystem” of a given gang. The ecosystem is of particular importance to law-
enforcement officers for many reasons. One, in particular, is the issue of
gang retaliation. Consider the following: in the aftermath of a violent
incident initiated by group A against group B, police enforcement will
increase patrols in the territory controlled by group B. As a result, group B
will rely on an allied organization (group C) to conduct retaliatory action
against group A. Identifying the ecosystem of a given gang helps provide
insight into organizations likely to conduct such activities.

Further, it identifies individuals who have connections to various other sub-
groups. The intuition behind these individuals is that they connect various
organizations together. We refer to these individuals as “connectors.”
Connectors are also important to law-enforcement personnel as they are
individuals who often connect geographic disparate sub-groups of a larger gang
organization. Another use for connectors is as a liaison between gangs that
cooperate. For instance, suppose group 1 wants to sell drugs in group 2’s
territory. An agreement is reached between the two groups that group 1 can
conduct these sales provided it pays a tax to group 2. Hence, a liaison
between the two groups may be a facilitator in such an arrangement.

##  III Evaluation

We evaluated ORCA on a police dataset of $5418$ arrests from a single police
district over a three period of time. There were $11,421$ relationships among
the arrests. From this data, ORCA assembled a social network consisting of
$1468$ individuals (who were members in one of $18$ gangs) and $1913$
relationships. ORCA was able to complete this assembly in addition to all
analysis (determining degree of membership, finding seed sets, and developing
ecosystems) took $34.3$ seconds on a commodity laptop (Windows 8, B960 2.2 GHz
processor with 4 GB RAM).

###  III-A Degree of Membership

First, we examined the performance of the system in identifying degree of
membership. As stated earlier, this was determined using our MANCaLog
framework that operated in two steps: learning the influence functions and
then applying the MANCaLog inference engine to compute degree of membership.
Figure 3 shows a plot of the influence functions derived for five of the
gangs. Specifically, it shows the number of neighboring individuals in the
given gang on the x-axis compared to the computed degree of membership for an
individual having that number of gang contacts on the y-axis. We note that our
degree of membership computation provides a similar result to [8] which
studies the related problem of social contagion.

In the second step required to determine the degree of membership, the
MANCaLog inference engine created logic programs that utilized the influence
functions to determine degree of membership for individuals who did not admit
to being in a gang. All of the $180$ unadmitted individuals connected to the
derived social network (individuals who did not admit to being in a gang but
were arrested with at least one other individual in the district) were
assigned a non-zero degree of membership based on learned rules. The majority
of these individuals could be assigned a degree of membership greater than
$0.5$ for at least one gang (see Figure 4). We note that many of these
individuals were assigned a degree of membership to multiple gangs. This may
primarily be due to the fact that the included arrest reports spanned a three-
year time-frame - which our police counterparts state is a relatively long
period of time in street gang culture. In such a time, gang members often
change allegiances and/or form new gangs. Exploring a longitudinal study where
time is explicitly considered is an important direction for future work. The
MANCaLog framework allows for temporal reasoning as well.

Figure 3: Influence functions for five of the gangs examined by ORCA. Figure
4: Histogram showing the number of individuals assigned a degree of membership
within a certain range.

###  III-B Seed Set Identification

We also examined the seed sets found by ORCA for each of the street gangs.
These are sets of individuals who can initiate a social cascade (under the
majority threshold tipping model) that will cause universal adoption in the
gang. The size of the seed set as a percentage of the population of each gang
is shown in Figure 5. The street gangs in the police district we examined were
racially segregated - belonging to one of two races. Police officers working
in the district have told us that gangs of Racial Group A are known for a more
centralized organizational structure while gangs of Racial Group B have
adopted a decentralized model. These groups are denoted in Figure 4. We were
able to quantify this observation as the gangs in Racial Group A had (on
average) seed sets $3.86\%$ smaller than those in Group B.

Figure 5: Seed size as a percentage of the total gang membership for the $18$
street gangs analyzed by ORCA organized into two different racial groups.

###  III-C Finding Communities and Identifying Ecosystems

As stated earlier, law enforcement personnel have much interest in identifying
sub-organizations of a given street gang. ORCA tackled this problem using the
Louvain algorithm as described in the previous section. For the $18$ street
gangs examined in this study, ORCA identified subgroups for each in an attempt
to optimize modularity. As modularity approaches $1.0$, the communities
produced by the algorithm become more segregated. We show the modularity
values in Figure 6. Again, we also noticed a difference in the organizational
structure based on the gang’s race (aligning with anecdotal police
observations). Gangs associated with Racial Group A (more centralized) had
lower modularity scores (by an average of $11.2\%$) than gangs associated with
Racial Group B (decentralized). This also corresponds with police observations
that gangs affiliated with Racial Group B tend to operate in smaller factions
(where the gang organizations server more as a confederation) as opposed to
gangs affiliated with Group A (which are typically more hierarchical).

Figure 6: The modularity of the partition found with the Louvain algorithm for
each gang.

The creation of a gang’s ecosystem (as described earlier) is derived directly
from the result of the Louvain partition. Sub-groups identified with the
partition are connected together based on social links between members or by
members who claim to be in more than one gang. Such an ecosystem is shown in
Figure 7. ORCA also provides an analytical report of an ecosystem listing all
sub-group relations and the strength of the relations (based on number of
social ties between organizations) (see Figure 8). Additionally, ORCA also
could identify individuals that connected various organizations. Example
output of this feature is shown in Figure 9.

Figure 7: An example ecosystem for one of the gangs analyzed by ORCA. Figure
8: An example ecosystem analysis for one of the gangs analyzed by ORCA. Figure
9: An example report of connectors within one of the gangs analyzed by ORCA.

###  III-D User Interface

Another priority in developing ORCA was to make the software accessible to a
variety of police personnel. We created a user interface using the TKinter
library version 8.5333http://wiki.python.org/moin/TkInter and provided network
visualizations using Matplotlib 1.2.0444http://matplotlib.org/. As displayed
in Figure 10, the users can utilize the interface to easily run new analysis
using arrest and relationship data taken from the result of a database query.
They can also view the reports generated by completed analyisis. When the
analysis is complete, the law enforcement personnel have the ability to view a
full report that is automatically generated in a portable document format
(PDF) using version 1.7 of the PyFPDF
library555https://code.google.com/p/pyfpdf/. This is useful because it is a
widely used format that is familiar to most users. In addition to the full
report, the users can also view specific reports or network visualizations
which allows them to easily view the results of the analysis that has been
conducted. Finally, we would like to note that this is a draft interface that
we will refine throughout the summer of 2013 by conducting studies in the
field as law enforcement personnel utilize the software. We plan to
incorporate the feedback that we gather into further versions of the software
so that it becomes a more useful and usable tool.

Figure 10: The ORCA user interface with network visualization and PDF report
output.

##  IV Related Work

There have been several other pieces of software introduced for both general
purpose network analysis as well as for law enforcement. The main contribution
of the software presented in this paper if the ORCA software. ORCA is that it
is tailor made for specific police requirements dealing with violent street
gangs. Further, it leverage new articiial intellignece and data mining
techniques such as MANCALog and TIP_DECOMP. Although ORCA is the first attempt
to combine these techniques in a piece of intelligence analysis software,
there have been other efforts to build software designed for social network-
based intelligence analysis that we describe below.

Previous software designed to support law-enforcement through social network
analysis include CrimeFighter [14] and CrimeLink [15]. However, these tools
provide complementary capabilities to ORCA. While ORCA is primarily designed
for intelligence analysis with respect to street gangs, CrimeLink is designed
to support investigations while CrimeFighter is directed more toward targeting
of individuals in criminal organizations - rather than understanding
membership in and linkages between sub-groups. Further, neither of these
software packages support degree of membership, seed set identification, or
the creation of ecosystems.

General-purpose social network analysis software also exist, including the
Organization Risk Analyzer (ORA) [16], and NodeXL (an add-on to Microsoft
Excel) [17]. While these tools are very powerful and contain many features,
they do not provide the degree of membership calculation or the identification
of seed sets as ORCA does. ORCA is also purposely designed for police
intelligence analysts to better understand street gangs, and features such as
community finding in ORCA are designed for with this application in mind
(hence the creation of ecosystems and identification of inter-group
connectors). This nuanced use of community detection for street gang
intelligence analysis is not ready for use “out of the box” in a general-
purpose social network analysis package such as ORA or NodeXL.

##  V Conclusion

In this paper we introduced ORCA - the Organizational, Relationship, and
Contact Analyzer - a tool designed from the ground-up to support intelligence
analysis to aide law enforcement personnel in combating violent street gangs.
We have shown how ORCA can meet various police analysis needs to include
determining degree for gang members of unknown affiliation, identifying sets
of influential individuals in a gang, and finding and analyzing the sub-
organizations of a gang to determine inter/intra-group relationships.

Our next step with regard to this work is to integrate geospatial and temporal
elements in the analysis - particularly with respect to community finding and
degree of membership. Currently we are working closely with a major
metropolitan police department to transition ORCA for use in law enforcement.
Throughout the summer of 2013, we are sending project assistants to work
closely with the police in order to identify additional police requirements
that can be addressed with techniques similar to those discussed in this
paper. Currently the police are employing this analysis for one district.
There are plans to expand to other districts in late 2013.

## Acknowledgment

The authors are supported under by the Army Research Office (project
2GDATXR042) and the Office of the Secretary of Defense. The opinions in this
paper are those of the authors and do not necessarily reflect the opinions of
the funders, the U.S. Military Academy, or the U.S. Army.

## References

  * [1] J. Bertetto, “Countering criminal street gangs: Lessons from the counterinsurgent battlespace,” _Law Enforcement Executive Forum_ , vol. 12, no. 3, p. 43, 2012.
  * [2] P. Shakarian, G. I. Simari, and R. Schroeder, “MANCaLog: A logic for multi-attribute network cascades,” in _Proc. of 26th Intl. Conf. on Autonomous Agents and Multiagent Systems (AAMAS)_ , May 2013.
  * [3] D. Kempe, J. Kleinberg, and E. Tardos, “Maximizing the spread of influence through a social network,” in _KDD ’03: Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining_. New York, NY, USA: ACM, 2003, pp. 137–146.
  * [4] P. Shakarian and D. Paulo, “Large social networks can be targeted for viral marketing with small seed sets,” in _Proc. 2012 IEEE/ACM Intl. Conf. on Advances in Social Networks Analysis and Mining (ASONAM-12)_ , Aug. 2012.
  * [5] M. E. J. Newman and M. Girvan, “Finding and evaluating community structure in networks,” _Phys. Rev. E_ , vol. 69, no. 2, p. 026113, Feb 2004.
  * [6] V. Blondel, J. Guillaume, R. Lambiotte, and E. Lefebvre, “Fast unfolding of communities in large networks,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2008, p. P10008, 2008.
  * [7] E. Goode, “Combating gang warfare with green beret tactics,” _The New York Times_ , Apr. 2012.
  * [8] D. Centola, “The Spread of Behavior in an Online Social Network Experiment,” _Science_ , vol. 329, no. 5996, pp. 1194–1197, Sep. 2010.
  * [9] M. Granovetter, “Threshold models of collective behavior,” _The American Journal of Sociology_ , no. 6, pp. 1420–1443.
  * [10] S. B. Seidman, “Network structure and minimum degree,” _Social Networks_ , vol. 5, no. 3, pp. 269 – 287, 1983.
  * [11] S. Carmi, S. Havlin, S. Kirkpatrick, Y. Shavitt, and E. Shir, “From the Cover: A model of Internet topology using k-shell decomposition,” _PNAS_ , vol. 104, no. 27, pp. 11 150–11 154, 2007. [Online]. Available: http://www.pnas.org/cgi/content/abstract/104/27/11150
  * [12] M. Kitsak, L. K. Gallos, S. Havlin, F. Liljeros, L. Muchnik, H. E. Stanley, and H. A. Makse, “Identification of influential spreaders in complex networks,” _Nat Phys_ , vol. 6, no. 11, pp. 888–893, Nov. 2010.

  * [13] G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, “Heterogeneous $k$-core versus bootstrap percolation on complex networks,” _Phys. Rev. E_ , vol. 83, May 2011.
  * [14] R. Petersen, C. Rhodes, and U. Wiil, “Node removal in criminal networks,” in _Intelligence and Security Informatics Conference (EISIC), 2011 European_ , Sep. 2011, pp. 360 –365.
  * [15] J. Schroeder, J. Xu, and H. Chen, “Crimelink explorer: Using domain knowledge to facilitate automated crime association analysis,” in _Intelligence and Security Informatics_ , ser. Lecture Notes in Computer Science, H. Chen, R. Miranda, D. Zeng, C. Demchak, J. Schroeder, and T. Madhusudan, Eds. Springer Berlin Heidelberg, 2003, vol. 2665, pp. 168–180.
  * [16] K. M. Carley, “Ora: Organization risk analyzer,” CMU CASOS, Tech. Rep. January, 2004.
  * [17] D. Hansen, B. Shneiderman, and M. A. Smith, _Analyzing Social Media Networks with NodeXL: Insights from a Connected World_. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 2010\.


# An online bin-packing problem with an underlying ternary structure

Helmut Prodinger  Helmut Prodinger, Mathematics Department, Stellenbosch
University, 7602 Stellenbosch, South Africa<EMAIL_ADDRESS>
###### Abstract.

Following an orginal idea by Knödel, an online bin-packing problem is
considered where the the large items arrive in double-packs. The dual problem
where the small items arrive in double-packs is also considered. The
enumerations have a ternary random walk flavour, and for the enumeration, the
kernel method is employed.

###### Key words and phrases:

Knödel walks, third-order recursion, kernel method, coefficient extraction,
state diagram

###### 2010 Mathematics Subject Classification:

05A15, 68R05

##  1\. Introduction

Walter Knödel introduced the following online bin-packing problem [3]: There
are bins of size 1, and random items of size $\frac{2}{3}$ (large items) and
of size $\frac{1}{3}$ (small items) appear and are put into the boxes. A
typical scenario is that a number $j$ of partially filled boxes exist, and the
number $j$ becomes $j+1$ resp. $j-1$, depending on whether a the new item is
of large resp. small type. “At random” means that both types appear with the
same probability $\frac{1}{2}$.

In my collection of examples [4], I showed how to deal with the Knödel problem
using the kernel method. I was, however, not the only author who was intrigued
by such questions; a notable paper is by Michael Drmota [1], which is of a
more probabilistic type, whereas I tried to emphasize the combinatorial point
of view.

The present paper has a certain ‘ternary’ flavour: the next section deals with
the instance of large items appearing in double-packs. The handler breaks off
the double-packs, and then treats the items as Knödel would have done.
Typically, the number of partially filled boxes increases by 2 or decreases by
1. In order to keep the system balanced, we assume that the small items appear
twice as often as the double-packs.

The last section deals with the dual problem, where the small items appear in
double-packs and the large items as single units.

The kernel method is used to obtain all the relevant enumerations. The recent
paper [5] served as an inspiration, but deals with a different issue. It must
be said that, when [4] was prepared, such ternary questions would have been
outside of my reach. Luckily, now, they are not.

We confine ourselves here just to enumerations, deriving explicit generating
functions in one or two variables. Questions of a more probabilistic nature
are not treated.

##  2\. The first model

The following items arrive at random: a double-pack of items, each of size
$\frac{2}{3}$, and an item of size $\frac{1}{3}$. We could equip the set-up
with general probabilities $p$ and $q=1-p$, but we restrict ourselves to the
‘balanced’ case where the single items are twice as likely as the double-
packs, so we set $p=\frac{1}{3}$ and $q=\frac{2}{3}$.

The following state diagram (we show only a finite part of it) describes the
situation. There are states representing ‘$i$ boxes filled to $\frac{2}{3}$’;
a double-pack pushes the $i$ to $i+2$, and a single item reduces it to $i-1$.
There is an exceptional state, called $\beta$, standing for one box, filled to
$\frac{1}{3}$. The red edges represent an arrival of a double-pack, and will
be labelled by $pz$; the black edges represent an arrival of a single item,
and will be labelled by $qz$.

$0$$1$$2$$3$$4$$5$$6$$7$$8$$\beta$

From the state diagram, we set off an infinite set of generating functions in
the variable $z$, where the coefficient of $z^{n}$ is the probability that $n$
random steps lead to state $i$, for $i\geq 0$ or $i=\beta$. Mostly, we just
write $f_{i}$ instead of $f_{i}(z)$. The following system of recursions can be
read off immediately:

$\displaystyle f_{0}$ $\displaystyle=1+qzf_{1},\quad f_{\beta}=qzf_{0},$
$\displaystyle f_{1}$ $\displaystyle=zf_{\beta}+qzf_{2}=qz^{2}f_{0}+qzf_{2},$
$\displaystyle f_{i}$ $\displaystyle=pzf_{i-2}+qzf_{i+1},\quad i\geq 2.$

Our method to solve this system is the kernel method. For that, we introduce a
bivariate generating function $F(u,z)$, but we mostly write just $F(u)$:

$\displaystyle F(u)$ $\displaystyle=\sum_{i\geq 0}u^{i}f_{i}(z)$
$\displaystyle=1+qzf_{1}+qz^{2}uf_{0}+qzuf_{2}+\sum_{i\geq
2}u^{i}\Big{[}pzf_{i-2}+qzf_{i+1}\Big{]}$
$\displaystyle=1+qzf_{1}+qz^{2}uf_{0}+pzu^{2}F(u)+\sum_{i\geq
1}u^{i}qzf_{i+1}$
$\displaystyle=1+qzf_{1}+qz^{2}uf_{0}+pzu^{2}F(u)+\frac{qz}{u}\Big{(}F(u)-f_{0}-uf_{1}\Big{)}$
$\displaystyle=1+qz^{2}uf_{0}+pzu^{2}F(u)+\frac{qz}{u}\Big{(}F(u)-f_{0}\Big{)}$

Note that $f_{0}=F(0)$. It is beneficial to introduce the new variable $u=zU$;
doing this, powers of $z$ that appear are multiples of $3$. As can be seen,
the numbers of steps leading to a state $i$ belong to just one residue class
modulo 3. We compute

$F(u)=\frac{-3U-2z^{3}U^{2}f_{0}+2f_{0}}{z^{3}U^{3}-3U+2}=\frac{-3U-2xU^{2}f_{0}+2f_{0}}{xU^{3}-3U+2}.$

As it is common using the kernel method, setting $U=0$ leads to a void
equation. However, factorizing the denominator is the method of choice. There
is ‘bad’ factor in the denominator, which must also appear in the numerator,
which allows us to compute $f_{0}$ and consequently the whole bivariate
generating function. In order to deal with the ternary equation successfully,
we further set $x=z^{3}=\frac{27}{4}t(1-t)^{2}$ and we find the 3 roots

$\displaystyle U_{1}$ $\displaystyle=\frac{2}{3(1-t)},\quad
U_{2}=\frac{1}{\sigma},\quad U_{3}=\frac{1}{\tau},$

with

$\displaystyle\sigma$
$\displaystyle=\frac{3}{4}(t-\sqrt{4t-3t^{2}}\,),\quad\tau=\frac{3}{4}(t+\sqrt{4t-3t^{2}}).$

Plugging $U=\frac{2}{3(1-t)}$ into the numerator (this is the bad factor, as
explained a little bit later), leads to

$f_{0}=\frac{1}{(1-t)(1-3t)}$

and furthermore to the simplified numerator

$\frac{-3U-2xU^{2}f_{0}+2f_{0}}{U-\frac{2}{3(1-t)}}=\frac{1}{1-3t}\Big{(}-3+\frac{27}{2}t(t-1)U\Big{)}.$

The variable $x$ is given in terms of $t$. The inverse relation is of
interest. It can be obtained by the Lagrange inversion formula or, as here, by
contour integration:

$\displaystyle[x^{k}]t$ $\displaystyle=\frac{1}{2\pi
i}\oint\frac{dx}{x^{k+1}}t=\frac{1}{2\pi
i}\frac{27}{4}\Big{(}\frac{4}{27}\Big{)}^{k+1}\oint\frac{dt(1-t)(1-3t)}{t^{k+1}(1-t)^{2k+2}}t$
$\displaystyle=\frac{1}{2\pi
i}\Big{(}\frac{4}{27}\Big{)}^{k}\oint\frac{dt(1-3t)}{t^{k}(1-t)^{2k+1}}=\Big{(}\frac{4}{27}\Big{)}^{k}[t^{k-1}]\frac{1-3t}{(1-t)^{2k+1}}$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{k}\bigg{[}\binom{3k-1}{k-1}-3\binom{3k-2}{k-2}\bigg{]},$

which, after simplification, gives us

$t=\sum_{k\geq 1}\frac{1}{k}\binom{3k-2}{k-1}\frac{2^{2k}}{3^{3k}}.$

A similar computation leads to

$\frac{1}{1-t}=\sum_{k\geq
0}\frac{1}{2k+1}\binom{3k}{k}\frac{2^{2k+1}}{3^{3k+1}}.$

From this we infer that $U\sim\frac{2}{3}$, or $u\sim\frac{2}{3}z$, explaining
why we are talking about the bad factor. We continue the computation:

$\displaystyle F(u)$
$\displaystyle=\frac{1}{1-3t}\Big{(}-3+\frac{27}{2}t(t-1)U\Big{)}\frac{1}{x(U-\frac{1}{\sigma})(U-\frac{1}{\tau})}$
$\displaystyle=\frac{1}{1-3t}\Big{(}-3+\frac{27}{2}t(t-1)U\Big{)}\frac{\frac{9}{4}t(t-1)}{x(1-\sigma
U)(1-\tau U)}$
$\displaystyle=\frac{1}{(1-3t)(1-t)}\Big{(}1-\frac{9}{2}t(t-1)U\Big{)}\frac{1}{(1-\sigma
U)(1-\tau U)}.$

Partial fraction decomposition leads to (we use the abbreviation
$W=\sqrt{4t-3t^{2}}\,$)

$\displaystyle\frac{1}{(1-\sigma U)(1-\tau U)}$
$\displaystyle=\frac{1}{2}\Big{(}1-\frac{t}{W}\Big{)}\frac{1}{1-\sigma
U}+\frac{1}{2}\Big{(}1+\frac{t}{W}\Big{)}\frac{1}{1-\tau U}$
$\displaystyle=\frac{1}{2}\Big{[}\frac{1}{1-\sigma U}+\frac{1}{1-\tau
U}\Big{]}+\frac{t}{2W}\Big{[}\frac{1}{1-\tau U}-\frac{1}{1-\sigma U}\Big{]}$
$\displaystyle=\frac{1}{2}\Big{[}\frac{1}{1-\sigma U}+\frac{1}{1-\tau
U}\Big{]}+\frac{3t}{4(\tau-\sigma)}\Big{[}\frac{1}{1-\tau U}-\frac{1}{1-\sigma
U}\Big{]}.$

For the further simplification we will resort to two identities going by the
name of Girard-Waring formula, see e. g. [2]:

$X^{m}+Y^{m}=\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}\frac{m}{m-k}(XY)^{k}(X+Y)^{m-2k};$
$\frac{X^{m}-Y^{m}}{X-Y}=\sum_{0\leq
k\leq(m-1)/2}(-1)^{k}\binom{m-1-k}{k}(XY)^{k}(X+Y)^{m-1-2k}.$

Of course, we will apply them with $X=\tau$ and $Y=\sigma$. Then

$\displaystyle[U^{m}]\frac{1}{2}\Big{[}\frac{1}{1-\sigma U}+\frac{1}{1-\tau
U}\Big{]}$ $\displaystyle=\frac{1}{2}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}\frac{m}{m-k}\Big{(}\frac{9}{4}t(t-1)\Big{)}^{k}\Big{(}\frac{3}{2}t\Big{)}^{m-2k}$
$\displaystyle=\frac{1}{2}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}\frac{m}{m-k}t^{k}(t-1)^{k}t^{m-2k}$
$\displaystyle=\frac{1}{2}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}\frac{m}{m-k}(t-1)^{k}t^{m-k}$

and

$\displaystyle[U^{m}]$
$\displaystyle\frac{3t}{4(\tau-\sigma)}\Big{[}\frac{1}{1-\tau
U}-\frac{1}{1-\sigma U}\Big{]}$ $\displaystyle=\frac{3t}{4}\sum_{0\leq
k\leq(m-1)/2}(-1)^{k}\binom{m-1-k}{k}\Big{(}\frac{9}{4}t(t-1)\Big{)}^{k}\Big{(}\frac{3}{2}t\Big{)}^{m-1-2k}$
$\displaystyle=\Big{(}\frac{3}{2}\Big{)}^{m-1}\frac{3t}{4}\sum_{0\leq
k\leq(m-1)/2}(-1)^{k}\binom{m-1-k}{k}(t-1)^{k}t^{m-1-k}$
$\displaystyle=\frac{1}{2}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq
k\leq(m-1)/2}(-1)^{k}\binom{m-1-k}{k}(t-1)^{k}t^{m-k}.$

Combining the two leads to a pleasant simplification:

$\displaystyle[U^{m}]$ $\displaystyle\frac{1}{2}\Big{[}\frac{1}{1-\sigma
U}+\frac{1}{1-\tau
U}\Big{]}+[U^{m}]\frac{3t}{4(\tau-\sigma)}\Big{[}\frac{1}{1-\tau
U}-\frac{1}{1-\sigma U}\Big{]}$
$\displaystyle=\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}(t-1)^{k}t^{m-k},$

or simpler

$\displaystyle[U^{m}]\frac{1}{(1-\sigma U)(1-\tau
U)}=\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}(t-1)^{k}t^{m-k}.$

We need a second similar term:

$\displaystyle[U^{m}]$
$\displaystyle\Big{(}-\frac{9}{2}t(t-1)U\Big{)}\frac{1}{(1-\sigma U)(1-\tau
U)}$ $\displaystyle=-\frac{9}{2}t(t-1)[U^{m-1}]\frac{1}{(1-\sigma U)(1-\tau
U)}$
$\displaystyle=-\frac{9}{2}t(t-1)\Big{(}\frac{3}{2}\Big{)}^{m-1}\sum_{0\leq
k\leq m/2}(-1)^{k}\binom{m-1-k}{k}(t-1)^{k}t^{m-1-k}$
$\displaystyle=-3\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-1-k}{k}(t-1)^{k+1}t^{m-k}.$

Putting everything together we found

$\displaystyle F(u)$ $\displaystyle=\frac{1}{(1-3t)(1-t)}\sum_{m\geq
0}\frac{u^{m}}{z^{m}}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}(t-1)^{k}t^{m-k}$
$\displaystyle-\frac{3}{(1-3t)(1-t)}\sum_{m\geq
0}\frac{u^{m}}{z^{m}}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-1-k}{k}(t-1)^{k+1}t^{m-k}.$

Reading off coefficients of powers of $z$ as well is now done with Cauchy’s
integral formula; the contours are always small circles (or equivalent) around
the origin. The starting point is

$\displaystyle[z^{n}u^{j}]F(u)$
$\displaystyle=\Big{(}\frac{3}{2}\Big{)}^{j}[z^{n+j}]\frac{1}{(1-3t)(1-t)}\sum_{0\leq
k\leq j/2}(-1)^{k}\binom{j-k}{k}(t-1)^{k}t^{j-k}$
$\displaystyle-\Big{(}\frac{3}{2}\Big{)}^{j}[z^{n+j}]\frac{3}{(1-3t)(1-t)}\sum_{0\leq
k\leq j/2}(-1)^{k}\binom{j-1-k}{k}(t-1)^{k+1}t^{j-k}$

and we will treat the two sums separately. There is only a contribution if
$n+j\equiv 0\bmod 3$. (This is also clear from the combinatorial context.)
Assume this and set $N:=\frac{n+j}{3}$.

Step 1:

$\displaystyle[x^{N}]$ $\displaystyle\frac{1}{(1-3t)(1-t)}\sum_{0\leq k\leq
j/2}(-1)^{k}\binom{j-k}{k}(t-1)^{k}t^{j-k}$ $\displaystyle=\frac{1}{2\pi
i}\oint\frac{dx}{x^{N+1}}\frac{1}{(1-3t)(1-t)}\sum_{0\leq k\leq
j/2}(-1)^{k}\binom{j-k}{k}(t-1)^{k}t^{j-k}$ $\displaystyle=\frac{1}{2\pi
i}\oint\frac{27}{4}\frac{dt}{\big{(}\frac{27}{4}t(1-t)^{2}\big{)}^{N+1}}\sum_{0\leq
k\leq j/2}(-1)^{k}\binom{j-k}{k}(t-1)^{k}t^{j-k}$ $\displaystyle=\frac{1}{2\pi
i}\oint\Big{(}\frac{4}{27}\Big{)}^{N}dt\sum_{0\leq k\leq
j/2}(-1)^{k}\binom{j-k}{k}(t-1)^{k-2N-2}t^{j-k-N-1}$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{N}[t^{N-j+k}]\sum_{0\leq k\leq
j/2}(-1)^{k}\binom{j-k}{k}(t-1)^{k-2N-2}$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{N}(-1)^{N-j}\sum_{0\leq k\leq
j/2}\binom{j-k}{k}\binom{k-2N-2}{N-j+k}.$

Step 2:

$\displaystyle[x^{N}]$ $\displaystyle\frac{3}{(1-3t)(1-t)}\sum_{0\leq k\leq
j/2}(-1)^{k}\binom{j-1-k}{k}(t-1)^{k+1}t^{j-k}$
$\displaystyle=3\Big{(}\frac{4}{27}\Big{)}^{N}[t^{N-j+k}]\sum_{0\leq k\leq
j/2}(-1)^{k}\binom{j-1-k}{k}(t-1)^{k-2N-1}$
$\displaystyle=3\Big{(}\frac{4}{27}\Big{)}^{N}(-1)^{N-j}\sum_{0\leq k\leq
j/2}\binom{j-1-k}{k}\binom{k-2N-1}{N-j+k}.$

We put all the results of this section together in a theorem.

######  Theorem 1.

The generating function $F(u)=F(u,z)$ has the following explicit form:

$F(u)=\frac{1}{(1-3t)(1-t)}\Big{(}1-\frac{9}{2}t(t-1)U\Big{)}\frac{1}{(1-\sigma
U)(1-\tau U)}$

Here, $u=zU$, $z^{3}=x=\frac{27}{4}t(1-t)^{2}$, and

$\displaystyle\sigma$
$\displaystyle=\frac{3}{4}(t-\sqrt{4t-3t^{2}}\,),\quad\tau=\frac{3}{4}(t+\sqrt{4t-3t^{2}}\,).$

Note that $(1-\sigma U)(1-\tau U)=1-\frac{3}{2}tU+\frac{9}{4}t(t-1)U^{2}$.
Written in the new variable $U$, only powers of $z$ that are multiples of $3$
appear. Further, we get the representation sorted by powers of $u$:

$\displaystyle F(u,z)$ $\displaystyle=\frac{1}{(1-3t)(1-t)}\sum_{m\geq
0}\frac{u^{m}}{z^{m}}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq k\leq
m/2}(-1)^{k}\binom{m-k}{k}(t-1)^{k}t^{m-k}$
$\displaystyle-\frac{3}{(1-3t)(1-t)}\sum_{m\geq
0}\frac{u^{m}}{z^{m}}\Big{(}\frac{3}{2}\Big{)}^{m}\sum_{0\leq
k\leq(m-1)/2}(-1)^{k}\binom{m-1-k}{k}(t-1)^{k+1}t^{m-k}.$

Reading off coefficients of $z^{N}u^{j}$, where $N=\frac{n+j}{3}$ leads to

$\displaystyle[z^{N}u^{j}]F(u,z)$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{N}(-1)^{N-j}\sum_{0\leq k\leq
j/2}\binom{j-k}{k}\binom{k-2N-2}{N-j+k}$
$\displaystyle-3\Big{(}\frac{4}{27}\Big{)}^{N}(-1)^{N-j}\sum_{0\leq
k\leq(j-1)/2}\binom{j-1-k}{k}\binom{k-2N-1}{N-j+k}.$

For the special state $\beta$, the following series representation holds:

$\displaystyle f_{\beta}(z)=\sum_{n\geq
0}\frac{2^{2n+1}}{{3}^{3n+1}}\binom{3n+1}{n}z^{3n+1}.$

∎

The computation for the special state was not shown yet:

$\displaystyle[z^{3n+1}]$ $\displaystyle
f_{\beta}=\frac{2}{3}[x^{n}]\frac{1}{(1-t)(1-3t)}=\frac{2}{3}\frac{1}{2\pi
i}\oint\frac{dx}{x^{n+1}}\frac{1}{(1-t)(1-3t)}$
$\displaystyle=\frac{2}{3}\frac{27}{4}\frac{1}{2\pi
i}\oint\frac{dt}{(\tfrac{27}{4})^{n+1}t^{n+1}(1-t)^{2n+2}}=\frac{2}{3}\Big{(}\frac{4}{27}\Big{)}^{n}\frac{1}{2\pi
i}\oint\frac{dt}{t^{n+1}(1-t)^{2n+2}}$
$\displaystyle=\frac{2}{3}\Big{(}\frac{4}{27}\Big{)}^{n}[t^{n}]\frac{1}{(1-t)^{2n+2}}=\frac{2^{2n+1}}{{3}^{3n+1}}\binom{3n+1}{n}.$

##  3\. The dual model

Now, the red edges mean the arrival of the large objects (size $\frac{2}{3}$)
and the black edges mean a double-pack of the small edges (size $\frac{1}{3}$
each). To keep the system balanced, the large objects should arrive twice as
often as the double-packs of small edges. Again, the generating function
$g_{i}$ refers to paths of length $n$ leading eventually into state $i$. After
$n$ steps, only a state $i$ can be reached with $n\equiv i\bmod 3$. The state
diagram and the recursions are immediate:

$0$$1$$2$$3$$4$$5$$6$$7$$8$$\beta$

We work only with $p=\frac{2}{3}$, $q=\frac{1}{3}$. Directly from the state
diagram,

$\displaystyle g_{0}$
$\displaystyle=1+zg_{\beta}+qzg_{2}=1+qz^{2}g_{1}+qzg_{2},$ $\displaystyle
g_{\beta}$ $\displaystyle=qzg_{1},\quad g_{1}=zg_{0}+qzg_{3},$ $\displaystyle
g_{i}$ $\displaystyle=pzg_{i-1}+qzg_{i+2},\quad i\geq 2.$

Summing the recursions,

$\displaystyle G(u)$ $\displaystyle=g_{0}+ug_{1}+\sum_{i\geq
2}u^{i}\Big{(}pzg_{i-1}+qzg_{i+2}\Big{)}$
$\displaystyle=g_{0}+uzg_{0}+qzug_{3}+pzu\sum_{i\geq
1}u^{i}g_{i}+\frac{qz}{u^{2}}\sum_{i\geq 4}u^{i}g_{i}$
$\displaystyle=g_{0}+uzg_{0}+pzuG(u)-pzug_{0}+\frac{qz}{u^{2}}\sum_{i\geq
3}u^{i}g_{i}$
$\displaystyle=g_{0}+uzg_{0}+pzuG(u)-pzug_{0}+\frac{qz}{u^{2}}(G(u)-g_{0}-ug_{1}-u^{2}g_{2})$
$\displaystyle=g_{0}+quzg_{0}+pzuG(u)+\frac{qz}{u^{2}}G(u)-\frac{qz}{u^{2}}g_{0}-\frac{qz}{u}g_{1}-qzg_{2}$
$\displaystyle=g_{0}+quzg_{0}+pzuG(u)+\frac{qz}{u^{2}}G(u)-\frac{qz}{u^{2}}g_{0}-\frac{qz}{u}g_{1}+1+qz^{2}g_{1}-g_{0}.$

Solving, we find with $V=uz$:

$G(u)=\frac{-{V}^{3}g_{0}-3{V}^{2}-g_{1}{V}^{2}{z}^{2}+{z}^{3}g_{0}+g_{1}V{z}^{2}}{2V^{3}-3V^{2}+x}.$

Now we factorize the denominator:

$2(V-\tfrac{3}{2}(1-t))(V-\sigma)(V-\tau)=2V^{3}-3V^{2}+x.$

This time, both, $(V-\sigma)$ and $(V-\tau)$ are bad factors. Plugging into
the numerator, we find two equations, and the solutions:

$g_{0}=\frac{4}{(1-3t)(4-3t)},\quad g_{1}=\frac{27t(1-t)}{z^{2}(1-3t)(4-3t)}.$

Dividing out the bad factors, we find

$\frac{-{V}^{3}g_{0}-3{V}^{2}-g_{1}{V}^{2}{z}^{2}+{z}^{3}g_{0}+g_{1}V{z}^{2}}{(V-\sigma)(V-\tau)}=\frac{12(t-1)-4V}{(1-3t)(4-3t)}.$

Altogether:

$\displaystyle G(u)$
$\displaystyle=\frac{6(t-1)-2V}{(1-3t)(4-3t)}\frac{1}{V-\frac{3}{2}(1-t)}$
$\displaystyle=\frac{6(1-t)+2V}{(1-3t)(4-3t)\frac{3}{2}(1-t)}\frac{1}{1-\frac{2}{3(1-t)}V}$
$\displaystyle=2\frac{2+\frac{2}{3(1-t)}V}{(1-3t)(4-3t)}\frac{1}{1-\frac{2}{3(1-t)}V}.$

Furthermore

$\displaystyle[V^{j}]G(u)$
$\displaystyle=\frac{2}{(1-3t)(4-3t)}\bigg{[}2\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}+\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}\bigg{]}$
$\displaystyle=\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j},\quad
j\geq 1,$

and

$\displaystyle[u^{j}]G(u)=z^{j}\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}.$

Now let us consider $j+3N$ steps to reach state $j$, and then

$\displaystyle[z^{j+3N}u^{j}]G(u)$
$\displaystyle=[x^{N}]\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}$
$\displaystyle=\frac{1}{2\pi
i}\oint\frac{dx}{x^{N+1}}\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}$
$\displaystyle=\frac{27}{4}\Big{(}\frac{4}{27}\Big{)}^{N+1}\frac{1}{2\pi
i}\oint\frac{dt(1-t)(1-3t)}{t^{N+1}(1-t)^{2N+2}}\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{N}\Big{(}\frac{2}{3}\Big{)}^{j}\frac{1}{2\pi
i}\oint\frac{dt}{t^{N+1}(1-t)^{2N+j+1}}\frac{6}{(4-3t)}$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{N}\Big{(}\frac{2}{3}\Big{)}^{j-1}[t^{N}]\frac{1}{(1-t)^{2N+j+1}}\frac{1}{(1-\frac{3}{4}t)}$
$\displaystyle=\Big{(}\frac{4}{27}\Big{)}^{N}\Big{(}\frac{2}{3}\Big{)}^{j-1}\sum_{i=0}^{N}\Big{(}\frac{3}{4}\Big{)}^{N-i}\binom{2N+j+i}{i}$
$\displaystyle=\sum_{i=0}^{N}\frac{2^{2i+j-1}}{3^{2N+i+j-1}}\binom{2N+j+i}{i}.$

The coefficients of $g_{0}$ are different:

$[z^{3N}]g_{0}=\sum_{i=0}^{N}\frac{2^{2i}}{3^{2N+i}}\binom{2N+i}{i}.$

Furthermore,

$[z^{3N+1}]g_{\beta}=\frac{1}{3}[z^{3N}]g_{1}=\sum_{i=0}^{N}\frac{2^{2i}}{3^{2N+i+1}}\binom{2N+1+i}{i}.$

Here are the main results of this section:

######  Theorem 2.

The generating function $G(u)=G(u,z)$ has the following explicit form:

$G(u)=2\frac{2+\frac{2}{3(1-t)}V}{(1-3t)(4-3t)}\frac{1}{1-\frac{2}{3(1-t)}V}.$

Here, $u=\frac{V}{z}$, $z^{3}=x=\frac{27}{4}t(1-t)^{2}$. Written in the new
variable $V$, only powers of $z$ that are multiples of $3$ appear. Further, we
get the representation sorted by powers of $u$:

$\displaystyle[V^{j}]G(u)$
$\displaystyle=\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j},\quad
j\geq 1,$

and

$\displaystyle[u^{j}]G(u)=z^{j}\frac{6}{(1-3t)(4-3t)}\Big{(}\frac{2}{3}\frac{1}{1-t}\Big{)}^{j}.$

Reading off coefficients of $z^{j+3N}u^{j}$ leads to

$\displaystyle[z^{j+3N}u^{j}]G(u,z)=\sum_{i=0}^{N}\frac{2^{2i+j-1}}{3^{2N+i+j-1}}\binom{2N+j+i}{i}.$

For the special cases, the following series representation holds:

$\displaystyle[z^{3N}]g_{0}$
$\displaystyle=\sum_{i=0}^{N}\frac{2^{2i}}{3^{2N+i}}\binom{2N+i}{i},$
$\displaystyle[z^{3N+1}]g_{\beta}$
$\displaystyle=\frac{1}{3}[z^{3N}]g_{1}=\sum_{i=0}^{N}\frac{2^{2i}}{3^{2N+i+1}}\binom{2N+1+i}{i}.$

∎

## References

  * [1] Michael Drmota.  Discrete Random Walks on One-Sided “Periodic” Graphs,  Discrete Mathematics and Theoretical Computer Science, January 1, 2003, DMTCS Proceedings vol. AC, Discrete Random Walks (DRW’03), 83–94.
  * [2] Henry Gould.  The Girard-Waring power sum formulas for symmetric functions and Fibonacci sequences.  The Fibonacci Quarterly, 37(2):135–140, 1999.
  * [3] Walter Knödel.  Über das mittlere Verhalten von online-Packungsalgorithmen,  Elektron. Informationsverarb. und Kybernetik (EIK), 19 (1983), 9, 427–433.
  * [4] Helmut Prodinger.  The Kernel Method: A Collection of Examples,  Séminaire Lotharingien de Combinatoire, B50f (2004), 19 pp.
  * [5] Helmut Prodinger.  Enumeration of S-Motzkin paths from left to right and from right to left — a kernel method approach,  PU. M. A. Vol. 29 (2020), No.1, pp. 28–38.