# CORE: Simple and Effective Session-based Recommendation within Consistent
Representation Space

Yupeng Hou†<EMAIL_ADDRESS>Gaoling School of Artificial Intelligence,
Renmin University of China ,  Binbin Hu<EMAIL_ADDRESS>Ant Group ,
Zhiqiang Zhang<EMAIL_ADDRESS>Ant Group and  Wayne Xin Zhao🖂
<EMAIL_ADDRESS>Gaoling School of Artificial Intelligence, Renmin
University of ChinaBeijing Key Laboratory of Big Data Management and Analysis
MethodsBeijing Academy of Artificial Intelligence

(2022)

###### Abstract.

Session-based Recommendation (SBR) refers to the task of predicting the next
item based on short-term user behaviors within an anonymous session. However,
session embedding learned by a non-linear encoder is usually not in the same
representation space as item embeddings, resulting in the inconsistent
prediction issue while recommending items. To address this issue, we propose a
simple and effective framework named CORE, which can unify the representation
space for both the encoding and decoding processes. Firstly, we design a
representation-consistent encoder that takes the linear combination of input
item embeddings as session embedding, guaranteeing that sessions and items are
in the same representation space. Besides, we propose a robust distance
measuring method to prevent overfitting of embeddings in the consistent
representation space. Extensive experiments conducted on five public real-
world datasets demonstrate the effectiveness and efficiency of the proposed
method. The code is available at: https://github.com/RUCAIBox/CORE.

session-based recommendation

$\dagger$ Work done during internship at Ant Group.

🖂 Corresponding author.

††journalyear: 2022††copyright: acmcopyright††conference: Proceedings of the
45th International ACM SIGIR Conference on Research and Development in
Information Retrieval; July 11–15, 2022; Madrid, Spain.††booktitle:
Proceedings of the 45th International ACM SIGIR Conference on Research and
Development in Information Retrieval (SIGIR ’22), July 11–15, 2022, Madrid,
Spain††price: 15.00††isbn: 978-1-4503-8732-3/22/07††doi:
10.1145/3477495.3531955††ccs: Information systems Recommender systems

Figure 1. (a) Encoder-decoder framework of most existing session-based
recommendation models and (b) Inconsistent prediction issue while measuring
the distance between embeddings for recommending.

##  1\. Introduction

Session-based Recommendation (SBR) aims to capture short-term and dynamic user
preferences according to user behaviors (_e.g.,_ clicked items) within
anonymous sessions, which is an important research topic in recommender
systems (Wang et al., 2019). Recent SBR methods usually follow the encoder-
decoder framework and focus on designing effective neural architecture as
encoder, such as Recurrent Neural Network (RNN) (Hidasi et al., 2016a; Li et
al., 2017; Ren et al., 2019), Transformer (Kang and McAuley, 2018; Sun et al.,
2019; Xie et al., 2020; Zhou et al., 2020; Bian et al., 2021; He et al., 2021;
Zhou et al., 2022) and Graph Neural Network (GNN) (Wu et al., 2019; Xu et al.,
2019; Chen and Wong, 2020; Qiu et al., 2019). As for the decoder, a widely
adopted way is to calculate dot product of session and item embeddings as the
interaction probability of the next item.

In line with the focus of the above-mentioned efforts, we notice that session
embeddings are not usually in the same representation space as item
embeddings. Generally speaking, user behaviors within a short session tend to
share a common focus (Chen and Wong, 2020). In session-based recommendation,
session embeddings are expected to reflect users’ short-term preference and
should be similar to embeddings of preferred items. However, as shown in
Figure 1, when item embeddings are encoded by non-linear neural networks, the
resulting session embedding doesn’t necessarily fall into the space spanned by
bases of item embeddings, and thus it may become unable to concisely represent
user’s preference. For example, as illustrated in Figure 1, for sessions with
a common objective (item $a$ clicked multiple times for simulation), we
observe that their embeddings that encoded by non-linear encoder (_e.g.,_
GRU4Rec (Hidasi et al., 2016a)) fall into different points in embedding space,
giving inconsistent prediction when calculating similarities to item
embeddings.

Considering the above issues, we strive to unify the representation space for
both encoding and decoding processes in session-based recommendation. The
basic idea is to represent sessions in the item embedding space, for example,
by directly summing up item embeddings within the session. However, in this
way, we fail to model the sequential nature, which can be captured by neural
encoders like RNNs. Thus we consider overcoming two key challenges. (1) How to
design a more suitable encoder, so that we can take advantage of deep non-
linear neural networks’ immense capability, given that session and item
embeddings share a consistent representation space. (2) Once the
representation space is unified and item embeddings are directly involved in
score calculation and model optimization, how to measure the distance between
embeddings to avoid overfitting of item embeddings?

To this end, we propose a rather simple yet effective framework for session-
based recommendation, where session embeddings and item embeddings are in
COnsistent REpresentation space, namely CORE. Firstly, we encode session
embeddings as the linear combination of item embeddings within the session,
ensuring the encoded embedding has a consistent representation space as item
embeddings. The weight of each item in the session is learned via a deep non-
linear neural network to incorporate various of inductive biases, such as
order and importance of items. Secondly, we revisit the widely adopt dot
product distance measuring in the perspective of optimizing tuplet loss and
improve the robustness of distance measuring from multiple aspects. Extensive
experiments on five public datasets demonstrate the effectiveness and
efficiency of the proposed framework. We also show that performances of
existing session-based recommendation models can be significantly improved by
the proposed techniques.

##  2\. Methodology

In this section, we introduce the proposed framework CORE, a simple and
effective framework for session-based recommendation within consistent
representation space. Figure 2 presents the overall architecture of CORE.

Now, we start with a brief description of the typical session-based
recommendation, which follows an encoder-decoder framework. Firstly, each item
is embedded into an unified embedding space.
$\bm{h}_{i}=\operatorname{Emb}(v_{i})\in\mathbb{R}^{d}$ denotes the item
embedding for item $v_{i}$, where $\operatorname{Emb}(\cdot)$ is the item
embedding look-up table and $d$ is the dimension of vectors. Then we have
$\bm{h}_{s}=\operatorname{Encoder}([\bm{h}_{s,1},\ldots,\bm{h}_{s,n}])\in\mathbb{R}^{d}$
to encode a session $s$ with $n$ items, where $\operatorname{Encoder}(\cdot)$
is usually a non-linear neural network. Finally, we can predict the
probability distribution for the next item, _i.e.,_
$\bm{\hat{y}}=\operatorname{Decoder}(\bm{h}_{s})\in\mathbb{R}^{m}$, where $m$
is the number of all items.

###  2.1. Representation-Consistent Encoding

Here we aim at encoding sessions into item embedding space to overcome the
inconsistent representation space issue. As most existing encoders adopt a
non-linear encoder (_e.g.,_ RNNs or Transformers) directly stacked over the
input item embeddings, the encoded session embeddings are not in the same
representation space as items. In order to keep embeddings within the same
space, a natural idea is whether we can remove non-linear activation functions
over item embeddings and encode sessions as the linear combination of item
embeddings.

Figure 2. Overall framework of CORE. Table 1. Overall performance comparison
on five datasets. “$*$” indicates the statistical significance for $p<0.01$
compared to the best baseline method with paired $t$-test. Sessions are split
into train/validation/test set in a ratio of 8:1:1 for fair evaluation. We
indicate performances of FPMC on Yoochoose as “$-$” due to the OOM issue.

Dataset | Metric | FPMC | GRU4Rec | NARM | SR-GNN | NISER+ | LESSR | SGNN-HN | SASRec | GC-SAN | CL4Rec | CORE-ave | CORE-trm | Improv.
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Diginetica | R@20 | 31.83 | 45.43 | 47.68 | 48.76 | 51.23 | 48.80 | 50.89 | 49.86 | 50.95 | 50.03 | 50.21 |  52.89* | +3.24%
M@20 | 8.79 | 14.77 | 15.58 | 16.93 | 18.32 | 16.96 | 17.25 | 17.19 | 17.84 | 17.26 | 18.07 |  18.58* | +1.42%
Nowplaying | R@20 | 10.18 | 13.80 | 14.17 | 15.28 | 16.55 | 17.60 | 16.75 | 20.69 | 18.30 | 20.59 | 20.31 |  21.81* | +5.41%
M@20 | 4.51 | 5.83 | 6.11 | 6.10 | 7.14 | 7.13 | 6.13 | 8.14 | 8.13 | 8.21 | 6.62 | 7.35 | $-$
RetailRocket | R@20 | 46.04 | 55.32 | 58.65 | 58.71 | 60.36 | 56.22 | 58.82 | 59.81 | 60.18 | 59.69 | 59.18 |  61.85* | +2.47%
M@20 | 21.95 | 33.18 | 34.69 | 36.42 | 37.43 | 37.11 | 35.72 | 36.03 | 36.85 | 35.95 |  37.52* |  38.76* | +3.55%
Tmall | R@20 | 20.30 | 23.25 | 31.67 | 33.65 | 35.97 | 32.45 | 39.14 | 35.82 | 35.32 | 35.59 |  44.67* |  44.48* | +14.13%
M@20 | 13.07 | 15.78 | 21.83 | 25.27 | 27.06 | 23.96 | 23.46 | 25.10 | 23.48 | 25.07 |  31.85* |  31.72* | +17.70%
Yoochoose | R@20 | $-$ | 60.78 | 61.67 | 61.84 | 62.99 | 62.89 | 62.49 | 63.55 | 63.24 | 63.61 | 58.83 |  64.61* | +1.57%
M@20 | $-$ | 27.27 | 27.82 | 28.15 | 28.98 | 28.59 | 28.24 | 28.63 | 29.00 | 28.73 | 25.05 | 28.24 | $-$

Along this line, we propose a Representation-Consistent Encoder (RCE), whose
output session embedding is the weighted summarization of item embeddings
within a session. Linear combination guarantees the session embeddings are
always in the same embedding space as items. Although non-linear layers are
removed between the encoded session embedding and input item embeddings, they
are essential to incorporating inductive biases and learning weights for input
item embeddings. Formally, we apply arbitrary deep neural networks to learn
weights for each item embedding in a session,

(1) $\displaystyle\bm{\alpha}=$
$\displaystyle\operatorname{DNN}([\bm{h}_{s,1};\bm{h}_{s,2};\ldots;\bm{h}_{s,n}]),$
(2) $\displaystyle\bm{h}_{s}=$
$\displaystyle\sum_{i=1}^{n}\alpha_{i}\bm{h}_{s,i},$

Then we give two detailed implementations of $\operatorname{DNN}$ without
carefully designed architecture.

####  2.1.1. Learning Weights via Mean Pooling.

This variant adopts a mean pooling layer as $\operatorname{DNNs}$, _i.e.,_
$\alpha_{i}=\frac{1}{n}$. This variant ignores the order of items in a
session, as well as the importance of each item.

####  2.1.2. Learning Weights via Transformer.

In this variant, we utilize $L$-layers self-attention blocks like SASRec (Kang
and McAuley, 2018) as $\operatorname{DNNs}$.

(3)
$\displaystyle\bm{F}=\operatorname{Transformers}([\bm{h}_{s,1};\bm{h}_{s,2};\ldots;\bm{h}_{s,n}]),$

where $\bm{F}\in\mathbb{R}^{n\times d^{\prime}}$ and $d^{\prime}$ is the
output dimension of feed forward network of the last layer of self-attention
blocks. Then we can obtain the normalized weights
$\bm{\alpha}\in\mathbb{R}^{n}$,

(4)
$\displaystyle\bm{\alpha}=\operatorname{softmax}(\bm{w}\cdot\bm{F}^{\top}),$

where $\bm{w}\in\mathbb{R}^{d^{\prime}}$ are learnable parameters. This
variant captures the sequential nature via positional encoding technique in
Transformer.

###  2.2. Robust Distance Measuring for Decoding

As sessions are encoded as linear combination of item embeddings and decoded
by measuring the distances to items in embedding space, item embeddings are
directly involved in the distance calculation between embeddings, leading to a
high risk of overfitting. Thus we seek a robust way to measure the distance in
the unified representation space to prevent overfitting. By reviewing the
widely adopted dot product distance, we have the following lemma:

###### Lemma 1.

Given a session embedding $\bm{h}_{s}$ and item embeddings
$\\{\bm{h}_{v}|v\in\mathcal{V}\\}$, when dot product is used to measure the
embedding distance, optimizing cross entropy loss is approximately
proportional to optimize $(N-1)$-tuplet loss (Sohn, 2016) with a fixed margin
of $2$.

###### Proof.

The lemma can be proved by rewriting the original loss function as,

$\displaystyle\ell_{\text{ori}}=$
$\displaystyle-\log\frac{\exp(\bm{h}_{s}\cdot\bm{h}_{v^{+}})}{\sum_{i=1}^{m}\exp(\bm{h}_{s}\cdot\bm{h}_{v_{i}})}$
$\displaystyle=$
$\displaystyle\log\left[1+(|\mathcal{V}|-1)\sum\limits_{v^{-}\in\mathcal{V}\backslash\\{v^{+}\\}}\exp(\bm{h}_{s}\bm{h}_{v^{-}}-\bm{h}_{s}\bm{h}_{v^{+}})\right],$
$\displaystyle\simeq$
$\displaystyle(|\mathcal{V}|-1)\sum\limits_{v^{-}\in\mathcal{V}\backslash\\{v^{+}\\}}\exp(\bm{h}_{s}\bm{h}_{v^{-}}-\bm{h}_{s}\bm{h}_{v^{+}})$
$\displaystyle\simeq$
$\displaystyle(|\mathcal{V}|-1)\sum\limits_{v^{-}\in\mathcal{V}\backslash\\{v^{+}\\}}(\bm{h}_{s}\bm{h}_{v^{-}}-\bm{h}_{s}\bm{h}_{v^{+}}+1),$
$\displaystyle\propto$
$\displaystyle\sum_{v^{-}\in\mathcal{V}\backslash\\{v^{+}\\}}\left(\|\bm{h}_{s}-\bm{h}_{v^{+}}\|^{2}-\|\bm{h}_{s}-\bm{h}_{v^{-}}\|^{2}+2\right).$

where $m$ denotes the number of items, $v^{+}$ denotes the ground-truth next
item for session $s$. ∎

Given the above lemma, we consider improving the robustness of distance
measuring in three aspects. Firstly, data distributions for different
recommendation scenarios may vary a lot, and a fixed margin is inappropriate.
Thus, we replace the fixed margin $2$ by a controllable hyper-parameter $\tau$
to suit different scenarios. Secondly, we utilize Dropout (Srivastava et al.,
2014), a widely adopted technique for robust training, directly over candidate
item embeddings. Thirdly, inspired by recent advances in Contrastive Learning
(Chen et al., 2020; He et al., 2020; Lin et al., 2022), we propose to measure
the distance via cosine distance for better alignment and uniformity of item
embeddings (Wang and Isola, 2020). Then we design the loss function with the
proposed Robust Distance Measuring (RDM) technique as,

(5)
$\displaystyle\ell=-\log\frac{\exp\left(\cos(\bm{h}_{s},\bm{h}^{\prime}_{v^{+}})/\tau\right)}{\sum_{i=1}^{m}\exp\left(\cos(\bm{h}_{s},\bm{h}^{\prime}_{v_{i}})/\tau\right)},$

where $\bm{h}^{\prime}$ denotes the item embeddings with dropout.

##  3\. Experiments

Table 2. Statistics of the datasets.

Dataset | # Interactions | # Items | # Sessions | Avg. Length
---|---|---|---|---
Diginetica | 786,582 | 42,862 | 204,532 | 4.12
Nowplaying | 1,085,410 | 59,593 | 145,612 | 9.21
RetailRocket | 871,637 | 51,428 | 321,032 | 6.40
Tmall | 427,797 | 37,367 | 66,909 | 10.62
Yoochoose | 1,434,349 | 19,690 | 470,477 | 4.64

Datasets and evaluation metrics. We conduct experiments on five public
datasets collected from real-world platforms: _Diginetica_ , _Nowplaying_ ,
_RetailRocket_ , _Tmall_ and _Yoochoose_ with their statistics shown in Table
2. We filter out sessions of length 1 and items appearing less than 5 times
across all datasets (Li et al., 2017; Wu et al., 2019; Chen and Wong, 2020),
and split the sessions in each dataset into train/validation/test set in
temporal order in a ratio of 8:1:1 (Zhao et al., 2020). To evaluate the
performance of different methods, we employ two widely-used metrics, top-$20$
Recall (R@$20$) and top-$20$ Mean Reciprocal Rank (M@$20$) (Wu et al., 2019;
Yu et al., 2020).

Baselines. To evaluate the performance of the proposed method, we compare it
with the following representative baselines: (1) Matrix factorization based
methods: FPMC (Rendle et al., 2010); (2) RNN-based methods: GRU4Rec (Hidasi et
al., 2016a) and NARM (Li et al., 2017); (3) GNN-based methods: SR-GNN (Wu et
al., 2019), NISER+ (Gupta et al., 2019), LESSR (Chen and Wong, 2020) and SGNN-
HN (Pan et al., 2020). (4) Transformer-based methods: SASRec (Kang and
McAuley, 2018), GC-SAN (Xu et al., 2019) and CL4Rec (Xie et al., 2020). Note
that we don’t take methods that introduce additional collaborative filtering
information (Wang et al., 2020; Jiang et al., 2020; Xia et al., 2021; Cho et
al., 2021) or other side features (Hidasi et al., 2016b) as baselines. For
CORE, we implement two simple variants using mean pooling (Sec. 2.1.1) and
Transformers (Sec. 2.1.2) as $\operatorname{DNNs}$ in Eqn. (1), namely _CORE-
ave_ and _CORE-trm_ , respectively.

Implementation details The proposed models and all the baselines are
implemented based on a popular open-source recommendation library
RecBole111https://recbole.io (Zhao et al., 2021) and its extension RecBole-
GNN222https://github.com/RUCAIBox/RecBole-GNN for easy development and
reproduction. The dimension of the latent vectors is fixed to $100$, and each
session is truncated within a maximum length of $50$. We optimized all the
compared methods using Adam optimizer (Kingma and Ba, 2015) with a learning
rate of $0.001$, and adopted early-stop training if the M@$20$ performance on
the validation set decreased for $5$ consecutive epochs. We use a batch size
of $2048$ for all methods (except TAGNN, for which we use $100$ due to the
large memory consumption). Other hyper-parameters of baselines are carefully
tuned following the suggestions from the original papers and we report each
performance under its optimal settings. For CORE, we apply a grid search for
controllable margin $\tau$ among $\\{0.01,0.05,0.07,0.1,1\\}$ and dropout
ratio $\rho$ among $\\{0,0.1,0.2\\}$. We finally report metrics on the test
set with models that gain the highest performance on the validation set.

Figure 3. Performances over training time relative to that of CORE-ave on
Diginetica.

Overall comparison. From the experimental results in Table 1, we can observe:
CORE outperforms all the baselines significantly over $8$ of $10$ metrics on
the adopted five datasets. Different from these baselines, CORE doesn’t apply
non-linear layers over item embeddings to encode sessions, but learns weights
for each item embedding and adopts a weighted sum for encoding session
embeddings and item embeddings in consistent representation space. Besides,
notice that CORE-ave considers neither item order in a session, nor item
importances. However, it still outperforms all the baselines on $3$ metrics,
while gaining comparable results with several strong baselines on other
metrics. Without carefully designed encoder architecture, CORE can achieve
impressive performance, which further confirms the importance of encoding
session embeddings and item embeddings in consistent representation space.

Analysis 1: Efficiency. In Figure 3, we plot the performance of several
popular session-based recommendation models as well as CORE-trm over their
training time per epoch relative to that of CORE-ave on the Diginetica
dataset. We measure the training time on an NVIDIA TITAN V GPU. As we can see,
CORE-ave is the fastest while achieving competitive performance. By applying
mean pooling over sessions, CORE-ave minimizes memory usage and only learns a
single item embedding table during training. CORE-trm has a similar training
time as SASRec while achieving the best performance among the compared
baselines.

Table 3. Ablation study of CORE’s variants on Diginetica and RetailRocket. Method | Diginetica | RetailRocket
---|---|---
R@20 | M@20 | R@20 | M@20
CORE | $52.89$ | $18.58$ | $61.85$ | $38.76$
w/o RCE | 49.82 | 17.41 | 59.59 | 36.27
w/o RDM | 52.31 | 18.38 | 60.93 | 37.72
SASRec | 49.86 | 17.19 | 59.81 | 36.03
Table 4. Performance comparison of different methods and their improved variants on two datasets. Method | Diginetica | RetailRocket
---|---|---
R@20 | M@20 | R@20 | M@20
NARM | 47.68 | 15.58 | 58.65 | 34.69
\+ RCE | 51.86 | 18.27 | 60.77 | 37.01
\+ RDM | 51.62 | 17.79 | 61.33 | 37.11
\+ All | $52.51$ | $18.58$ | $62.19$ | $38.84$
SR-GNN | 48.76 | 16.93 | 58.71 | 36.42
\+ RCE | 49.51 | 17.53 | 57.05 | 35.70
\+ RDM | 51.36 | 18.57 | 61.41 | 38.27
\+ All | $52.38$ | $18.95$ | $61.43$ | $38.38$

Analysis 2: Ablation study. CORE involves several components (_i.e.,_ RCE and
RDM) and we now analyze how each part contributes to the performance. As
SASRec and CORE-trm (CORE for simplify) shares the same Transformer
architecture, we select SASRec as the base model to compare. We mainly
consider the following variants: CORE w/o RCE means replacing representation-
consistent encoder of CORE to SASRec’s encoder; CORE w/o RDM means replacing
the proposed robust distance measuring techniques to the traditional dot
product distance; In Table 3, we can observe that the performance order can be
summarized as CORE $>$ CORE w/o RDM $>$ CORE w/o RCE $\simeq$ SASRec. These
results indicate that all the parts are useful to improve the final
performance.

Analysis 3: Improving existing methods with RCE & RDM. Here we slightly modify
several popular existing session-based models, showing how the proposed
components RCE and RDM improve the performances of existing methods.

Take _SR-GNN(Wu et al., 2019) \+ RCE_ as an example, the original encoder can
be formulated as $\bm{h}_{s}=\bm{W}[\bm{h}_{s,n};\bm{h}_{g}]$, where
$\bm{W}\in\mathbb{R}^{d\times 2d}$ are learnable parameters and $\bm{h}_{g}$
can be seen as the linear combination of item embeddings within a session.
Then we remove $\bm{W}$ and change the encoder to
$\bm{h}_{s}=(\bm{h}_{s,n}+\bm{h}_{g})/2$ and the modified SR-GNN encoder is a
variant of the proposed RCE, where the session embedding is linear combination
of item embeddings.

Performance comparison between existing works (NARM (Li et al., 2017) and SR-
GNN (Wu et al., 2019)) and their improved variants are shown in Table 4. We
can see that generally the combination of RCE and RDM can gain dramatic
performance improvements compared to the original models. As for only adding
one of the proposed techniques, RDM consistently improve the performance,
while RCE has a positive effect in most cases.

(a) GRU4Rec

(b) SASRec

(c) CORE (ours)

Figure 4. Visualization of learned session embeddings.

Analysis 4: Visualization of session embeddings. To show the effectiveness of
the proposed model CORE, we visualize the learned session embeddings using
t-SNE (Van der Maaten and Hinton, 2008) algorithm in Figure 4. Detailed,
sessions with the same next item are viewed as in the same class, and are
marked the same color. We randomly sample $15$ items as ground truth next
items. Then we extract all the corresponding sessions from our test set in
Diginetica. We can see that different classes’ session embeddings learned by
CORE are well separated from each other compared to those learned by GRU4Rec
and SASRec.

Analysis 5: Parameter tuning. At last, we exime the impact of several
importance parameters towards CORE, _i.e.,_ the margin $\tau$ and the dropout
ratio $\rho$. In particular, we vary $\tau$ in $\\{0.01,0.02,\ldots,0.1\\}$
and present the results in Figure 5 (a) and Figure 5 (b). Our method
constantly outperforms the best baseline, and achieves the best performance
when $\tau=0.07$ on Diginetica and $\tau=0.08$ on RetailRocket. Overall, the
performance is stable around $0.04\leq\tau\leq 0.08$. Next, we vary $\rho$ in
the range $0$ and $0.5$ with steps of $0.1$. As shown in Figure 5 (c) and
Figure 5 (d), when $\rho\leq 0.3$, the item dropout increases the robustness
of item embedding learning, while a large ratio of item dropout may hurt the
performance, as we can see, the performance decreases sharply when $\rho>0.3$.

(a) Diginetica

(b) RetailRocket

(c) Diginetica

(d) RetailRocket

Figure 5. Parameter tuning of CORE on Diginetica and RetailRocket datasets.

##  4\. Conclusion

In this paper, we propose CORE, a simple and effective framework for session-
based recommendation within consistent representation space, which unifies the
representation space throughout encoding and decoding to overcome the
inconsistent prediction issue. Different from stacking multiple non-linear
layers over the item embeddings, we propose to just apply weighted sum for
item embeddings to encode sessions in the consistent representation space as
items. Besdies, we propose robust distance measuring techniques from multiple
aspects to prevent overfitting of item embeddings in the proposed framework.
Extensive experiments on five public datasets have shown the effectiveness and
efficiency of the proposed approach, as well as how the proposed techniques
can help existing methods.

For future work, we will consider studying the expressive ability of the
proposed representation-consistent encoder both theoretically and empirically.
Besides, we will explore how to introduce side features and useful inductive
biases to the proposed framework.

###### Acknowledgements.

This work was partially supported by the National Natural Science Foundation
of China under Grant No. 61872369 and 61832017, Beijing Outstanding Young
Scientist Program under Grant No. BJJWZYJH012019100020098, and CCF-Ant Group
Research Fund. This work was partially supported by Beijing Academy of
Artificial Intelligence (BAAI). Xin Zhao is the corresponding author.

## References

  * (1)
  * Bian et al. (2021) Shuqing Bian, Wayne Xin Zhao, Kun Zhou, Jing Cai, Yancheng He, Cunxiang Yin, and Ji-Rong Wen. 2021.  Contrastive Curriculum Learning for Sequential User Behavior Modeling via Data Augmentation. In _CIKM_.
  * Chen et al. (2020) Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey E. Hinton. 2020.  A Simple Framework for Contrastive Learning of Visual Representations. In _ICML_. 1597–1607.
  * Chen and Wong (2020) Tianwen Chen and Raymond Chi-Wing Wong. 2020.  Handling Information Loss of Graph Neural Networks for Session-based Recommendation. In _KDD_.
  * Cho et al. (2021) Junsu Cho, SeongKu Kang, Dongmin Hyun, and Hwanjo Yu. 2021\.  Unsupervised Proxy Selection for Session-based Recommender Systems. In _SIGIR_. 327–336.
  * Gupta et al. (2019) Priyanka Gupta, Diksha Garg, Pankaj Malhotra, Lovekesh Vig, and Gautam M Shroff. 2019.  NISER: Normalized Item and Session Representations with Graph Neural Networks.  _arXiv preprint arXiv:1909.04276_ (2019).
  * He et al. (2020) Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross B. Girshick. 2020.  Momentum Contrast for Unsupervised Visual Representation Learning. In _CVPR_. 9726–9735.
  * He et al. (2021) Zhankui He, Handong Zhao, Zhe Lin, Zhaowen Wang, Ajinkya Kale, and Julian McAuley. 2021\.  Locally constrained self-attentive sequential recommendation. In _CIKM_.
  * Hidasi et al. (2016a) Balázs Hidasi, Alexandros Karatzoglou, Linas Baltrunas, and Domonkos Tikk. 2016a.  Session-based Recommendations with Recurrent Neural Networks. In _ICLR_.
  * Hidasi et al. (2016b) Balázs Hidasi, Massimo Quadrana, Alexandros Karatzoglou, and Domonkos Tikk. 2016b.  Parallel Recurrent Neural Network Architectures for Feature-rich Session-based Recommendations. In _RecSys_. 241–248.
  * Jiang et al. (2020) Xunqiang Jiang, Binbin Hu, Yuan Fang, and Chuan Shi. 2020\.  Multiplex memory network for collaborative filtering. In _SDM_.
  * Kang and McAuley (2018) Wang-Cheng Kang and Julian J. McAuley. 2018.  Self-Attentive Sequential Recommendation. In _ICDM_. 197–206.
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba. 2015.  Adam: A Method for Stochastic Optimization. In _ICLR_.
  * Li et al. (2017) Jing Li, Pengjie Ren, Zhumin Chen, Zhaochun Ren, Tao Lian, and Jun Ma. 2017\.  Neural Attentive Session-based Recommendation. In _CIKM_. 1419–1428.
  * Lin et al. (2022) Zihan Lin, Changxin Tian, Yupeng Hou, and Wayne Xin Zhao. 2022\.  Improving Graph Collaborative Filtering with Neighborhood-enriched Contrastive Learning. In _TheWebConf_.
  * Pan et al. (2020) Zhiqiang Pan, Fei Cai, Wanyu Chen, Honghui Chen, and Maarten de Rijke. 2020.  Star Graph Neural Networks for Session-based Recommendation. In _CIKM_.
  * Qiu et al. (2019) Ruihong Qiu, Jingjing Li, Zi Huang, and Hongzhi Yin. 2019\.  Rethinking the Item Order in Session-based Recommendation with Graph Neural Networks. In _CIKM_.
  * Ren et al. (2019) Pengjie Ren, Zhumin Chen, Jing Li, Zhaochun Ren, Jun Ma, and Maarten de Rijke. 2019\.  RepeatNet: A Repeat Aware Neural Recommendation Machine for Session-Based Recommendation. In _AAAI_. 4806–4813.
  * Rendle et al. (2010) Steffen Rendle, Christoph Freudenthaler, and Lars Schmidt-Thieme. 2010\.  Factorizing personalized Markov chains for next-basket recommendation. In _WWW_. 811–820.
  * Sohn (2016) Kihyuk Sohn. 2016\.  Improved Deep Metric Learning with Multi-class N-pair Loss Objective. In _NIPS_. 1849–1857.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey E. Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. 2014.  Dropout: a simple way to prevent neural networks from overfitting.  _JMLR_ 15, 1 (2014), 1929–1958.
  * Sun et al. (2019) Fei Sun, Jun Liu, Jian Wu, Changhua Pei, Xiao Lin, Wenwu Ou, and Peng Jiang. 2019.  BERT4Rec: Sequential Recommendation with Bidirectional Encoder Representations from Transformer. In _CIKM_. 1441–1450.
  * Van der Maaten and Hinton (2008) Laurens Van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-SNE.  _JMLR_ 9, 11 (2008).
  * Wang et al. (2019) Shoujin Wang, Longbing Cao, and Yan Wang. 2019.  A Survey on Session-based Recommender Systems.  _arXiv preprint arXiv:1902.04864_ (2019).
  * Wang and Isola (2020) Tongzhou Wang and Phillip Isola. 2020.  Understanding Contrastive Representation Learning through Alignment and Uniformity on the Hypersphere. In _ICML_ , Vol. 119. 9929–9939.
  * Wang et al. (2020) Ziyang Wang, Wei Wei, Gao Cong, Xiao-Li Li, Xianling Mao, and Minghui Qiu. 2020\.  Global Context Enhanced Graph Neural Networks for Session-based Recommendation. In _SIGIR_.
  * Wu et al. (2019) Shu Wu, Yuyuan Tang, Yanqiao Zhu, Liang Wang, Xing Xie, and Tieniu Tan. 2019\.  Session-Based Recommendation with Graph Neural Networks. In _AAAI_. 346–353.
  * Xia et al. (2021) Xin Xia, Hongzhi Yin, Junliang Yu, Qinyong Wang, Lizhen Cui, and Xiangliang Zhang. 2021\.  Self-Supervised Hypergraph Convolutional Networks for Session-based Recommendation. In _AAAI_.
  * Xie et al. (2020) Xu Xie, Fei Sun, Zhaoyang Liu, Shiwen Wu, Jinyang Gao, Bolin Ding, and Bin Cui. 2020.  Contrastive Learning for Sequential Recommendation.  _arXiv preprint arXiv:2010.14395_ (2020).
  * Xu et al. (2019) Chengfeng Xu, Pengpeng Zhao, Yanchi Liu, Victor S. Sheng, Jiajie Xu, Fuzhen Zhuang, Junhua Fang, and Xiaofang Zhou. 2019\.  Graph Contextualized Self-Attention Network for Session-based Recommendation. In _IJCAI_. 3940–3946.
  * Yu et al. (2020) Feng Yu, Yanqiao Zhu, Qiang Liu, Shu Wu, Liang Wang, and Tieniu Tan. 2020.  TAGNN: Target Attentive Graph Neural Networks for Session-based Recommendation. In _SIGIR_.
  * Zhao et al. (2020) Wayne Xin Zhao, Junhua Chen, Pengfei Wang, Qi Gu, and Ji-Rong Wen. 2020.  Revisiting Alternative Experimental Settings for Evaluating Top-N Item Recommendation Algorithms. In _CIKM_.
  * Zhao et al. (2021) Wayne Xin Zhao, Shanlei Mu, Yupeng Hou, Zihan Lin, Kaiyuan Li, Yushuo Chen, Yujie Lu, Hui Wang, Changxin Tian, Xingyu Pan, Yingqian Min, Zhichao Feng, Xinyan Fan, Xu Chen, Pengfei Wang, Wendi Ji, Yaliang Li, Xiaoling Wang, and Ji-Rong Wen. 2021.  RecBole: Towards a Unified, Comprehensive and Efficient Framework for Recommendation Algorithms. In _CIKM_.
  * Zhou et al. (2020) Kun Zhou, Hui Wang, Wayne Xin Zhao, Yutao Zhu, Sirui Wang, Fuzheng Zhang, Zhongyuan Wang, and Ji-Rong Wen. 2020\.  S3-Rec: Self-Supervised Learning for Sequential Recommendation with Mutual Information Maximization. In _CIKM_.
  * Zhou et al. (2022) Kun Zhou, Hui Yu, Wayne Xin Zhao, and Ji-Rong Wen. 2022\.  Filter-enhanced MLP is All You Need for Sequential Recommendation. In _TheWebConf_.


tocchapterList of Publications

[LE] [RO]
[RE] [LO]


# Belief propagation generalizes backpropagation

Frederik Eaton
<EMAIL_ADDRESS>
(September 29, 2022)

###### Abstract

The two most important algorithms in artificial intelligence are
backpropagation and belief propagation. In spite of their importance, the
connection between them is poorly characterized. We show that when an input to
backpropagation is converted into an input to belief propagation so that
(loopy) belief propagation can be run on it, then the result of belief
propagation encodes the result of backpropagation; thus backpropagation is
recovered as a special case of belief propagation. In other words, we prove
for apparently the first time that belief propagation generalizes
backpropagation. Our analysis is a theoretical contribution, which we motivate
with the expectation that it might reconcile our understandings of each of
these algorithms, and serve as a guide to engineering researchers seeking to
improve the behavior of systems that use one or the other.

##  1 Introduction

In this paper we consider a connection between two algorithms which could be
said to have the status of being the two most fundamental algorithms in the
various fields of computer science concerned with the numerical modeling of
real world systems, these fields being sometimes known as artificial
intelligence or machine learning, sometimes called control theory or
statistical modeling or approximate inference. The two algorithms that form
our subject matter are usually referred to as backpropagation and belief
propagation, respectively, although these are modern terms for concepts that
go back several hundred years in Western thought [1, 2, 3].

Back-propagation is another name for the chain rule of differential calculus
[4], applied iteratively to a network of functions, or in other words to a
function of functions of multiple independent variables and other such
functions; the input to backpropagation may also be known in the field as a
“deep network” or a “neural network” [5].

Belief propagation, by contrast, takes as its input a network of probability
distributions, also called a probabilistic network. Belief propagation is
equivalent to an iterative application of Bayes’ rule, which is the rule for
inferring the posterior distribution $P(X|Y)$ of a random variable $X$ from
its prior $P(X)$ and some table of conditional probabilities $P(Y|X)$
describing the possible observations of some related variable $Y$: $P(X|Y)$ =
${1\over Z(Y)}P(X)P(Y|X)$, often stated as $P(X|Y)={P(X)P(Y|X)\over P(Y)}$ [6,
7]. Belief propagation forbids the sharing of variables among multiple tables
of conditionals in its input, implying that there must be no loops in the
network, but the symmetry between the posteriors $P(X|Y)$ and the possible
observations $P(Y|X)$ allows this dynamic programming algorithm to be recast
as a sequence of “message updates” which apply equally well to networks with
variable sharing and therefore loops, as was observed by Pearl in the 1980s
[8]. The generalized, message-based form of belief propagation is called
“loopy belief propagation” [9] (it is also called the “sum-product algorithm”
[10]), and while this generalization can no longer be said to compute exact
posteriors for its variables, its numerical behavior and its usefulness in
approximation have been the object of much study, and it enjoys widespread
application in certain specialized domains such as error correcting codes [11,
12]. Loopy belief propagation, to restate the above definition, is “exact on
trees”, a tree being a network with no loops, in which case loopy belief
propagation reduces to belief propagation. Its theoretical properties are
otherwise difficult to characterize [13], and much research has been directed
towards improving the accuracy of belief propagation by accounting for the
presence of loops in the input model [14, 15, 16, 17, 18].

Figure 1: Example network

This paper considers another class of inputs for which loopy belief
propagation computes exact quantities, namely probabilistic networks that
arise through a straightforward “lifting”111 Our terminology. For a similar
use of the term ”lifting” in probabilistic inference see [19]; this is
connected to type-theoretic lifting [20]; but not to be confused with lifted
inference [21, 22], type lifting in compilers [23], or von Neumann’s concept
of lifting in measure theory [24]. of function networks. It is simple to show
that the “delta function” posteriors computed by loopy belief propagation on
these networks are exact, as they are just a probabilistic representation of
the deterministic computation embodied in the original function network.
Moreover, we show for the first time that when the output node of the lifted
network is attached to a Boltzmann distribution [25, 26] prior, the messages
that propagate backwards through the network encode a representation of the
exact derivatives of the output variable with respect to each other variable,
making loopy belief propagation on the lifted network an extended or lifted
form of backpropagation on the original function network. This second result
is the main contribution of the paper, establishing that belief propagation is
a generalization of backpropagation.

##  2 Example model

We find it useful to introduce the concepts of this paper through a small
example function network. We define a network containing a single loop, and a
single shared variable $x$, from the following equations:

$\displaystyle z=f(u,v)\quad u=g(w,x)\quad v=h(x,y)\quad x=j(t)$ (1)

The network is illustrated in figure 1.

##  3 Backpropagation

Running backpropagation on this function network means calculating the
derivative of $z$ with respect to the six other variables recursively using
the chain rule of differential calculus, starting with the variable $u$

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}u}={\partial z\over\partial
u}\equiv{\partial f\over\partial u}\equiv f^{(u)}(u,v)$ (2)

where we have used “$\equiv$” to show the equivalence of alternate notations
for partial derivatives. Then for $v$ we have

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}v}=f^{(v)}(u,v)$ (3)

and for $w$

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}w}=\frac{\mathrm{d}z}{\mathrm{d}u}{\partial
u\over\partial w}\equiv\frac{\mathrm{d}z}{\mathrm{d}u}g^{(w)}(w,x)$ (4)

and so on. The term “adjoint” is used as a shorthand: since $\mathrm{d}z$
always appears in the numerator in backpropagation, rather than write “the
derivative of $z$ with respect to $w$” we call this quantity “the adjoint of
$w$ [with respect to $z$]” [27]. Calculating the adjoint of $x$ requires our
first addition:

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}x}=\frac{\mathrm{d}z}{\mathrm{d}u}{\partial
u\over\partial x}+\frac{\mathrm{d}z}{\mathrm{d}v}{\partial v\over\partial
x}\equiv\frac{\mathrm{d}z}{\mathrm{d}u}g^{(x)}(w,x)+\frac{\mathrm{d}z}{\mathrm{d}v}h^{(x)}(x,y)$
(5)

The last adjoint to be calculated is $t$:

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}t}=\frac{\mathrm{d}z}{\mathrm{d}x}j^{(t)}(x)$
(6)

For a general function network, the chain rule would be written [4, 28]

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}x_{i}}=\sum_{k\succ
i}\frac{\mathrm{d}z}{\mathrm{d}x_{k}}{\partial f_{k}\over\partial x_{i}}$ (7)

where “$k\succ i$” means iterating over the parents $k$ of $i$, and where the
general network is defined as a collection of functions and variables

$\displaystyle x_{k}=f_{k}(\\{x_{i}|i\prec k\\})$ (8)

The adjoints of parent variables are calculated and recorded before their
children in a backward pass over the network, giving rise to the term
“backpropagation” [5].

##  4 Lifting

We are interested in knowing what happens when we try to run belief
propagation on our network, but first we have to convert the function network
into a probabilistic network with continuous real-valued variables. To use
belief propagation in this setting, we must represent the variables in our
network as probability density functions. This requires that we first define a
probability distribution over the real numbers which places all of its mass on
a single value:

$\displaystyle P(x=0)=1,\quad P(x\neq 0)=0$ (9)

The density for this distribution is called the “Dirac delta function” [29],
written $\delta(x)$. This is not a true function since it is infinite at
$x=0$, but we can think of it as a limit of functions, for example a limit of
Gaussians whose standard deviation tends to zero (see figure 2):

$\displaystyle\delta(x)=\lim_{\sigma\to 0}{1\over\sigma\sqrt{2\pi}}e^{-{1\over
2}({x\over\sigma})^{2}}$ (10)

Although this limit itself is not well-defined, it tells us symbolically how
to treat the delta function when it appears inside an integral, namely by
doing the integral first and then taking the limit:

$\displaystyle\int f(x)\delta(x)\mathrm{d}x\equiv\lim_{\sigma\to 0}\int
f(x){1\over\sigma\sqrt{2\pi}}e^{-{1\over
2}({x\over\sigma})^{2}}\mathrm{d}x=f(0)$ (11)

It may be that an algorithmic implementation of our proposed lifting would
approximate delta functions with very narrow Gaussians, in which case we still
expect belief propagation to be well-behaved, but we do not go into an
analysis of that behavior here.

Figure 2: Gaussian distributions getting narrower

The lifting operation simply replaces each function node $z=f(u,v)$ with a
positive-valued “factor” defined on all three variables:

$\displaystyle F(u,v,z)=\delta(f(u,v)-z)$ (12)

which encodes the functional relationship as a density. When working with such
expressions, one must remember that the distinction between input and output
variables hasn’t been entirely lost; it is not the case that
$\delta(y-f(x))=\delta(f^{-1}(y)-x)$, because there is a Jacobian scaling
factor:

$\displaystyle\delta(f(x))$
$\displaystyle=\left|f^{(x)}(x)\right|^{-1}\delta(x),\textrm{\quad hence}$
(13) $\displaystyle\delta(y-f(x))$
$\displaystyle=\left|f^{(x)}(x)\right|^{-1}\delta(f^{-1}(y)-x)$ (14)

##  5 Belief propagation

The message updates for (loopy) belief propagation can be written concisely by
defining two types of messages, messages going from variables to factors, and
messages going from factors to variables [10]. Messages only go between
variables and the factors to which they are immediately connected; both types
of messages are represented as positive functions of the variable involved.
The message from a variable to a factor is simply the product of all the
messages coming from the other factors to that variable. For example,
referring to figure 3, which shows the lifted form of the example network, the
message from $x$ to $J$ is updated as:

$\displaystyle m_{(x,J)}(x):=m_{(G,x)}(x)m_{(H,x)}(x)$ (15)

The message from a factor to a variable is calculated by multiplying the
factor by all of the messages coming into the factor from other variables, and
then integrating (or summing) over the other variables. For example, the
message from $H$ to $x$ is updated as follows:

$\displaystyle m_{(H,x)}(x):=\int
H(x,y,v)m_{(v,H)}(v)m_{(y,H)}(y)\mathrm{d}v\mathrm{d}y$ (16)

Convergence of the message updates is usually independent of their initial
values, but for simplicity we assume that they are initialized to a constant:

$\displaystyle m^{0}_{(x,J)}(x):=1\quad m^{0}_{(H,x)}(x):=1\quad\textrm{etc.}$
(17)

With ordinary (non-loopy) belief propagation, for efficiency the different
messages are updated in a single forward and single backward pass over the
network; any further updates would leave them unchanged, so they can be said
to have converged at this point. Readers who have encountered Hidden Markov
Models [30, 31], and their continuous, real-valued counterpart the Kalman
filter [32], will be familiar with these forward and backward passes, which
are examples of belief propagation on these specialized probabilistic
networks.

With loopy belief propagation, the messages may be updated in any order. The
order of message updates may affect the rate of convergence, but not the final
values to which the messages converge, as long as convergence is achieved.

After the messages have converged, the posterior of each variable is estimated
as the product of the messages coming into it:

$\displaystyle
P(x)\approx{1\over\int\mathrm{d}x}m_{(G,x)}(x)m_{(H,x)}(x)m_{(J,x)}(x)$ (18)

where ${1\over\int\mathrm{d}x}$ represents a normalization constant.

###  5.1 General form of belief propagation messages

For reference, we now give the message updates of belief propagation for a
general probabilistic network, consisting of a set of factors
$\\{F_{\alpha}\\}$ and variables $\\{x_{i}\\}$ (see [10]). The message from a
factor $F_{\alpha}$ to a variable $x_{i}$ is updated as:

$\displaystyle m_{(F_{\alpha},x_{i})}(x_{i}):=\int
F_{\alpha}(x_{\alpha})\prod_{j\sim\alpha\setminus
i}m_{(x_{j},F_{\alpha})}(x_{j})\>\mathrm{d}x_{\alpha\setminus i}$ (19)

where the subscripts $i$ and $j$ index variables in the network, the subscript
$\alpha$ which indexes factors also represents a set of variables neighboring
the respective factor, and $j\sim\alpha\setminus i$ denotes any variable $j$
neighboring $\alpha$ except $i$.

The update for a message from a variable to a factor is similarly written:

$\displaystyle m_{(x_{i},F_{\alpha})}(x_{i}):=\prod_{\beta\sim
i\setminus\alpha}m_{(F_{\beta},x_{i})}(x_{i})$ (20)

##  6 Running the lifted model

In order to simulate evaluating the original function network on a given set
of inputs, we assign priors to the input variables of the probabilistic
network, which we do by attaching single-variable factors to them. These
factors are just delta functions encoding the input values to the original
function network. For example, if the input values are

$\displaystyle w=w^{*},\quad t=t^{*},\quad y=y^{*}$ (21)

then we introduce factors

$\displaystyle W^{*}(w)=\delta(w-w^{*}),\quad T^{*}(t)=\delta(t-t^{*}),$
$\displaystyle Y^{*}(y)=\delta(y-y^{*})$ (22)

These factors will cause messages to propagate upwards through the network
which consist of delta functions that encode the computation of the original
function network.

Finally, we introduce a factor $B$ assigning a Boltzmann prior to the output
node. We omit the arbitrary temperature constant, which can be recovered by
replacing $e$ with $\exp({1/kT})$ in our notation.

$\displaystyle B(z)=e^{z}$ (23)

The use of this prior could be seen as representing our desire to maximize the
output of the function network. We will show that it causes messages to be
propagated downward through the network that make it possible for derivatives
to be calculated locally at each node. The Boltzmann prior is not a true
probability distribution, since it is not normalizable, but this is not a
concern for messages.

The lifted example network is shown in figure 3. The upward delta function
messages have been omitted, but example downward messages are illustrated with
plots next to each edge.

Figure 3: Lifted example network, with messages

##  7 Behavior of the lifted model

We are now interested in understanding the behavior of belief propagation on
the lifted model. This behavior is specified to a great extent by the
topological structure of the probabilistic network, which inherits certain
properties from the fact that it derives from a function network. Each factor
was originally a function with one or more input variables and only one
output, and with each variable occurring as the output of at most one
function. Therefore although there is a loop in the network, the network’s
structure is not entirely general, and interactions between messages are
channeled in such a way that certain invariants are maintained by the message
updates. In addition to the distinction defined earlier between variable-
factor and factor-variable messages, it is possible to assign an “upwards” or
“downwards” direction to each message. We find that when running loopy belief
propagation on a network that was produced by our lifting transformation,
messages propagating upwards through the network are delta functions, but the
downward messages can take an arbitrary form and are not converted into delta
functions by any of the upward messages. When the algorithm converges, the
posteriors at each variable node are delta functions centered at the value of
the variable in the original function network computation, but the downward
messages are nevertheless able to encode additional information from which the
values of derivatives may be obtained.

To see why the upward messages are allowed to take a separate form from the
downward messages, note that each variable has only one downward message
leaving it, towards the factor that it represents the output of, and that this
message is calculated only as the product of other downward messages coming
from factors for which it had served as an input. Thus, although the product
of any function with a delta function is another delta function, no delta
functions enter into the product when downward messages are updated according
to the variable$\to$factor message updates (equation 20).

The first downward message in the example network comes from the Boltzmann
prior $B$:

$\displaystyle m_{(B,z)}(z)=e^{z}$ (24)

This is propagated unchanged from $z$ to $F$:

$\displaystyle m_{(z,F)}(z)=e^{z}$ (25)

We next calculate the message from $F$ to $U$:

$\displaystyle m_{(F,u)}(u)=\int
F(u,v,z)m_{(z,F)}(z)m_{(v,F)}(v)\mathrm{d}z\mathrm{d}v$ (26)

Now $m_{(v,F)}(v)$ is an upward message and therefore a delta function,
$\delta(v-v^{*})$. Substituting $F$ according to our lifting, we have

$\displaystyle m_{(F,u)}(u)$ (27)
$\displaystyle=\int\delta(z-f(u,v))m_{(z,F)}(z)\delta(v-v^{*})\mathrm{d}v\mathrm{d}z$
(28) $\displaystyle=m_{(z,F)}(f(u,v^{*}))$ (29)
$\displaystyle=\exp(f(u,v^{*}))$ (30)

This is propagated without change to $G$, the only other neighbor of $u$:

$\displaystyle m_{(u,G)}(u)=\exp(f(u,v^{*}))$ (31)

Similarly we have

$\displaystyle m_{(v,H)}(v)=m_{(F,v)}(v)=\exp(f(u^{*},v))$ (32)

Similarly to the message from $F$ to $u$, the message from $G$ to $x$ must
incorporate the upward message $m_{(w,G)}(w)$, which is a delta function
$\delta(w-w^{*})$:

$\displaystyle m_{(G,x)}(x)$ (33)
$\displaystyle=\int\delta(u-g(w,x))m_{(u,G)}(u)\delta(w-w^{*})\mathrm{d}u\mathrm{d}w$
(34) $\displaystyle=m_{(u,G)}(g(w^{*},x))$ (35)
$\displaystyle=\exp\Bigl{(}f(g(w^{*},x),v^{*})\Bigr{)}$ (36)

And we see that, correspondingly for $H$,

$\displaystyle m_{(H,x)}(x)=\exp\Bigr{(}f(u^{*},h(x,y^{*}))\Bigl{)}$ (37)

The message from $x$ to $J$ is simply the product of these two messages:

$\displaystyle
m_{(x,J)}(x)=\exp\Bigl{(}f(g(w^{*},x),v^{*})+f(u^{*},h(x,y^{*}))\Bigr{)}$ (38)

Notice that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}\log m_{(x,J)}(x)={\partial
f\over\partial u}\frac{\mathrm{d}u}{\mathrm{d}x}+{\partial f\over\partial
v}\frac{\mathrm{d}v}{\mathrm{d}x}=\frac{\mathrm{d}f}{\mathrm{d}x}$ (39)

when evaluated at $x=x^{*}$ (and with $w=w^{*}$ and so on) which gives
$\frac{\mathrm{d}f}{\mathrm{d}x}$ according to the chain rule. The
relationship only holds when the derivative is evaluated at $x=x^{*}$, because
$u^{*}$ and $v^{*}$ depend on $x^{*}$, and they appear as constants in the two
terms.

It remains to show that the relationship of equation 39 holds more generally.
Setting aside the example network, let us assume we are given an arbitrary
function network and its lifted counterpart, a probabilistic network on which
we have executed belief propagation. We want to establish two invariants which
hold for the messages in the network. These invariants apply to downward
messages of both types and relate them to the variable adjoints, which is to
say the derivatives of an objective variable, $z$, with respect to each
variable.

######  Theorem 1.

The following invariant holds for the downward message from any variable $x$
to the factor $F$ adjacently below it:

$\displaystyle\left.\frac{\mathrm{d}}{\mathrm{d}x}\right|_{x=x^{*}}\log
m_{(x,F)}(x)=\frac{\mathrm{d}z}{\mathrm{d}x}$ (a)

And the following invariant holds for the downward message from a factor $F$
with output $y$ to one of its neighboring (input) variables, $x$:

$\displaystyle\left.\frac{\mathrm{d}}{\mathrm{d}x}\right|_{x=x^{*}}\log
m_{(F,x)}(x)=\frac{\mathrm{d}z}{\mathrm{d}y}\left({\partial y\over\partial
x}\equiv{\partial f\over\partial x}\right)$ (b)

We prove invariants (a) and (b) using induction, by assuming that they already
hold for all the downward-directed messages above the current edge in the
network, and expanding the current message using the message update rules of
equations 19 and 20. After substituting equation 20 into invariant (a), we get

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}\log\left(m_{(x,F)}(x)=\prod_{G\succ
x}m_{(G,x)}(x)\right)$ (40)

where $G\succ x$ represents any factor $G$ above the variable $x$ in the
network. Since $x$ must be the output node of $F$, this product iterates over
all the neighbors of $x$ not equal to $F$, as specified by the message update
rule. This becomes

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}\log\prod_{G\succ x}m_{(G,x)}(x)$
$\displaystyle=\sum_{G\succ x}\left(\frac{\mathrm{d}}{\mathrm{d}x}\log
m_{(G,x)}(x)\right)$ (41) $\displaystyle=\sum_{G\succ
x}\frac{\mathrm{d}z}{\mathrm{d}g}{\partial g\over\partial
x}=\frac{\mathrm{d}z}{\mathrm{d}x}$ (42)

the second equality following from the induction hypothesis and invariant (b).
The lower-case $g$ stands for the function encoded by the factor $G$ and its
output node. This summation is $\frac{\mathrm{d}z}{\mathrm{d}x}$ by the chain
rule, which establishes invariant (a) for the message from $x$ to its child
$F$.

To prove the second invariant, we substitute equation 19 into (b), which is to
say

$\displaystyle m_{(F,x)}(x)=\int$
$\displaystyle\delta\left(\hat{f}-f(x,\\{y\\})\right)$ (43) $\displaystyle
m_{(\hat{f},F)}(\hat{f})\prod_{y}m_{(y,F)}(y)\mathrm{d}\hat{f}\mathrm{d}\\{y\\}$
(44)

where $\hat{f}$ signifies the output variable associated with the function
$f$, and $y$ represents all the inputs of $f$ except $x$. As with equation 28
above (in our analysis of the example model), the messages $m_{(y,F)}(y)$ are
all delta functions, so the substitution becomes

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}\log$ $\displaystyle\
m_{(F,x)}(x)=\frac{\mathrm{d}}{\mathrm{d}x}\log
m_{(\hat{f},F)}(f(x,\\{y^{*}\\}))$ (45)
$\displaystyle=\left(\frac{\mathrm{d}}{\mathrm{d}\hat{f}}\log
m_{(\hat{f},F)}(\hat{f})\right){\partial f\over\partial
x}=\frac{\mathrm{d}z}{\mathrm{d}f}{\partial f\over\partial x}$ (46)

where the last equality follows from the induction hypothesis and invariant
(a).

We must finally prove the “base case” of the induction, namely that invariant
(a) holds for the message from $z$ to the function node $F$ directly below it.
Since the only other neighbor of $z$ is the Boltzmann factor $B$, this message
is equal to the message from $B$ to $z$:

$\displaystyle m_{(z,F)}(z)=m_{(B,z)}(z)=e^{z}$ (47)

Invariant (a) then becomes

$\displaystyle\left.\frac{\mathrm{d}}{\mathrm{d}z}\right|_{z=z^{*}}\log\left(m_{(z,F)}(z)=e^{z}\right)=\frac{\mathrm{d}}{\mathrm{d}z}\log
e^{z}=1=\frac{\mathrm{d}z}{\mathrm{d}z}$ (48)

and so it is satisfied for the base case. This completes the proof by
induction.

We have described running belief propagation on a network where the
independent variables are assigned delta function priors:

$\displaystyle X^{*}(x)=\delta(x-x^{*})$ (49)

In the case where these delta functions are approximated using a more general
distribution such as a narrow Gaussian, it may be more useful to estimate the
adjoints of our variables using a form of invariant (a) that does not depend
on choosing a specific value $x^{*}$ at which to evaluate the derivative.

######  Corollary 1.

$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}x}=-\int\left({\partial\over\partial
x}X^{*}(x)\right)\log m_{(x,X^{*})}(x)\mathrm{d}x$ (50)

For the case of delta functions, this is equivalent to invariant (a) by
analogy to the following integration by parts identity:

$\displaystyle\int\left({\partial\over\partial
x}\delta(x)\right)f(x)\mathrm{d}x=-{\partial f\over\partial x}(0)$ (51)

But when $X^{*}$ is a Gaussian, for example, the above expression 50 for
$\frac{\mathrm{d}z}{\mathrm{d}x}$ is equivalent to a kind of smoothed
numerical differentiation.

Our proof by induction makes it clear that as with backpropagation, the
converged belief propagation messages in our model can be calculated in a
single forward and single backward pass over the network.

##  8 Conclusions

###  8.1 Motivation

Most papers in machine learning seek to introduce a new computer algorithm to
the field. The purpose of this paper is rather to shed light on a connection
between two well-established algorithms, to provide groundwork for a better
theoretical understanding of both algorithms, and to eliminate some of the
mystery surrounding them for students.

Belief propagation and backpropagation apply to different classes of input
model. Belief propagation applies to probabilistic models and is used in
domains where there is a need to model uncertainty directly, and
backpropagation applies to deterministic models, where it is used to provide
gradients to support the fitting of such models to data. Because all real-
world data contains some measure of uncertainty, there is considerable overlap
between these two domains, and it could be said that any essential difference
between them is only a matter of engineering philosophy; based on the
engineer’s decision about whether to model uncertainty directly or indirectly,
and at which level of the system to do so [33, 34].

There has been recent interest in extensions of backpropagation that
incorporate uncertainty more directly into the algorithm; some of these, such
as stochastic gradient descent [35] or drop-out [36], apply backpropagation to
inputs which change at random; others, such as probabilistic backpropagation
[37], extend backpropagation by replacing deterministic quantities with
probabilistic representations of the same quantities, somewhat related to the
“lifting” we refer to in this paper. We hope that it would be possible to
assist these investigations by clarifying the mathematical relationship
between backpropagation and belief propagation.

The problem of characterizing the behavior of belief propagation on a lifted
function network whose inputs have been initialized with distributions other
than delta functions remains an open question. In this case, we can expect in
general that the converged messages will not produce exact posteriors and will
not lead to exact adjoints being calculated, because the downward messages
will have the effect of slightly changing the variable locations specified in
the upward messages (figure 4) and these effects will be compounded as the
messages propagate around loops. We do not know whether a loop-corrected form
of belief propagation would be necessary to make this more general scenario
useful.

Figure 4: A Gaussian message being shifted to the right after multiplication
by an exponential

However, before learning the inputs to a function network using gradient
descent, the precise value of the input variables is in general unknown. Being
able to make this uncertainty explicit at a more basic level, by running some
form of backpropagation on a “lifted” probabilistic version of the network,
with inputs that are not delta functions, could be desirable for a number of
reasons, for example because it allows the convergence rate of the input
variables to be reflected in their posterior distributions, or because it
allows some of the input variables to be specified with less certainty than
others, which could provide an evolving indicator of where the training
algorithm should focus its attention.

Readers who are interested in probabilistic approaches to the problem of
training “neural networks” could start with David MacKay’s thesis [38] which
proposes approximations that could be used to model uncertainty at the level
of variables in the network. Extensions to this idea are explored in for
example [39] and [40]; more recently, [41] points out that by replacing
backpropagation with message passing, it becomes easier to train networks that
have discrete weights, which can be useful for hardware-based network
implementations with limited numerical precision. “Probabilistic
backpropagation” [37] is the name given to an approach that combines a forward
pass that approximates the distribution at each network node as a Gaussian,
with a backwards pass that backpropagates adjoints of these distribution
parameters. Experiments show that the method compares favorably with plain
backpropagation and with Hamiltonian Monte-Carlo, a probabilistic training
method based on sampling [42], although there appear to be many details in the
implementation. Our paper is less concerned with experimental results, and
more concerned with making a sea of different ideas more navigable by pointing
out some overlooked connections that exist within it.

Belief propagation and backpropagation are both useful for analyzing large
models because they have the same time complexity as running the model itself.
Like strict (non-loopy) belief propagation, backpropagation is a dynamic
programming algorithm that requires only two passes over the input network,
the first pass serving to compute the value of the objective or output
variable, and the second pass serving to compute the derivatives. Loopy belief
propagation, on the other hand, exchanges numerical messages locally on the
network for some usually small number of iterations, typically until the
messages converge; for error-correcting codes and certain other applications,
convergence is fast enough that the algorithm does not add significant time
complexity [12, 11, 43, 44]. While belief propagation and backpropagation both
distinguish between input and output variables, loopy belief propagation
requires no such distinction to be made.

Framing an algorithm in terms of locally-exchanged messages can be useful for
distributing it across multiple computers, and there may be some value derived
from being able to rethink backpropagation in terms of iterative local
message-passing. Another contribution of this paper is to show that by placing
backpropagation in the framework of loopy belief propagation, the input-output
relationships of backpropagation become part of the messages rather than being
hard-coded through the functions of the network, and the original function
network can be inverted with respect to one of the input variables, simply by
moving the Boltzmann prior onto this variable while leaving the rest of the
network unchanged.

###  8.2 Generality

It is desirable to point out that the form of the “lifting” of a function
network to a probabilistic network which we describe here is a straightforward
requirement of the problem of converting from one class of inputs to the
other. The “Dirac delta function” is a well-understood formalism for
specifying a probability distribution that takes only a single value, and it
is used to lift both variables and functions into the domain of probabilities.

The use of the Boltzmann distribution is motivated as follows: backpropagation
is most commonly used to solve optimization problems; the most natural way of
converting an optimization problem to a probabilistic inference problem is to
place a Boltzmann distribution over the objective: $p(E)\propto\exp{E\over
kT}$, where $E$ is the objective or output variable of the function network,
and $kT$ is a constant specifying the tightness of the distribution around the
optimum. This distribution has its origins in thermodynamics, where it
describes the distribution over the states of a system with energy $E$ and
temperature $T$ [45]. Also called the Gibbs measure in mathematical contexts,
the Boltzmann distribution has widespread use in machine learning, for example
in stochastic neural networks, see for example the “Boltzmann machine” [46,
47]; and arises almost universally in probability theory in a less
recognizable form, the exponential family model, which appears whenever data
consist of exchangeable observations [48].

There are a few hurdles to overcome in attempting to unify belief propagation
and backpropagation. The first is that the domains of each algorithm are
different, one being probabilistic and the other being deterministic. This is
addressed by our “lifting” transformation, but the transformation produces
models that are considered less tractable than a typical input to belief
propagation: first of all, a typical lifted function network will contain many
loops; and secondly, a function network operates on real-valued rather than
discrete variables. Computing the message updates of belief propagation on the
lifted network requires some difficult modeling decisions: how to represent
distributions over real variables, whether to represent delta functions
specially or as a limit of narrow Gaussians, how to perform the numerical
integration required by the message updates, and how to represent the
Boltzmann distribution and other messages which may be unnormalizable. All of
these hurdles can be surmounted in various ways. There is a relatively long
history of the successful use of belief propagation and related message-
passing algorithms to perform efficient probabilistic inference in real-valued
probability networks, see for example “assumed density filtering” and
“expectation propagation” [49].

Finally, this work has relevance to researchers seeking to invent novel ways
to improve the training phase of models based on function networks, which are
seeing increasingly widespread application in computer science. Rather than
changing the structure or mathematical relationships of the network to make it
behave more tractably under a backpropagation-based training method, one could
instead consider tuning its independent variables by applying belief
propagation (or one of its many extensions) to a lifted, probabilistic,
version of the network where it is possible to reason about uncertainty more
directly.

To this end, it would seem helpful to observe that the original
backpropagation algorithm is recovered exactly by loopy belief propagation in
the case where the network is initialized with delta functions, just as it has
often been helpful in the analysis of loopy belief propagation on general
probabilistic networks to observe that the algorithm is exact in the case
where the input network is a tree.

## References

  * [1] Isaac Newton.  Philosophiæ Naturalis Principia Mathematica.  Royal Society Press, 1687.
  * [2] G.W. Leibniz.  The Early Mathematical Manuscripts of Leibniz.  The Open Court Publishing Company, 1920.
  * [3] H.A. Bethe.  Statistical Theory of Superlattices.  Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 150(871):552–575, 1935.
  * [4] J.L. Lagrange.  Théorie des fonctions analytiques.  Courcier, 1797.
  * [5] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  Nature, 323(6088):533–536, 1986.
  * [6] T. Bayes.  An Essay towards Solving a Problem in the Doctrine of Chances. By the Late Rev. Mr. Bayes, F. R. S. Communicated by Mr. Price, in a Letter to John Canton, A. M. F. R. S.  Philosophical Transactions of the Royal Society of London, 53:370–418, 1763.
  * [7] P.S. Laplace.  Théorie analytique des probabilités.  Courcier, 1812.
  * [8] Judea Pearl.  Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference.  Morgan Kaufmann, 1988.
  * [9] Brendan J. Frey, Ralf Koetter, and Nemanja Petrovic.  Very loopy belief propagation for unwrapping phase images.  Advances in Neural Information Processing Systems, 14, 2001.
  * [10] F.R. Kschischang, B.J. Frey, and H.A. Loeliger.  Factor graphs and the sum-product algorithm.  IEEE Transactions on information theory, 47(2):498–519, 2001.
  * [11] N. Wiberg, H.A. Loeliger, and R. Kotter.  Codes and iterative decoding on general graphs.  European Transactions on telecommunications, 6(5):513–525, 1995\.
  * [12] David J.C. MacKay and Radford M. Neal.  Good codes based on very sparse matrices.  In IMA International Conference on Cryptography and Coding, pages 100–111. Springer, 1995.
  * [13] Alexander T. Ihler, John W. Fisher III, Alan S. Willsky, and David Maxwell Chickering.  Loopy belief propagation: convergence and effects of message errors.  Journal of Machine Learning Research, 6(5), 2005.
  * [14] Kyomin Jung and Devavrat Shah.  Inference in binary pair-wise markov random fields through self-avoiding walks.  Preprint on http://arxiv. org/abs/cs. AI/0610111v2, 2006.
  * [15] M. Chertkov and V.Y. Chernyak.  Loop series for discrete statistical models on graphs.  Journal of Statistical Mechanics: Theory and Experiment, 2006, 2006\.
  * [16] Yusuke Watanabe and Kenji Fukumizu.  Loop series expansion with propagation diagrams.  Journal of Physics A: Mathematical and Theoretical, 42(4):045001, 2008.
  * [17] M.J. Wainwright, T.S. Jaakkola, and A.S. Willsky.  Tree-reweighted belief propagation algorithms and approximate ML estimation by pseudomoment matching.  In Workshop on Artificial Intelligence and Statistics, volume 21, 2003.
  * [18] J.M. Mooij, B. Wemmenhove, H.J. Kappen, and T. Rizzo.  Loop corrected belief propagation.  In Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, 2007.
  * [19] Oleg Kiselyov and Chung-chieh Shan.  Embedded probabilistic programming.  In IFIP Working Conference on Domain-Specific Languages, pages 360–384. Springer, 2009.
  * [20] Ralf Hinze.  Lifting operators and laws, 2010.
  * [21] David Poole.  First-order probabilistic inference.  In IJCAI, volume 3, pages 985–991, 2003.
  * [22] Rodrigo de Salvo Braz, Eyal Amir, and Dan Roth.  Lifted first-order probabilistic inference.  Statistical relational learning, page 433, 2007.
  * [23] Bratin Saha and Zhong Shao.  Optimal type lifting.  In International workshop on types in compilation, pages 156–177. Springer, 1998.
  * [24] Alexandra Ionescu Tulcea and Cassius Ionescu Tulcea.  On the lifting property.  Technical report, Yale Univ Dept of Mathematics, 1961.
  * [25] Ludwig Boltzmann.  Über die mechanische Bedeutung des zweiten Hauptsatzes der Wärmetheorie.  Staatsdruckerei, 1866.
  * [26] Ludwig Boltzmann.  Weitere studien über das wärmegleichgewicht unter gasmolekülen (1872).  In Kinetische Theorie II, pages 115–225. Springer, 1970.
  * [27] Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  Journal of Marchine Learning Research, 18:1–43, 2018.
  * [28] Jerrold E. Marsden, Anthony J. Tromba, and Alan Weinstein.  Basic multivariable calculus.  Springer, 1993.
  * [29] Paul Adrien Maurice Dirac.  The principles of quantum mechanics.  Oxford university press, 1981.
  * [30] Leonard E. Baum and Ted Petrie.  Statistical inference for probabilistic functions of finite state markov chains.  The annals of mathematical statistics, 37(6):1554–1563, 1966.
  * [31] X.D. Huang, Y. Ariki, and M.A. Jack.  Hidden Markov Models for Speech Recognition.  Edinburgh University Press, 1991.
  * [32] R. E. Kalman.  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82(1):35–45, 03 1960.
  * [33] Edwin T. Jaynes.  Probability theory: The logic of science.  Cambridge university press, 2003.
  * [34] Kevin P. Murphy.  Machine learning: a probabilistic perspective.  MIT press, 2012.
  * [35] Jack Kiefer and Jacob Wolfowitz.  Stochastic estimation of the maximum of a regression function.  The Annals of Mathematical Statistics, pages 462–466, 1952.
  * [36] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research, 15(56):1929–1958, 2014.
  * [37] Jose Miguel Hernandez-Lobato and Ryan Adams.  Probabilistic backpropagation for scalable learning of bayesian neural networks.  In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37, pages 1861–1869, 07–09 Jul 2015.
  * [38] David John Cameron Mackay.  Bayesian methods for adaptive models.  PhD thesis, California Institute of Technology, 1992.
  * [39] Sara A Solla and Ole Winther.  Optimal perceptron learning: an online bayesian approach.  On-Line Learning in Neural Networks. Cambridge University Press, Cambridge, 1998.
  * [40] Alfredo Braunstein and Riccardo Zecchina.  Learning by message passing in networks of discrete synapses.  Physical review letters, 96(3):030201, 2006.
  * [41] Daniel Soudry, Itay Hubara, and Ron Meir.  Expectation backpropagation: Parameter-free training of multilayer neural networks with continuous or discrete weights.  Advances in neural information processing systems, 27, 2014.
  * [42] Radford M Neal.  Bayesian Learning for Neural Networks.  PhD thesis, University of Toronto, 1995.
  * [43] Tom Richardson.  The geometry of turbo-decoding dynamics.  IEEE Transactions on Information Theory, 46(1):9–23, 2000.
  * [44] Kevin Murphy, Yair Weiss, and Michael I. Jordan.  Loopy belief propagation for approximate inference: An empirical study.  arXiv preprint arXiv:1301.6725, 2013.
  * [45] Enrico Fermi.  Thermodynamics.  Blackie, 1936.
  * [46] David Sherrington and Scott Kirkpatrick.  Solvable model of a spin-glass.  Physical Review Letters, 35:1792–1796, Dec 1975.
  * [47] David Ackley, Geoffrey Hinton, and Terrence Sejnowsky.  A learning algorithm for boltzmann machines.  Cognitive Science, 1985.
  * [48] M.J. Wainwright and M.I. Jordan.  Graphical models, exponential families, and variational methods.  New Directions in Statistical Signal Processing, 2003.
  * [49] Tom Minka.  A family of algorithms for approximate Bayesian inference.  PhD thesis, MIT, 2001.