bag model nodes, bag mappings $f_{B}$ were implemented by dense neural
networks as well, however, the form of aggregation $g$ is yet to be specified.
As shown in the last chapter, simple $\operatorname{mean}$ aggregation is
sufficient in theory. Nevertheless, in our experience, models using more types
of aggregations than only $\operatorname{mean}$ perform better in practice.
Aggregation in HMill bag model node takes form of a concatenation
$g=\mathopen{}\mathclose{{}\left[g_{1}|g_{2}|\ldots g_{q}}\right]$ of $q$ base
aggregation functions
$g_{i}\colon\cup_{k=1}^{\infty}{\mathopen{}\mathclose{{}\left(\mathbb{R}^{m}}\right)}^{k}\to\mathbb{R}^{m}$.
Different choices of base aggregation $g_{i}$, or their combinations, are
suitable for different problems. Nevertheless, because the input is
interpreted as an unordered bag of instances, every aggregation function
should be invariant to permutation of its parameters and also should not scale
with increasing size of the bag. In this section, several options for $g_{i}$
are proposed. All base aggregations proposed in this thesis are element-wise
operations—each element in the resulting vector is computed using only entries
at the same position in input. Given a set of vectors
$\\{\bm{x}_{1},\ldots,\bm{x}_{k}\\}$, $g_{i}$ would first select from all
vectors their element at the first index and compute one number from the
resulting set of real numbers
$\mathopen{}\mathclose{{}\left\\{x_{1},\ldots,x_{k}}\right\\}$, then all
elements at the second index, and so forth. This way, we obtain $m$ numbers,
or in other words, a vector from $\mathbb{R}^{m}$, a space of the same
dimension as is the size of vectors on input. For simplicity, all the
following formulas below represent computation of only one output dimension of
each base aggregation. We repeat the procedure $q$ times for all base
aggregations, which leads to the final result in $\mathbb{R}^{qm}$.

####  4.2.1 Non-parametric aggregations

The most simple choices for base aggregation $g_{i}$ are $\max$:
$g_{\max}(\\{x_{1},\ldots,x_{k}\\})=\max_{i=1,\ldots,k}x_{i}$ (4.3) and
$\operatorname{mean}$:
$g_{\operatorname{mean}}(\\{x_{1},\ldots,x_{k}\\})=\frac{1}{k}\sum_{i=1}^{k}x_{i}$
(4.4) The motivation for either of this form of base aggregation function
comes from the standard assumption and the collective assumption introduced in
Sections 2.1.1 and 2.1.2, respectively. However, let us emphasize that in
HMill instances may be represented not only by vectors in real vector spaces,
but by entire sample (sub)trees. The $\max$ aggregation is suitable for cases
when one instance in the bag may give evidence strong enough to predict the
label. This is oftentimes the case in the computer security domain, where
signals indicating the maliciousness are rather rare, yet usually strongly
incriminating. On the other side of the spectrum lies the
$\operatorname{mean}$ aggregation function, which detects well trends
identifiable globally over the whole bag.

####  4.2.2 Parametric aggregations

Whereas non-parametric aggregations do not use any parameter, parametric
aggregations represent an entire class of functions parametrized by a real
vector of parameters, which is learned during training. In this thesis, we use
two parametric aggregations—$\operatorname{lse}$ (log-sum-exp) [Kraus et al.,
2015] and $\operatorname{\mathnormal{p}-norm}$ [Gulcehre et al., 2014]. Log-
sum-exp aggregation is parametrized by a vector of positive numbers
$\bm{r}\in{(\mathbb{R}^{+})}^{m}$ that specifies one real parameter for
computation in each output dimension and is computed as follows:
$g_{\operatorname{lse}}(\\{x_{1},\ldots,x_{k}\\};r)=\frac{1}{r}\log\mathopen{}\mathclose{{}\left(\frac{1}{k}\sum_{i=1}^{k}\exp({r\cdot
x_{i}})}\right)$ (4.5) Here, $r$ denotes an element of $\bm{r}$, corresponding
to the dimension, where the result is being computed. With different values of
$r$, log-sum-exp behaves differently and in fact both $\max$ (4.3) and
$\operatorname{mean}$ (4.4) are limiting cases for $\operatorname{lse}$:
$\displaystyle\lim_{r\to
0}g_{\operatorname{lse}}(\mathopen{}\mathclose{{}\left\\{x_{i}}\right\\}_{i=1}^{k};r)$
$\displaystyle=g_{\operatorname{mean}}(\mathopen{}\mathclose{{}\left\\{x_{i}}\right\\}_{i=1}^{k})$
(4.6)
$\displaystyle\lim_{r\to\infty}g_{\operatorname{lse}}(\mathopen{}\mathclose{{}\left\\{x_{i}}\right\\}_{i=1}^{k};r)$
$\displaystyle=g_{\max}(\mathopen{}\mathclose{{}\left\\{x_{i}}\right\\}_{i=1}^{k})$
(4.7) If $r$ is very small, the output of the function approaches simple
average (4.6), on the other hand, if $r$ is a large number, log-sum-exp
becomes a smooth approximation of the $\max$ function (4.7). Since computing
the value of $g_{\operatorname{lse}}$ according to (4.5) may lead to numerical
instabilities, the following equality can be derived:
$g_{\operatorname{lse}}(\\{x_{i}\\}_{i=1}^{k};r)=\alpha+g_{\operatorname{lse}}(\\{x_{i}-\alpha\\}_{i=1}^{k};r)$
(4.8) which holds for any $r>0$ and $\alpha\in\mathbb{R}$.333Proofs of all
three claims can be found in Appendix A (Statements 3, 5 and 6). Therefore, to
avoid overflow in summation of exponents, we use $\alpha=\max_{i}x_{i}$.
Finally, to satisfy the constraint that parameters $\bm{r}$ are strictly
positive, we represent them as a vector of (unconstrained) real values
$\bm{\rho}\in\mathbb{R}^{m}$. To evaluate $g_{\operatorname{lse}}$, we re-
parametrize with $\rho\to\log(1+\exp(\rho))$ (softplus
function444Specifically, its stable implementation
$\rho\to\max(\rho,0)+\operatorname{log1p}1p(\exp(-|\rho|))$.) or any other
transformation that maps from $\mathbb{R}$ to $(0,\infty)$.

(a) $r\to 0$

(b) $r=1$

(c) $r=2$

(d) $r\to\infty$

Figure 4.1: Contours of $g_{\operatorname{lse}}(x_{1},x_{2};r)$ base
aggregation defined in (4.5) for the case when $k=2$ and $r$ takes different
values. In case (4.1(a)), where $r\to 0$, $\operatorname{lse}$ aggregation is
equivalent to the $\operatorname{mean}$ aggregation and in case (4.1(d)),
$\operatorname{lse}$ is equivalent to $\max$. Here, $x$ axis represents values
of $x_{1}$ and $y$ axis values of $x_{2}$, however, as all aggregations are
invariant to permutations, contours would be the same even with swapped axes.
Figure 4.2: A visual representation of
$g_{\operatorname{\mathnormal{p}-norm}}(x_{1},x_{2};p,c)$ when $k=2$, and
colors represent different values of $p$. Here
$\color[rgb]{0.05078125,0.30859375,0.546875}\definecolor[named]{pgfstrokecolor}{rgb}{0.05078125,0.30859375,0.546875}p=1$,
$\color[rgb]{0.625,0.109375,0.05078125}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.109375,0.05078125}p=2$,
$\color[rgb]{0.02734375,0.63671875,0.125}\definecolor[named]{pgfstrokecolor}{rgb}{0.02734375,0.63671875,0.125}p=5$
and
$\color[rgb]{0.484375,0.00390625,0.63671875}\definecolor[named]{pgfstrokecolor}{rgb}{0.484375,0.00390625,0.63671875}p\to\infty$
. For clarity, $c$ is set to $0$ and also normalization $1/k$ is not used. The
plot represents a set of points $(x_{1},x_{2})$ for which
$g_{\operatorname{\mathnormal{p}-norm}}(x_{1},x_{2};p,0)=1$.  Last aggregation
is (normalized) $\operatorname{\mathnormal{p}-norm}$, which is parametrized by
a vector of real numbers $\bm{p}$, where $\forall
i\in\mathopen{}\mathclose{{}\left\\{1,\ldots,m}\right\\}\colon p_{i}\geq 1$,
and another vector $\bm{c}\in\mathbb{R}$. It is computed with formula:
$g_{\operatorname{\mathnormal{p}-norm}}(\\{x_{1},\ldots,x_{k}\\};p,c)=\mathopen{}\mathclose{{}\left(\frac{1}{k}\sum_{i=1}^{k}|x_{i}-c|^{p}}\right)^{\frac{1}{p}}$
(4.9) Due to the parametric interpretation of both $p$ and $c$, the
$\operatorname{\mathnormal{p}-norm}$ aggregation (4.9) also generalizes on
several well-known cases (for simplicity we assume here that $c=0$):
$g_{\operatorname{\mathnormal{p}-norm}}(\mathopen{}\mathclose{{}\left\\{x_{i}}\right\\}_{i=1}^{k};p,0)=\begin{cases}\frac{1}{k}\sum_{i=1}^{k}\lvert
x_{i}\rvert,&\text{for }p=1\\\
\sqrt{\frac{1}{k}\sum_{i=1}^{k}(x_{i})^{2}},&\text{for }p=2\\\
\max_{i=1,\ldots,k}\lvert x_{i}\rvert,&\text{for }p\to\infty\\\ \end{cases}$
(4.10) Hence, for different values of $p$,
$\operatorname{\mathnormal{p}-norm}$ (4.9) interpolates between the
(normalized) $\ell_{1}$ norm and the $\ell_{\infty}$ norm.555For a proof for
the case where $p\to\infty$, see Statement 7 in Appendix A. Parameter $c$
allows to shift the origin for each coordinate. Figure 4.2 represents the
behavior of the $\operatorname{\mathnormal{p}-norm}$ aggregation. Similarly to
the case of $\operatorname{lse}$, for stable implementation we make use of
equality:666A proof can be found in Appendix A, Statement 4.
$g_{\operatorname{\mathnormal{p}-norm}}(\\{x_{i}\\}_{i=1}^{k};p,c)=\beta\cdot
g_{\operatorname{\mathnormal{p}-norm}}(\\{x_{i}/\beta\\}_{i=1}^{k};p,c/\beta)$
(4.11) which is valid for any $\beta>0,p\geq 1,c\in\mathbb{R}$. Therefore, we
set $\beta=\max\\{1,\max_{i}\lvert x_{i}-c\rvert\\}$. Finally, to satisfy the
constraint that all elements of vector $p$ are greater or equal to 1, we
employ a similar trick as in the case of the $\operatorname{lse}$ aggregation,
this time using transformation $\rho\to 1+\log(1+\exp(\rho))$.

####  4.2.3 Weighted aggregations

Following the motivation in 4.1, we define the weighted versions of
$\operatorname{mean}$ and $\operatorname{\mathnormal{p}-norm}$ aggregations:
$\displaystyle g_{\operatorname{mean}}(\\{(x_{i},w_{i})\\}_{i=1}^{k})$
$\displaystyle=\frac{1}{\sum_{i=1}^{k}w_{i}}\sum_{i=1}^{k}w_{i}\cdot x_{i}$
(4.12) $\displaystyle
g_{\operatorname{\mathnormal{p}-norm}}(\\{x_{i},w_{i}\\}_{i=1}^{k};p,c)$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{1}{\sum_{i=1}^{k}w_{i}}\sum_{i=1}^{k}w_{i}\cdot|x_{i}-c|^{p}}\right)^{\frac{1}{p}}$
(4.13) Note that since one weight is defined for each instance in the bag,
weights influence result identically across output dimensions of $g$, whereas
aggregation parameters $\bm{r}$, $\bm{p}$ and $\bm{c}$ define one parameter
for each dimension and therefore may differ. The former case corresponds to
altering a contribution of the given instance to the aggregation output. The
latter case on the other hand allows to aggregate representations obtained
from instance models differently according to dimension. For instance, some
elements of the instance representation are best aggregated with the
$\operatorname{mean}$ function and others with the $\max$ function. By
considering $\bm{r}$ a learnable vector of parameters, a suitable $\bm{r}$ can
be inferred for each output dimension during training. To finish, since $\max$
and its smooth approximation $\operatorname{lse}$ identify extreme phenomena
across the bag, it is reasonable not to use weights in their computation.
Thus, weights are ignored: $\displaystyle
g_{\max}(\\{(x_{i},w_{i})\\}_{i=1}^{k})$
$\displaystyle=g_{\max}(\\{x_{i}\\}_{i=1}^{k})$ (4.14) $\displaystyle
g_{\operatorname{lse}}(\\{(x_{i},w_{i})\\}_{i=1}^{k};r)$
$\displaystyle=g_{\operatorname{lse}}(\\{x_{i}\\}_{i=1}^{k};r)$ (4.15)

###  4.3 Missing data

One detail that was deliberately left out in the description of the HMill
framework is how it handles incomplete or missing data. This phenomenon is
nowadays ubiquitous in many data sources and occurs due to various reasons—a
high price of obtaining an observation, information being unreachable due to
privacy reasons, a gradual change in the definition of data being gathered, or
a faulty collection process, to name some examples. Regardless of the origin
of missing data, it is wasteful to throw away the incomplete observations
altogether. Possibly missing labels gave rise to fields such as supervised,
semi-supervised and unsupervised learning. The HMill framework helps to
leverage knowledge hidden in incomplete observations, not labels. Thanks to
the hierarchical structure of both samples and models, we can still represent
samples with missing information fragments at various levels of abstraction.
Problems of this type can be categorized into three not necessarily separate
types: Figure 4.3: Example of a sample tree following schema in Figure 3.7
with missing data. In this particular case all three possible types of missing
data appear. 1. Missing information fragments in array nodes 2. Empty bag
nodes 3. Whole data source missing in product nodes One example of a sample
tree with missing data is in Figure 4.3. In this particular example,
information fragment in an array node responsible for modelling stems (olive)
is not observed, which is denoted by the $\varnothing$ sign. The subtree
modelling blooms (navy) is a bag node with no instances observed. The third
source of data—buds (green)—is missing altogether. According to our
experience, the latter two cases are more common than the first one. To
incorporate missing data into the framework, we provide model nodes with
additional functionality to handle these situations and also declare sample
trees with missing data to match the same schema as their counterparts without
incomplete information so that they can be processed with same structured
models. Although this requires a slight modification of Definitions 1, 2, 3
and 4, we will not delve into details and will instead explain, how models can
be extended to deal with missing data. Recall that any HMill model node can be
interpreted as a function from a set of sample trees following a certain
schema $\mathcal{T}_{s}$ to $\mathbb{R}^{m}$, where $m\in\mathbb{N}$. This
means that during evaluation, each model node always receives vectors of the
same size from all of its children. To deal with missing data, we introduce a
vector of real parameters $\bm{\psi}$ into each model node. For an array model
node implementing a mapping $f\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$,
$\bm{\psi}$ has $m$ elements. During both training and evaluation of the
model, whenever the array model node encounters a missing information fragment
$\varnothing$, instead of attempting to evaluate $f$, it returns vector
$\bm{\psi}$. For bag model nodes, $\bm{\psi}$ has the same dimensionality as
the output of aggregation $g$. If a bag is empty, we feed to bag model $f_{B}$
the vector of default values $\bm{\psi}$ instead of aggregated instances. In
the case of product model node, we keep one vector $\bm{\psi}_{i}$ for each of
its children. Whenever is the $i$-th child missing, we provide the
corresponding $\bm{\psi}_{i}$ as input to concatenation and further
processing. Since all these parameters are learned, this technique enables
effortless processing of both samples containing all information and samples
with missing data, all at the price of a tiny increase in the overall number
of parameters. Note that even though the three distinguished cases of missing
data may seem equivalent, in fact, each of them represents different
situations and results in a different treatment. For example, consider the
sample from Figure 4.3, but at this time missing the middle branch altogether.
This would result in filling in default values from $\bm{\psi}_{2}$ at the
level of the root. However, in the current situation, we use $\bm{\psi}$ from
the bag model node, which means that different default values are used and
also in this case they are futher processed by bag model $f_{B}$. Therefore, a
small change in the sample results in an entirely different computation.
Parameters $\bm{\psi}$ can be either fixed or even learned during training.
This requires a procedure for computing gradients of the loss function for all
learnable parameters in the model. This is the main topic of the following
section.

###  4.4 Gradient computation

As all functions composing models are defined to be (piecewise)
differentiable, and the computational graph is a tree, gradients with respect
to all parameters can be computed by adaptation of the original
backpropagation algorithm [Rumelhart et al., 1986], which we describe in
Algorithm 1. It proceeds in a preorder fashion and for each model node in the
tree computes a derivative of loss with respect to the parameters of the node
and also to its input. The derivative with respect to the input is used
subsequently by the children of the node. In step 8, we use a combination sum
and chain rule for aggregation functions and the standard backpropagation
algorithm for neural network layers to backpropagate through a given model
node. Note that this pseudocode is supposed to give a rough idea of how
gradients are computed and some details are omitted. For instance, because all
instances in a bag data node are processed with the same instance model in a
corresponding bag model node, the parameters are shared. As a result,
gradients of all submodels that instance model consists of have to be
appropriately aggregated over all instances in the bag, for example, averaged
in the case of $\operatorname{mean}$ aggregation. Also, if default vectors
$\bm{\psi}$ are used, Algorithm 1 has to be modified to change the control
flow according to whether the default values have been filled in or not.
Delving further into details is beyond the scope of this thesis and we refer
the reader to the literature about differentiable programming and our
implementation.

[1] procedure Backprop($s$, $F$, $\mathcal{L}$, $t$)  a schema $s$, a model
$F\colon\mathcal{T}_{s}\to\mathbb{R}^{n}$ constructed from $s$, a loss
function $\mathcal{L}\colon\mathbb{R}^{n}\to\mathbb{R}^{+}_{0}$, a sample
$t\in\mathcal{T}_{s}$ $\frac{\partial\mathcal{L}}{\partial\theta}$ for all
model parameters $\theta$ evaluate $F(t)$ and save inputs $i_{f}(t)$ and
outputs $o_{f}(t)$ and of every model node $f$ in $F$ compute
$\frac{\partial\mathcal{L}}{\partial F(t)}$ $q=\\{\operatorname{root}(F)\\}$
$\triangleright$ push the root model node to the queue  while $\lvert
q\rvert>0$ do $f\leftarrow\operatorname{pop}(q)$ $\triangleright$ pop one
model node from the queue  $\Theta_{f}\leftarrow$ all parameters of $f$
compute $\frac{\partial\mathcal{L}}{\partial i_{f}(t)}$ and
$\frac{\partial\mathcal{L}}{\partial\theta}$ for all $\theta\in\Theta_{f}$
using $\frac{\partial\mathcal{L}}{\partial
o_{f}\mathopen{}\mathclose{{}\left(t}\right)}$ computed in the parent
$\forall c\in\operatorname{children}(t)\colon\operatorname{push}(q,c)$ end
while end procedure

Algorithm 1 Adaptation of original backpropagation algorithm [Rumelhart et
al., 1986] to a tree-structured computational graph.

###  4.5 Batching procedure

Nowadays, full Batch gradient descent, which computes gradients of model
parameters using the whole dataset in each iteration, is no longer employed
for many problems. The reasons are twofold. Firstly, iterating over all
training samples requires a lot of time and is also considered inefficient
given the fact that many samples are very similar. Secondly, it is known that
stochasticity in training procedure improves performance. Instead of full
Batch gradient descent, its stochastic variants are used—Stochastic gradient
descent, which selects777The selection is performed either by sampling or some
deterministic procedure, such as iteration over the dataset. one observation
for every update, or Mini-batch gradient descent, which interpolates between
these two, by selecting a predefined number of observations from the dataset.
We have already defined how a single HMill sample looks like. In this section,
we describe how to group samples into (mini)batches and efficiently compute
gradients with respect to the whole (mini)batch. Two prominent guiding
principles in modern scientific computing are array programming and a static
computational graph. Array programming refers to writing computer programs in
a way that facilitates the application of operations to multiple values at
once, so that specialized vector instructions in CPUs can be utilized.
Regardless of whether the data are processed on GPU or on vector processors,
the focus is to group (possibly independent) data together and perform
operations on it simultaneously. For instance, computation of an output of a
neural network can be expressed as a composition of matrix-vector
multiplications and elementwise application of activation functions. If
matrix-matrix multiplications are performed instead, a whole batch of samples
can be run through the network at once. This is much faster than feeding
samples to the network one by one. The second principle is to keep the
computational graph as static as possible so that both high-level
optimizations at the level of computational graph and low-level optimizations
in the compiler can be utilized. A naïve approach to compute gradients with
respect to each sample tree from the batch first and summing the results
violates both of these principles. Nevertheless, HMill offers a way to group
sample trees following the same schema into a single tree representing the
whole batch in a manner respecting both principles. In other words, as in the
standard neural network case, all matrix-vector operations can be transformed
into matrix-matrix operations, and the computational graph is identical. The
high-level overview of the whole procedure is that all sample trees from the
batch are first horizontally overlapped and the data nodes in the same
positions in the tree are merged into one. This is possible due to the fact
that all samples match the same schema and requires small extensions of all
data node types so that they can carry more than one sample. All model nodes
will instead of one vector per sample return a matrix per multiple samples, in
which outputs for samples are stored in columns.

####  4.5.1 Array data nodes

In the case of an array data node, the modification is done in the same way as
in the aforementioned example, that is, instead of one vector representing the
information fragment of the sample in the given position in the tree, it
carries multiple vectors concatenated into a matrix for multiple samples. This
is illustrated in Figure 4.4, where three array nodes representing three
different samples are merged together into one. Given an array data node and
its matching model node, evaluation is performed in the same way as in Figure
3.9, however, this time we input a matrix into function $f$. As a result, we
obtain not one vector, but a whole matrix, where each column corresponds to
one sample.

####  4.5.2 Bag data nodes

Figure 4.4: Illustration of merging three array data nodes representing three
different samples (drawn by different colors) into a single array node.
Merging several bag data nodes together is possible thanks to the following
observation. Because all samples in the batch follow the same schema, all
instances from all samples in the batch are represented by trees that also
follow the same (sub)schema. This implies that all trees representing
instances can be collapsed into a single tree as well recursively by following
the same algorithm, and this observation holds for any position of the bag
node in the schema tree. Hence, given a batch of samples, to merge a
corresponding bag data nodes together, we collect all instances from all
samples in the batch and apply the same algorithm. Put differently, the bag
node now has only one child, which is a result of merging all instances of all
samples it stores. Moreover, we remember the indices of instances in order to
assign them to bags during evaluation. For a specific example, refer to Figure
4.5, where a process of merging three bag nodes with two, one and three
instances in their bags is illustrated. As all instance trees (in this example
trees consisting of only one array node) follow the same schema, they can be
collapsed into one tree. In this case, this means that instead of six
different array data nodes representing six different instances in three bags,
we have one array node, which becomes the only child and stores all
information. Note, that this child can also be another arbitrarily complex
subtree—bag data node or product data node, depending on the schema, however,
thanks to the observation we made, it is always possible to merge all
instances of all samples into a single subtree. Therefore, the resulting bag
data node always has only one child after merge. The evaluation of a merged
node with a given bag model node is performed analogically to the single
sample case (Figure 3.10). We first use the instance model $f_{I}$ on the only
child obtaining a matrix, where each column stores a result of $f_{I}$ on the
corresponding instance. In example in Figure 4.5, we would get a matrix with
the number of rows equal to the dimensionality of target space of $f_{I}$ and
six columns. Then, aggregation $g$ is performed using the indices stored in
the merged bag node, so that instances belonging to the same bags are
correctly aggregated together. We store the results into a narrower matrix
having as many columns as there are samples (in this case three) and apply bag
model $f_{B}$ on it, which may change the number of rows (the dimensionality
of target space of $f_{B}$), but the number of columns stays the same. Figure
4.5: Merging three bag nodes corresponding to different samples (drawn by
different colors). See text for details.

####  4.5.3 Product data nodes

Finally, if all samples in the batch follow the same schema, this also means
that all product data nodes at the same position in different sample trees
have children that can be merged into one subtree for all samples. In other
words, even though children of a single product data node may match different
(sub)schemata and therefore they may not be combinable as is the case with bag
nodes, we can merge children across multiple sample trees provided that they
follow the same (sub)schema. This is done by merging all first children into a
single tree, then all second children into a single tree, and so forth. See
Figure 4.6, where this procedure is applied to product data nodes from three
different samples with three children each. Children are grouped according to
the (sub)schema they follow888or equivalently by their order in product data
node, since we have defined product nodes to impose ordering on their children
and all groups are recursively collapsed into one subtree. The evaluation of a
product model node is done in a similar way as in Figure 3.11. We apply
submodels $f_{1},\ldots,f_{l}$ to children and obtain a matrix with as many
columns as there are samples for each child. In our specific example, we would
obtain three matrices with three columns and possibly different number of
rows. Then, all matrices are vertically concatenated and mapping $f$ is
applied, which leads to one column for each sample on the output. Figure 4.6:
Procedure of merging three different product nodes matching same schema in
one. Here, the fill colors indicate different (sub)schema that children follow
(see Figure 3.5) and the same border colors represent subtrees coming from the
same sample.

Applying this recursive procedure on any number of samples leads to a single
tree, where all intermediate representations are stored in one matrix. Static
schema across the batch guarantees static computational graph, and therefore
the computation of gradients on the resulting tree is effective. To simplify
the explanation, we have assumed that there are no missing data present in any
sample. In practice, this is not true, and slightly more complex functionality
is required. For instance, if any of the bags is empty when merging bag data
nodes together, we use the methodology described in Section 4.3. Hence, the
resulting bag data node remembers that this particular bag has zero instances
and when performing aggregation, the resulting column is filled with the
default values.

Having read this and the previous chapter, the reader should now be familiar
with sample representation and model definition in HMill, as well as the
technical details of its implementation. Should the reader wish to learn even
more, we recommend referring to our implementation [Pevný and Mandlík, 2018].
The source code of version 1.2.0 of the framework is attached to the
submission. In the remaining chapters, we delve into the practical adoption of
HMill.

##  Chapter 5 Towards automated processing of structured data

In this chapter, we demonstrate how the HMill framework lends itself to the
processing of structured hierarchical data. We mentioned in Chapter 1 that a
majority of ML algorithms is applicable only to fixed-size input, even though
real-world data is hardly ever found in a form in which each observation is
described by a vector of fixed length. On the contrary, plenty of data appears
naturally structured as graphs or trees, with information being specified not
only in nodes but also in the structure itself. This brings higher flexibility
when opting for formats in which the data is stored, however brings new
challenges to methods processing such input. We focus on such hierarchical
formats in this part, and in the remaining Chapters 6 and 7 on the processing
of graphs. Figure 5.1: Syntax diagram of a JSON value. Specifically, we
consider hierarchically structured data formats ubiquitous in today’s Internet
communication, such as JSON and XML. We will demonstrate how these formats
seamlessly integrate into the framework and therefore how data stored in this
way can be learned from. We will illustrate this on JSON format in particular,
however, everything applies to XML s as well. JSON is a lightweight, human-
readable format designed for data exchange independently of the programming
language used. The standard111https://json.org defines three main structures
(data types): Values are the most general type of JSON structure and can be
either primitive data types, such as strings, numbers, booleans or nulls, or
composite object or array types. Objects are unordered collections of zero,
one, or more key-value pairs, which abstract on dictionaries, hash tables, or
associative arrays in a modern programming language. Keys are defined as
arbitrary strings and values are JSON values. Arrays are, contrary to objects,
ordered lists of values, corresponding to standard arrays in programming
languages. Elements of JSON arrays can be any JSON values.

(a)

(b)

Figure 5.2: Syntax diagram of a JSON object (5.2(a)) and a JSON array
(5.2(b)). Figure 5.3: Example of a JSON document. Apparently, JSON format
provides higher expressive power for objects of interest, as the single-vector
representation can be encoded as JSON array with real numbers as its elements
and constant size across observations. Whereas JSON arrays are most of the
time used to store elements with the same interpretation222JSON standard
though allows elements to be any JSON value., JSON object is a data type
intended to accommodate multiple possibly heterogeneous attributes with an
arbitrary representation. Needless to say, thanks to the recursive definitions
of JSON data types, both arrays and objects can be further nested, creating
complex composite structures capable of capturing many real-world phenomena.
An example of a valid JSON file storing information about an IoT device is in
Figure 5.3 and it is based on the task we attempt to solve in the experiment
section of this chapter. The whole top-level object consists of data entries
specifying the MAC address and the IP address of the device, services it
provides and the class of the device, which in this case serves as a label. As
one device can expose multiple ports and use different protocols for
communication, we model this as an array of objects, each specifying one
service. Every JSON value can be regarded as a node in a tree representation
of the document, with arrays and objects having several children and primitive
values, such as numbers or strings, being leaves of the tree. Schemata
representing the graphical structure can be defined in a similar way as in
HMill, and therefore, JSON documents can be validated against schemata. In
this thesis, we focus mainly on the classification task with valid JSON
documents or similar data formats as input. Traditional approaches include
rule-based systems or flattening. Rule-based systems combine manually designed
simple rules based on domain knowledge with transformations of suitable
subsets of data. Referring to Figure 5.3, one example of such rule could be
that the observed device is a phone as long as it uses the UDP protocol on
port $5353$. Hence, rules, such as logical AND or OR help with dealing with
JSON arrays of arbitrary length. To increase the space of such rule-based
hypotheses a transformation on the appropriate subsets of the document can be
defined or even learned. For instance, in the previous example, we could
collect all port-protocol pairs in documents describing phones and infer a
detector returning a score specifying how likely this pair appears in the
descriptions of phone devices. Then, a device would be classified as a phone
if all of its elements in "services" array reach score above a certain
threshold. The rule-based approach to handling tree-structured data in many
ways resembles the instance-space paradigm for solving MIL problems from
Section 2.1 with all of its benefits and drawbacks. The main advantage of the
rule-based approach is that rules are easily interpretable, however, they
require deep domain knowledge and can be very difficult to design properly.
Moreover, in the presence of the concept drift, rules may become quickly
obsolete, and sometimes, the language of simple rules may not even be
expressive enough to reasonably discriminate samples, especially for more
complex structures. The flattening approach deals with hierarchical data by
defining a procedure that yields a vectorized representation of any given
tree. This is convenient since a lot of well-researched ML algorithms can be
used with the resulting vector representation. One major disadvantage of this
approach is that it does not deal with complex variable structures in a
principled way—as different tree-structured observations may contain different
amount of information (for instance by storing JSON arrays of different
length), transformation to a fixed-size vector can be either lossy if the size
is too small, or wasteful if it is too high. The hierarchical nature of tree-
structured data presents an opportunity to increase the capacity of models by
leveraging structural information from the data while keeping the number of
parameters sufficiently low. The prior art on models learning from both raw
data in leaves and topologies of such structures is scarce and existing
methods use limiting assumptions, such as all leaves being present in the same
depth [Hammer et al., 2004, Zhang et al., 2018a], or the tree having a very
specific structure, such as syntactic parse trees in natural language
processing [Irsoy and Cardie, 2014, Kokkinos and Potamianos, 2017, Zhang et
al., 2018b]. To the best of our knowledge, the only fully general approach to
end-to-end learning on hierarchical data is proposed in [Woof and Chen, 2020].

###  5.1 Modelling JSONs with HMill

To use HMill for learning from JSON documents, we need to make one mild
assumption that all documents in the dataset must be structured in the same
way (in HMill parlance, they must follow the same schema). Given a set of JSON
documents $D$, this assumption, loosely speaking, translates to the following
conditions: • Keys of every JSON object located at a certain position in the
tree representation of any document $d\in D$ must come from a finite ‘super’
set. In other words, given a fixed position in the tree, we know which keys
can be encountered in this particular position before training procedure
starts. • Values corresponding to the same key of JSON objects located at the
same position must have the same structure for every document $d\in D$. • JSON
Arrays must be either empty or contain elements with the same structure. For a
certain position in the tree, elements of arrays of different documents from
$D$ must also have the same structure. The whole situation is analogical to
the notion of HMill schema, which is defined and discussed in Section 3.2.1.
To demonstrate this condition on the example from Figure 5.3, if we wanted to
train our model on a set of multiple similar JSON documents $D$, we would
require each one of them to contain at the top level only key-value pairs with
keys from set $\\{$ "mac", "ip", "services", "device_class" $\\}$, with every
value corresponding to any key apart from services being represented as a
string of characters. Note that this does not imply that all documents from
$D$ must have the same key sets—some keys might be missing (see Section 4.3).
All documents in $D$ would further have either no services key present, or its
value would be an array of another type of JSON objects with keys from set
$\\{$ "port", "protocol" $\\}$ with values represented as an integer for ports
and string for the protocol. Even though the assumption of the common schema
for all documents might seem limiting, according to our experience, this is
rarely the case. Despite the JSON standard being rather forgiving in how the
documents can be structured, when we observe documents from one specific data
source or a collection of data sources, they usually follow the same schema
directly, or only a trivial transformation is required. Figure 5.4: A
transformation of JSON document from Figure 5.3 into a valid HMill sample. See
text for details. Translating a JSON document to an HMill sample is
straightforward. Each leaf of the tree with primitive data type entries such
as numbers, booleans, or strings is transformed into an array data node
$\operatorname{an}(h(\ldots))$, where $h$ is a suitable mapping to some real
vector space (see Section 3.2.1). Therefore, the entry is converted to a
vector of real numbers. Note that mapping $h$ can be defined differently for
leaves at different positions in the tree, however, it is the same for leaves
at the same position in the tree across all documents. For instance, mapping
$h$ in an array node modelling IP addresses is the same for all documents but
may differ from a mapping used for MAC addresses. To model JSON objects, we
use product data nodes. Recall that given a particular position in the tree, a
finite set $K$ of all keys used by JSON objects at this position is known.
First, an arbitrary ordering $\prec$ is imposed on all elements $k\in K$, for
instance, trivial lexicographical order of key strings. This way, we can
assign a unique number from $\\{1,\ldots,|K|\\}$ to each key, according to its
rank with respect to $\prec$. This number is regarded as the index of the
child (subtree) corresponding to the JSON value stored with this key in the
product node. This is merely a technical detail, stemming from the fact that
children of product data nodes can be semantically heterogeneous, however,
key-value pairs from JSON objects are originally unordered. Therefore, we need
to define this arbitrary, yet precise mapping from keys to children indices.
As children of product data nodes do not have to match the same schema, we can
model heterogeneous values stored under different keys with an arbitrary
substructure. For instance, each of the mac, ip and services entries could be
modelled by a different subtree and therefore processed by submodels with
variable topology. Finally, JSON arrays are modelled as HMill bag data nodes
by interpreting each homogeneous entry in the array as an instance in the
corresponding bag. Note that HMill handles empty arrays as well as arrays of
any length. The whole transformation is depicted in Figure 5.4. The root of
the resulting sample tree is a product data node, with the first two children
being subtrees for modelling IP and MAC addresses. The last child is a bag
data node representing an array of elements from the services field, each of
which is another JSON object and is therefore modelled as another product data
node. Four distinct mappings $h_{1}$, $h_{2}$, $h_{3}$ and $h_{4}$ are defined
for modelling different information fragments (IP and MAC addresses, port and
protocol). For specific examples of various mappings $h$, refer to Figure 3.3.
Also, as JSON object is unordered, the order of keys in the second element of
services arrays is changed so that each product node modelling entries in the
array has submodel for ports as the first child, and submodel for protocols as
the second child. The resulting pipeline for learning from JSON documents is
not much different from the pipeline in the current machine learning. We first
collect the data, in this case, a collection of documents $D$ and deduce the
schema. This can be done in an automated way, for instance with JsonGrinder.jl
library [Pevný and Račinský, 2019]. Then, a suitable model is defined
according to the schema and trained in an end-to-end fashion on both primitive
types in leaves of the documents and the whole topology. As models learn from
the structure as well, a need for the tedious, trial-and-error design of a
flattening procedure to encode the structure in a fixed-size vector is
alleviated. Also, learning HMill models from the data is more general than
rule-based approaches, and a more complex hypothesis can be learned, as
aggregation functions in bag model nodes can easily substitute hand-written
rules. We believe that HMill will enable out-of-the-box learning from tree-
structured data expressed as JSON s, XML s and other similar formats in tasks
across variable domains and will help other researches with processing data so
far considered incompatible with current methods.

###  5.2 IoT device identification (use case)

The first of three use cases we will present throughout the thesis concerns
IoT domain. IoT is a term referring to a network of billions of existing
devices in the Internet that communicate with each other and collect and share
data through special sensors. Thanks to advances in multiple disciplines and
cheaper manufacturing, the IoT network saw immense growth in the last decade
and, for instance, nearly half of the North American households have an
Internet-connected television or a streaming device [Kumar et al., 2019]. It
is believed that IoT will benefit society across sectors, ranging from health
services and home automation for the elderly and disabled, to energy savings
in ‘smart homes’. However, The progress is hindered by security and privacy
concerns and rightfully so—some of the examples include up to $600\,000$
infected devices being used in massive Distributed Denial-of-Service (DDoS)
attacks by the infamous Mirai botnet [Antonakakis et al., 2017], or smart home
devices being compromised [Hernandez and Buentello, 2014, Kumar et al., 2018].
Thus, it is of utmost interest of security companies and device vendors to
turn their attention to the security of the IoT network and develop suitable
solutions. In [Kumar et al., 2019], researchers from University of Illinois
Urbana-Champaign, Stanford University and Avast Software cooperated to conduct
the first large-scale investigation and empirical analysis of IoT devices in
homes that are otherwise invisible in Internet scanning. Data was collected by
local subnet scanning manually initiated by users. One of the problems that
emerge when the network is scanned is an identification of the device type
(TV, wearable, voice assistant, …) using only information obtained from the
scan. In the following sections, we demonstrate how HMill helps with solving
this task.

####  5.2.1 Data description

The scan is performed by a tool called WiFi Inspector built into the antivirus
client.333All data about devices is collected only for research purposes,
sensitive information is anonymized or aggregated and users are informed and
can opt out anytime. See [Kumar et al., 2019] for further details about
privacy policies. First, the list of candidates is generated using the ARP
table and ARP, SSDP and mDNS scans of the network. Then, additional knowledge
about each candidate is gathered from both the network and the application
layer. Specifically, ports used for listening and protocols for communication
are further scanned, as well as DHCP, mDNS, SSDP and UPnP broadcasts made by
the device. Two specific examples are shown in Apendix B (Figures B.3 and
B.2). Figure 5.5: Distribution of class labels in the training set of the
dataset obtained from Kaggle. All ratios are stated in percents. We worked
with two specific datasets. The first dataset is publicly accessible and comes
from the Kaggle
challenge444https://www.kaggle.com/c/cybersecprague2019-challenge/overview. It
contains $57\,906$ documents for training and $77\,777$ documents for testing
grouped into $\lvert\mathcal{C}\rvert=13$ different classes. The schema of the
training data as obtained from JsonGrinder.jl can be found in Figure B.1 in
Appendix B. All fields except of device_class, which contains the label, and
device_id, which is a unique identifier of the sample, can be used as input to
the model. Describing the data in full detail is beyond the scope of this
thesis, and we refer interested readers to the Kaggle competition description
or the cited paper. A pie chart of the class distribution in the training data
is in Figure 5.5. The second dataset containing $8\,816$ documents for
training, $99719$ documents for testing, and $13$ classes was kindly provided
by Avast. There is more information in its documents, such as TLS cipher
suites list or HTTP user agents, therefore the schema of this dataset is
slightly different. Similarly to other problems tackled in the computer
security domain, identifying a device type from JSON documents of this schema
is trivial for humans most of the time, however, designing a system for
automatic recognition is a challenging task. Device labels can often be simply
inferred from information in a single string field in a JSON document, for
instance, manufacturers often inform about the model name and device
specifications in responses to the aforementioned probes. Nevertheless, this
is not always the case and for some devices, their type can be only deduced by
inspection of their overall communication interface and behavior. Moreover,
the domains of information fragments are typically discrete (strings,
categorical variables, or natural numbers) and many fields may contain an
empty string or a key-value pair may not be specified at all. This is due to
protocols being unsupported on the device side, timeouts or failures in the
network or other reasons. Finally, most of the information from JSON document
is usually redundant and not required to make the correct decision.

####  5.2.2 Experimental results

Learning to classify IoT device type from JSON s is rather straightforward
using the process described in Section 5.1. We inferred a schema of the
available documents and defined a hierarchical HMill model. Mapping $h$ in
each array data node, where fragment space $\mathcal{F}$ contained strings of
characters, was defined to compute histograms of tri-grams, in cases, where
there was only a handful of unique string values (for instance ports,
protocols, or manufacturers), we interpreted the field as a categorical
variable and used one-hot encoding, and in the remaining cases (integers,
reals), $h$ was set to identity. To promote effortless processing of sample
trees following an arbitrary schema, Mill.jl provides an auxiliary procedure
that helps with defining complex composite models. As arguments, it accepts a
schema and a prescription specifying topology of neural network components in
HMill models (mapping $f$ in array model nodes, $f_{B}$ in bag model nodes,
and $f$ in product model nodes) and types of aggregation functions $g$. For
simplicity, we opted for all neural networks components to consist of one
feedforward layer with output dimension $50$, dropout regularization
[Srivastava et al., 2014] with probability $0.5$ and the $\tanh$ activation
function. In every bag model node, we used a concatenation of all four
functions discussed in 4.2. Note that specifying output dimensions for
mappings $h$ in array data nodes as well as for neural network components is
enough to specify the whole model tree, since the input dimension of every
component can be deduced. The resulting model therefore mapped every JSON
document to a vector with dimension $50$, which was followed by one final
layer to obtain a vector of size $\lvert\mathcal{C}\rvert$ and
$\operatorname{softmax}$ transformation to obtain class probabilities
conditioned on the input.

In each experiment, the model was trained for $30$ epochs (iterations over the
whole training set) using minibatches of size $100$ and Adam555with the
default parameters $\alpha=0.001$, $\beta_{1}=0.9$, $\beta_{2}=0.999$ and
$\epsilon=10^{-8}$ optimizer [Kingma and Ba, 2014]. For each neural network
layer, the Glorot normal initialization [Glorot and Bengio, 2010] was employed
for weights, and all biases were initialized to zero. For a loss function, we
used the standard multiclass cross entropy (also known as multinomial logistic
loss) on the conditional probabilities output by the model:
$\mathcal{L}(\theta)=-\frac{1}{n}\sum_{i=1}^{n}\sum_{c\in\mathcal{C}}\mathbbm{1}\mathopen{}\mathclose{{}\left\\{c=c_{i}}\right\\}\log
p(c|d_{i};\theta)$ (5.1) where $n$ is the number of samples (documents) in the
batch, $d_{i}$ a representation of $i$-th document as an HMill sample, $c_{i}$
its class label, and $p(c|d_{i};\theta)$ the conditional probability of $i$-th
document having label $c$ according to the model parametrized by $\theta$.
After the training was finished, we measured the accuracy on the training set
as well as averaged F1 score over all classes. Measured over $10$ runs with a
different random seed, we reached accuracy $0.9709\pm 0.003$ 666$\mu\pm\sigma$
notation specifies the average and the standard deviation over all runs. and
F1 score $0.9582\pm 0.007$ for the Kaggle dataset, which is on par with the
winning entry. Metrics measured on the second dataset are in Table 5.1. We
compared the results to the system described in [Kumar et al., 2019], which is
deployed in Avast production. It consists of more than a thousand rules
written by experts and a machine learning-based predictor, each of which
returns a class prediction together with a confidence score. The overall
prediction is then derived by combining predictions of individual components.
The ML-based part is implemented as an ensemble of classifiers as described in
the paper. This ensemble was recently replaced by a single model—one-vs-rest
multiclass adaptation of soft-margin SVM classifier [Cortes and Vapnik, 1995],
with Radial Basis Function (RBF) kernel ($\gamma=1/(2\sigma^{2})=10^{-4}$) and
regularization constant $C=10^{5}$. It uses a manually designed set of
features extracted from JSON document, mainly from UPnP and mDNS. In the
table, we present the accuracy of the whole system, that is, expert rules
together with an ML classifier, for both ensemble and SVM-based
implementation. Table 5.1: Metrics measured on the testing set from the
dataset obtained from Avast. method accuracy F1 score ensemble $0.878$ - SVM
$0.925$ - HMill $\mathbf{0.969\pm 0.004}$ $\mathbf{0.924\pm 0.006}$ It needs
to be emphasized that due to the fact that datasets represent merely a
fraction of devices in the world, the results presented here are likely
slightly skewed in comparison to the performance measured in the production
environment. Furthermore, accuracy may not be the most suitable metric if
there are different costs of misclassification for labels, or if we are most
interested in a subset of high-confidence predictions and do not want to
measure accuracy on ‘hard’ samples. The fact that simple baseline HMill models
trained on raw data without any fine-tuning achieve competitive performance on
both datasets shows great promise. We assume that by further experimenting
with model components (for example adding layers to submodels in model nodes),
even better results are attainable. This, however, is not the goal of this
thesis in which we aim to demonstrate the versatility of the framework on
different problems.

##  Chapter 6 HMill for graph inference

The previous chapter showed how the HMill framework can be used to learn from
tree-structured data together with one particular use case concerning solving
the device identification problem. Yet, one of the main requirements outlined
in Chapter 1—dealing with the violated i.i.d. assumption—remains unaddressed.
As it was discussed before, independence between individual samples is in some
tasks rather unrealistic assumption, especially in the computer security
domain. Two remaining use cases of the framework in this thesis are prime
examples of this. The first one concerns malicious binary file classification
at the end-user’s machine from a ‘snapshot’ of the state of the operating
system. In the last use case, we explore possibilities for second-level domain
classification based on its behavior in the network. Even though both tasks
are essentially different problems with different data inputs, they have one
property in common—in the system, there are many objects of interest of the
same type (binaries in the former case and domains in the latter), all of
which contribute to the developments in the surrounding system. To give a
specific example, in the operating systems many binary files are being run at
once, some of them multiple times in different processes. Since processes may
further install other binary files and spawn or terminate other processes, the
whole operating system is interconnected in a complex way. To infer the
maliciousness of a binary, we should take into account the maliciousness of
other binaries or processes in the system it interacted with. Nevertheless,
the maliciousness of other objects is unknown as well and deducing it
corresponds to essentially the same task. Thus, the i.i.d. assumption does not
hold. In this chapter, we first briefly revise statistical frameworks and
methodologies for representing such complex real-world systems and then
demonstrate how HMill is used to perform inference in these systems. The final
part of this chapter presents experiment results from the task of binary file
classification based on inference in a graph representation of the operating
system.

###  6.1 Prior art

In the following text, we denote
$X=\mathopen{}\mathclose{{}\left\\{X_{1},X_{2},\ldots,X_{n}}\right\\}$ a
system of $n$ objects we want to model, where each $X_{i}$ is regarded as a
random variable taking values from some domain $\mathcal{X}$, which depends on
the problem. For instance, if we want to classify objects into several
classes, we define all variables to be binary
$X_{i}\in\mathopen{}\mathclose{{}\left\\{0,1}\right\\}$ in the case of binary
classification, discrete-valued
$X_{i}\in\mathopen{}\mathclose{{}\left\\{1,\ldots,\lvert\mathcal{C}\rvert}\right\\}$,
where $\mathcal{C}$ is a set of classes, for multi-class classification, and
real-valued $X_{i}\in\mathbb{R}$ for regression tasks. Under the independence
assumption, it is enough to consider only marginals $p_{i}(X_{i})$ since the
joint probability can be eventually computed using the product rule as
$p(X)=\prod_{i=1}^{n}p_{i}(X_{i})$. When this assumption does not hold, we
have no choice but to model the whole joint probability $p(X)$, so that all
dependences between variables can be captured. It is often not the case that
all random variables in $X$ are mutually dependent. This happens when $X$
consists of several independent subsets of variables, or when there are pairs
of variables that are conditionally independent111Random variables $X$ and $Y$
are said to be conditionally independent given a set of variables $Z$ if
$p(X|Y,Z)=p(X|Z)$, or equivalently, $p(X,Y|Z)=p(X|Z)p(y|Z)$. of each other
given some (smaller) subset of variables in the model. Consequently, the joint
probability distribution $p(X)$ can be factorized into a product of several
terms (usually called factors), each computed from a subset of all variables:
$p(X)=\prod_{S\subseteq X}f(S)$ (6.1) This is a generalization of the case
when i.i.d. assumption holds, as we can have one factor for each of the
variables equal to its marginal probability $p_{i}(X_{i})$.

####  6.1.1 Graphical models

Graphical models, or Structured probabilistic models, is an umbrella term for
all models, where the possible factorization is encoded into a graph, whose
set of vertices are random variables and an edge between two random variables
indicates dependence. Therefore, the more (conditionally) independent random
variables are, the more efficient their representation as a graphical model
is, compared to a full description of probabilities in a table. According to
the exact definition of factors $f$ and the type of the graph, graphical
models can be further split into directed graphical models, referred to as
Bayesian networks [Pearl, 2013], and undirected models, known as Markov random
fields [R. Kindermann, 1980, Jordan, 2004]. As opposed to Bayesian networks,
where factors $f$ are given by conditional distributions and reflected in
edges of a directed graph as depicted in Figure 6.1, an analogical procedure
for undirected graphs is not that straightforward. Intuitively, two random
variables not connected by an edge should be conditionally independent of each
other given any separating subset.222The separating subset between two
vertices $x$ and $y$ in the graph is a minimal set of other vertices such that
$x$ becomes unreachable from $y$ after the removal of the separating subset
from the graph. The corollary is that a random variable is conditionally
independent of any other variable in the graph, given all variables in its
neighborhood. This idea is theoretically justified in the Hammersley-Clifford
[Hammersley, 1971] theorem, which states that every strictly positive
probability distribution can be factorized into a product of potential
functions defined over maximal cliques333Maximal clique is a maximal complete
subgraph (adding another node from graph breaks the completeness property). in
the graph and vice versa, every function represented as this product describes
a valid probability distribution. Therefore, a factorization of the
probability distribution in a Markov random field is proportional to the
product of non-negative potential functions $\psi_{C}$ defined on maximal
cliques of the graph: $p(X)\propto\prod_{C}\psi_{C}(X_{C})$ (6.2) where
$X_{C}$ is the set of variables in clique $C$. Potential functions do not have
to integrate to one and therefore are not proper probability distributions
themselves, however, their normalized product corresponds to a valid joint
distribution: $p(X)=\frac{1}{Z}\prod_{C}\psi_{C}(X_{C})$ (6.3) where
$Z=\int_{\mathcal{X}}\prod_{C}\psi_{C}(X_{C})\mathrm{d}X$ is a normalizing
constant, also referred to as a partition function.

(a)

(b)

Figure 6.1: Example of a Bayesian network and a Markov random field
representing a system of five variables. In (6.1(a)), there is a Bayesian
network representing a probability distribution factorizable as
$p(a,b,c,d,e)=p(a)p(b\>|\>a)p(c)p(d\>|\>b,c)p(e\>|\>d)$. In Bayesian networks,
edges signify direct conditioning, that is, $b$ is dependent on
$\mathopen{}\mathclose{{}\left\\{a}\right\\}$, $d$ is dependent on
$\mathopen{}\mathclose{{}\left\\{b,c}\right\\}$, and $e$ is dependent on
$\mathopen{}\mathclose{{}\left\\{d}\right\\}$. Additionally, $e$ is
conditionally independent of all other vertices given its parent $d$, and $d$
is conditionally independent of $a$ given $b$ and $c$. A representation of a
(different) probability distribution of five variables as a Markov random
field is in (6.1(b)). This distribution can be factorized as
$p(a,b,c,d,e)=\frac{1}{Z}\psi(a,b,c)\psi(b,c,d)\psi(b,e)$, where
$\\{a,b,c\\}$, $\\{b,c,d\\}$ and $\\{b,e\\}$ are maximal cliques in the graph.
For example, $e$ is conditionally independent of every other node given $b$,
and $a$ and $d$ are conditionally independent given their separating subset
$\\{b,c\\}$.  Assume we are given a joint probability distribution $p(X)$
represented as a graphical model and a set of observed variables with known
values
$X_{O}=\mathopen{}\mathclose{{}\left\\{x_{1},\ldots,x_{l}}\right\\}\subseteq
X$. Denoting the remaining variables to be inferred $Y=X\setminus
X_{O}=\mathopen{}\mathclose{{}\left\\{Y_{1},\ldots,Y_{k}}\right\\}$, one of
the relevant tasks is to find the most probable state
$Y^{*}=\operatorname*{arg\,max}_{Y}p(Y|X_{O})=\operatorname*{arg\,max}_{Y}p(Y,X_{O})$,
that is, to perform maximum a posteriori (MAP) inference. Another problem is
to compute marginals $p(Y_{i}|X)$ for each $i=1,\ldots,k$. For a specific loss
function $\mathcal{L}\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}$,
marginals can be then used to obtain the optimal Bayes classifier for
unobserved variables $Y$. Note that $X_{O}$ can be empty. For graphs without
loops, the canonical Belief propagation [Pearl, 1988] algorithm exists.
Unfortunately, exact computation of either of these probabilities and even
evaluation of constant $Z$ is a computationally intractable problem in the
general case when graphs may contain loops. In this case, approximation
methods such as graph cuts [Kohli and Torr, 2005] or Monte Carlo Markov chains
[Wang et al., 2000b] can be utilized. A prominent family of approximation
algorithms for graphical models is the family of so-called message-passing
algorithms containing the standard Loopy belief propagation [Frey and MacKay,
1998, Kschischang et al., 2001] algorithm also known as the sum-product
algorithm. In message-passing inference, messages representing local
statistical information are iteratively exchanged between vertices, which
update their beliefs accordingly, until some convergence criterion is met. The
key idea is that even though messages are computed from a small neighborhood
of each vertex, after a sufficient number of iterations, information is spread
across the entire graph and with appropriate message definitions, marginals
can be computed in each vertex separately. There are several drawbacks in
practical modelling and inference with graphical models. Namely, one needs to
suitably construct the graph and define potential functions, and algorithms
for approximate inference usually lack theoretical guarantees. In [Lin et al.,
2015, Zheng et al., 2015, Schwing and Urtasun, 2015] authors attempted to
learn the definition of potential functions, however, this required an
introduction of another assumption, for instance, fixed number of neighbors of
each vertex. Finally, it has been shown in [Ross et al., 2011] that Loopy
belief propagation and its variants can be regarded as special cases of a
general iterative scheme in which a sequence of interdependent local
classifications is performed. This suggests that not only potential functions
but also messages themselves, which are in Loopy belief propagation computed
using explicit formulas, can be learned. The notion of message-passing
inference, where information encoded in the graph is repeatedly collected from
neighbors, aggregated and further propagated, together with the key insight
that we can learn message predictors, are main motivations for performing
inference with HMill as explained later.

####  6.1.2 Graph neural networks

Since graphs are a natural representation of many real-world phenomena
occurring in biology, sociology, computer security, or finance, research
efforts have been turned to development of methods for effective learning from
such structures [Bronstein et al., 2016]. Graph neural network (GNN) s, first
proposed in [Scarselli et al., 2009], are an extension of classical
feedforward and recurrent architectures to graph input. The goal is to learn
an embedding vector $\bm{h}_{\mathcal{G}}$ for the whole graph or a set of
embeddings
$\mathopen{}\mathclose{{}\left\\{\bm{h}_{v_{1}},\bm{h}_{v_{2}},\ldots,\bm{h}_{v_{n}}}\right\\}$
for each vertex $v_{i}$ of the graph. Depending on the application, this
embedding is further processed. Graph Convolutional Networks (GCNs) [Kipf and
Welling, 2016] reflect the structural dependences together with features from
vertices by generalizing a convolution operation to the graph domain, leading
to the following update rule:
$\bm{H}^{(l+1)}=\sigma\mathopen{}\mathclose{{}\left(\bm{\hat{D}}^{-\frac{1}{2}}\bm{\hat{A}}\bm{\hat{D}}^{-\frac{1}{2}}\bm{H}^{(l)}\bm{W}^{(l)}}\right)$
(6.4) where activations $\bm{H}^{(0)}$ in the first iteration are initialized
with a matrix of features $\bm{X}$, $\bm{\hat{A}}=\bm{A}+\bm{I}$ is the
adjacency matrix of the graph with added loops, $\bm{\hat{D}}$ is a diagonal
vertex degree matrix, $\bm{W}^{(l)}$ is a matrix of parameters in $l$-th
layer, and $\sigma$ is a non-linear activation function. By stacking multiple
layers of (6.4), information is aggregated from the graph, and hidden states
of vertices are gradually refined. GCNs belong to the family of spectral
approaches, which utilize spectral representation of the graph as proposed in
[Bruna et al., 2013]. The disadvantage of spectral methods is that the
parameters are learned using a specific graph structure, and the model cannot
be applied once the input graph is even slightly perturbed. Spatial approaches
leverage spatiality in the graph and define graph convolution based on vertex
relations. GraphSAGE [Hamilton et al., 2017] is one of the first
representatives of spatial methods, which defines graph convolution as
follows:
$\mathbf{h}_{v}^{t}=\sigma\mathopen{}\mathclose{{}\left(\mathbf{W}^{t}\cdot
a\mathopen{}\mathclose{{}\left(\mathbf{h}_{v}^{t-1},\mathopen{}\mathclose{{}\left\\{\mathbf{h}_{u}^{t-1},\forall
u\in\mathcal{N}(v)}\right\\}}\right)}\right)$ (6.5) where $a$ is a
permutation-invariant aggregation function, and $\mathcal{N}(v)$ denotes the
neighborhood of vertex $v$. Note that in spatial approaches the central idea
is to iteratively update hidden state $\mathbf{h}_{v}$ at each node $v$ based
on its neighborhood, which corresponds to the message-passing scheme. Indeed,
in [Gilmer et al., 2017], authors reformulated a handful of existing GNN
architectures into a common general message-passing-based framework. In [Yoon
et al., 2019], it has been shown how GNN s can be utilized to learn the
messages from the original Loopy Belief Propagation algorithm and eventually
outperform it. Lately, research efforts have been put into modelling a
heterogeneous network, a graph with multiple vertex and edge types. As it is
the case with vanilla GNN s, to deal with this intrinsic property, many
different methods were proposed [Chang et al., 2015, Gui et al., 2016, Fu et
al., 2017, Shi et al., 2017, Dong et al., 2017, Zhang et al., 2018a, Sun et
al., 2018, Fan et al., 2019, Wang et al., 2019]. Prior art on GNN s is vast
and specialized approaches have been proposed for a particular task being
solved (vertex classification, graph embedding, et cetera), different graph
types and sizes, and whether the model should generalize to changes in the
graph structure. We merely scratched the surface of the family of GNN s and
heterogeneous network modelling to provide motivation for the HMill-based
approach to inference described in the next section. We refer the reader to
[Bronstein et al., 2016] for a review of methods operating on non-Euclidean
data and to [Zhang et al., 2018b, Zhou et al., 2018, Zhou et al., 2018] for
comprehensive reviews of GNN s.

###  6.2 Modelling heterogeneous networks with HMill

In this section, we demonstrate how HMill models can be utilized to perform
graph inference. In our examples, we focus mainly on the task of vertex
classification, however, the general idea can be easily extended to other
tasks. We show that models implemented in bag and product nodes are
appropriate for simulating a message-passing procedure, effectively adapting
the framework for modelling of not only tree-structured data but also complex
interactions between samples encoded in a graph. Ultimately, both vertices and
edges can be described not by a fixed vector of features, which is assumed in
GNN s, but by any of the data types processable by HMill. Again, this is
consistent with the fact that the main contribution of the HMill framework is
in its flexibility, which enables modelling of a wide range of data sources.
In the following text, we assume a graph
$\mathcal{G}=\mathopen{}\mathclose{{}\left(\mathcal{V},\mathcal{E}}\right)$ to
be a tuple consisting of a finite set of vertices (nodes) $\mathcal{V}$, and a
finite set of edges $\mathcal{E}\subseteq\binom{\mathcal{V}}{2}$ or
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ for undirected and directed
graphs, respectively. For a vertex $v\in\mathcal{V}$, we will denote its
neighboring vertices by
$\mathcal{N}(v)=\mathopen{}\mathclose{{}\left\\{u\mid\mathopen{}\mathclose{{}\left\\{u,v}\right\\}\in\mathcal{E}}\right\\}$
in an undirected graph, and in a directed graph, we will distinguish two types
of neighborhoods based on the orientation of edges:
$\mathcal{N}^{+}(v)=\mathopen{}\mathclose{{}\left\\{u\mid(u,v)\in\mathcal{E}}\right\\}$
and
$\mathcal{N}^{-}(v)=\mathopen{}\mathclose{{}\left\\{u\mid(v,u)\in\mathcal{E}}\right\\}$.
Moreover, we will call ‘$k$-step neighborhood’ the set of all vertices whose
(unoriented) distance from $v$ is less than or equal to $k$. Furthermore, we
can have multiple types of vertices and edges specified in finite sets
$T_{\mathcal{V}}$ and $T_{\mathcal{E}}$, respectively. Mapping
$\operatorname{type}(v)\in T_{\mathcal{V}}$ attributes the type to every
vertex $v\in\mathcal{V}$ and, analogically, mapping $\operatorname{type}(e)\in
T_{\mathcal{E}}$ returns the type of each edge $e\in\mathcal{E}$. Given a
vertex $v$ and the graph it comes from, the task is to classify the vertex
into one of the classes from set $\mathcal{C}$. This task is the most general
one—for instance, we may be interested in classifying only vertices of a
certain type $t_{v}\in T_{\mathcal{V}}$, or given labels of a subset of
vertices in the graph inferring the labels of the rest, or a combination of
both. Examples are provided later. As opposed to Euclidean domains, graph
domains pose new challenges in data transformation. For example, generalizing
a convolution operation for graphs is not straightforward. In an image, the
neighborhood of each pixel contains a constant number of surrounding pixels,
which can be ordered by their relative position to the given pixel, whereas in
graphs, number of neighboring vertices of a vertex can be an arbitrary number,
and neighbors are not ordered. Therefore, formulas for message computation in
a sound inference procedure should accept any number of arguments and should
be invariant to the permutation of the arguments. The key idea of HMill-based
message computation is that neighborhood $\mathcal{N}(v)$ of a node
$v\in\mathcal{V}$ can be considered a bag and each of vertices neighboring $v$
can be regarded as an instance of this bag. A message predictor can then be
implemented using MIL paradigm, which guarantees that all requirements for a
message computation are fulfilled since a bag in a MIL problem is defined to
be an unordered, arbitrarily large set of instances. Thus, all existing
knowledge about MIL models can be applied. Formally, for a node
$v\in\mathcal{V}$, we construct an HMill bag data node
$\operatorname{bn}(\mathcal{N}(v))$ and define a message predictor using the
framework. As we shall demonstrate in the next paragraphs, using the HMill
framework for graph inference should be regarded not as a specific method, but
rather a modelling toolkit for diverse types of input graphs, which due to its
flexibility opens a new realm of modelling choices.

####  6.2.1 Unrolling the procedure

The standard message-passing procedure in GNN s goes as follows: 1. Hidden
states $\mathbf{h}_{v}^{0}$ are initialized for each vertex $v$ in the graph
2. In each iteration $t$, hidden states $\mathbf{h}_{v}^{t}$ of each vertex
$v$ are updated using a function taking hidden states of vertices in
$\mathcal{N}(v)$ from iteration $t-1$ as input. 3. After performing $k$
iterations, hidden states $\mathbf{h}_{v}^{k}$ are processed in a ‘readout’
phase. This usually involves classification of vertices based on their states,
or somehow aggregating hidden states of all vertices in the graph to produce a
representation of the whole graph. Each iteration is computed synchronously in
the sense that hidden states of all vertices are updated at once. The main
advantage of this approach is that intermediate results can be reused and,
therefore, no redundant computation is performed. Rather than this graph-
centric point of view, we advocate a different, node-centric perspective in
this work. We assume that a single vertex is one sample (observation) in
machine learning terms, and the graph is only queried for information
characterizing the vertex. Let us emphasize that ‘unrolling’ a computation on
a single vertex in this way does not contradict the dependence claim from the
beginning of this chapter. We will still be able to infer knowledge (for
instance a class label) about all vertices in the graph taking into account
their interactions with other vertices. The aforementioned perspective only
demonstrates that inference on different vertices in the graph can be
considered a completely independent procedure and computed as such at the
price of limited reuse of intermediate results. See Figure 6.2 for a specific
example, where for simplicity, we show only one step of message-passing
inference. On the left, there is an illustration of how one step of message-
passing inference is carried out in an undirected graph. Messages, which we
universally denote by $\mu$, are sent in both directions across each edge.
This way, we can compute a new hidden state in each vertex in the graph by
aggregating information from $1$-step neighborhood. Equivalently, the new
hidden states can be computed for each vertex ‘independently’ from smaller
subgraphs of the original graph, which is illustrated for three specific
vertices on the right. Figure 6.2: A depiction of how a message-passing-based
inference can be unrolled. Arrows $\mu$ signify the flow of information
(messages). Three random vertices $u$, $v$ and $w$ from the graph on the left
are highlighted, and the $1$-step neighborhood of vertex $v$ is shown as well.
The equivalent and independent procedures for the hidden state computation are
illustrated on the right. See text for details. Assume a vertex $v$, which we
will also refer to as the ‘current vertex’, for which we want to run the
inference procedure using the graph as a source of information. This
information may include features defined in the vertex itself or in its
incident edges, or similar knowledge derived from other vertices in the graph,
which are locally connected to the node, directly with an edge or longer path.
In this manner, we can consider a ‘hidden state’ $\mathbf{h}_{u}^{t}$ of a
node $u$ located close to $v$ in the graph, as a ‘feature’. The fact that
these types of features derived from the vertex or its neighborhood are
inherently heterogeneous, differently structured, and some of them recursively
computed, directly motivates the use of HMill for their processing. In the
following text, we explain how this can be done, and two more specific use
cases are presented afterwards. First, assume that $v$ and all vertices in its
neighborhood $\mathcal{N}(v)$ are described as fixed size vectors
$\mathbf{h}$. We construct a bag data node
$\operatorname{bn}(\mathopen{}\mathclose{{}\left\\{\operatorname{an}(\mathbf{h}_{u})\mid
u\in\mathcal{N}(v)}\right\\})$ containing as many instances as is the degree
of $v$. Due to the properties of bags in MIL paradigm, we can obtain such bag
data node for every node in the graph and all of them will follow the same
schema. Thus, we can define a single model for processing all of them. In the
HMill framework, bag model nodes consist of an instance model $f_{I}$, an
aggregation function $g$ and a bag model $f_{B}$. In this case, $f_{I}$ would
be responsible for the initial transformation of vector representation of
neighboring vertices, $g$ would output a single vector, and $f_{B}$ would
perform the final transformation. The first modelling choice is how to include
representation of vertex $v$ itself to the computation. We elaborate on this
later. The described HMill-based inference procedure using $v$ and its
$1$-step neighborhood $\mathcal{N}(v)$ is equivalent to performing one step in
the message-passing scheme. Recall that an instance model $f_{I}$ can be any
valid HMill model. By a recursive argument, if we define $f_{I}$ to be itself
a bag model node with the same architecture as in the $1$-step case, the
overall model performs two-step message-passing procedure. Since every HMill
model returns a vector, $f_{I}$ would also return a vector for each neighbor
(instance) and the rest of the computation in $g$ and $f_{B}$ remains the
same. Therefore, with HMill we can implement a message-passing procedure with
an arbitrary number of steps. The depth of the model and resulting tree sample
representations will be $k+1$, where $k$ is the number of steps.

####  6.2.2 Modelling vertices

Vertex $v$ itself and other nodes in the graph may be described by a vector of
features or raw data encoded in different data types, or a combination of
both. One of the main advantages of using HMill-based inference is that a node
description can be provided as an arbitrarily complex structure, for instance,
a JSON document. The inference procedure with a complex description of
vertices can be split broadly into two abstract parts—in the first one,
important knowledge is distilled from each node, and in the second part, it is
distributed across the graph and combined with data in other vertices. We have
already demonstrated how hierarchical structures are modelled with HMill in
Chapter 5. Thanks to the fact that message-passing inference itself can be
implemented with HMill models as well, the whole computation seamlessly
integrates into the framework. With the reasonable assumption that node
representations all match the same schema, the whole inference procedure for
one node can be implemented as a single HMill model, which can be trained in
an end-to-end manner. Encoding of the $\operatorname{type}(v)$ into the node’s
representation is also possible either trivially as a one-hot encoded vector,
or as some more complex encoding. We have shown that information from the
neighborhood of vertex $v$ can be aggregated with a bag model node to compute
a message, however, how is a representation of the vertex itself added to
computation is yet to be specified. This is done with an HMill product node,
whose children are one bag node describing the neighborhood and one other data
node for vertex representation, whose type depends on the data type of $v$’s
description—array node if $v$ is described by only a vector of features, or
more complex data nodes in other cases. In other perspective, knowledge stored
in the neighborhood of $v$ can be regarded as just another deep subtree in
product node representation of $v$. In JSON parlance, this is equivalent to
adding another key-value pair into the JSON object that describes $v$. Refer
to Figure 6.3 for a depiction of this process. Node $v$, together with its
$1$-step neighborhood of six vertices, is drawn on the far left and right
sides. In both cases, the neighborhood of $v$ is modelled by a bag node with
six instances in this particular case, each of which is described further by
an arbitrarily complex HMill sample tree, depicted by a blue subtree. This
subtree can describe not only the neighboring node itself (for instance its
features) but its neighborhood as well. In the left picture, $v$ is described
by a feature vector, whereas in the right one, we used a red subtree
signifying that in fact any valid HMill (sub)tree can be used to describe $v$.
Finally, direct feedback is a phenomenon that occurs when an output of a
vertex is influenced by the vertex itself due to information being propagated
in several iterations. For instance, consider an edge
$\mathopen{}\mathclose{{}\left\\{v,u}\right\\}\in\mathcal{E}$. If we want to
infer a label for $v$, according to the inference procedure we first evaluate
results of $f_{I}$ on all vertices in $\mathcal{N}(v)$ and then aggregate
them. However, if more steps are performed and therefore an intermediate
result on vertex $u$ is computed from its neighborhood as well, it holds that
$v\in\mathcal{N}(u)$, and output from vertex $v$ will be influenced twice by
its representation. Depending on the domain, this may be desirable or not. In
HMill-based inference, this can be easily remedied by redefining a bag
representing $u$ to contain vertices $\mathcal{N}(u)\setminus v$.

(a) feature vector

(b) HMill subtree

Figure 6.3: An illustration of how information about current node $v$, which
receives the message, is added to its representation using a product data
node. The information (features) describing $v$ can be in the form of a
feature vector (6.3(a)) or a more complex HMill representation (6.3(b)). The
fact that a vertex in the graph is described by an HMill compatible structure
is illustrated with an arrow . The types of arrowheads used for different
purposes is consistent with the rest of the thesis—for graph edges, we use ,
for illustrating parent-children relation in sample and model trees  or , and
are used as generic pointers. See text for details.

####  6.2.3 Modelling edges

In some applications, most useful information is stored not in vertices
themselves, but in relationships between them signified by edges. Analogically
to vertices, each edge can be represented by a more or less complex HMill
structure. Instead of regarding each neighbor $u\in\mathcal{N}(v)$ as an
instance, we can consider the relationship
$\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$ an instance and include
information about the relationship
$\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$ instead of about $u$ only. If
both $u$ and $\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$ are important,
bag instances can be defined as a tuple
$\mathopen{}\mathclose{{}\left(u,\mathopen{}\mathclose{{}\left\\{u,v}\right\\}}\right)$.
The most straightforward way to represent such instance is by a product data
node with two subtrees, one for $u$ and one for
$\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$. Graph edges of different
types, or directed edges, which are just a special case of the former, or
both, can be modelled in different ways. The more naïve option is to use one-
hot encoding on the type of an edge and add this to the representation of
$\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$. Another way is to model each
type of edges separately. This corresponds to constructing $\lvert
T_{\mathcal{E}}\rvert$ bags ($2\lvert T_{\mathcal{E}}\rvert$ if graph edges
are both typed and directed), one for each type of edges. For each edge type
$\tau\in T_{\mathcal{E}}$, we will include only neighbors connected with the
edge of type $\tau$ into the bag. In other words, for each $\tau\in
T_{\mathcal{E}}$, we define one bag
$b_{\tau}=\mathopen{}\mathclose{{}\left\\{u\mid
u\in\mathcal{N}(v),\operatorname{type}(\mathopen{}\mathclose{{}\left\\{u,v}\right\\})=\tau}\right\\}$.
Which way is more suitable again depends on application and specific examples
are provided later. In Figure 6.4, we sketch the latter option. The situation
is similar to one depicted in Figures 6.3(a) and 6.3(b), however, in this case
the graph being modelled is oriented. This means that we deal with two types
of edges with a possibly different interpretation. Thus, we use two different
bag nodes, one for incoming and one for outcoming edges. Recall, this implies
that each incoming edge will be modelled with the same (sub)model tree and the
same applies to outcoming edges. However, the parameters will not be shared
between incoming and outcoming edges. In the picture, we use golden and blue
subtrees to illustrate a representation of incoming and outcoming relations in
a general way, nevertheless, both representations can be defined as equal.
Note, that these subtrees can be used to describe not only the neighboring
node $u$ but also the relationship $\mathopen{}\mathclose{{}\left(v,u}\right)$
or $\mathopen{}\mathclose{{}\left(u,v}\right)$, which is depicted
schematically in the figure. Figure 6.4: One possible way of dealing with
multiple edge types. In this particular case, edge types correspond to
incoming and outcoming edges with respect to current node $v$. See text for
details.

####  6.2.4 Training the model

Refer to Figure 6.5 for a specific example of how message-passing inference
can be translated into HMill. For simplicity, we assume that no information
about edges themselves is known, however, each node is described by a vector
of features. The figure shows how HMill sample for implementing two steps of
message-passing inference looks like. First, each of the vertices $a$, $b$,
and $c$ is described by a product node merging feature vector representing the
vertex and information obtained from its neighborhood, which is stored in a
bag node. Then, we store all three product data nodes into another bag data
node, which is merged with feature vector of $v$. Note that all bag node
definitions are valid, as the product data node representation of vertices in
$\mathcal{N}(v)$ follow the same HMill schema.

(a) Original subgraph

(b) Node representations

(c) Resulting sample

Figure 6.5: HMill sample tree extracted from the neighborhood of vertex $v$
for performing $2$-step message-passing inference in a graph, where vertices
are described by vectors of features. The original subgraph induced by
vertices in the $2$-step neighborhood of $v$ is depicted in (6.5(a)). The
picture (6.5(b)) in the middle shows the same subgraph, where feature vectors
of dimension $4$ describing vertices are drawn, together with arrows
illustrating, how information is transferred in each of two steps. The
uppermost horizontal edge is not present in (6.5(b)), since more than $2$
inferential steps would be required to reach it. The final sample tree, drawn
aligned with the original subgraph to facilitate comparison, is in the third
figure (6.5(c)). For bag node representations of neighborhoods of vertices
$a$, $b$, and $c$, we used the batch version to save space (see Section 4.5).
See text for details. As the whole model is a single HMill model, it can be
trained in an end-to-end fashion. Thus, processing of the data in vertices and
edges is optimized jointly with a message-passing part. If the node-centric
perspective on message-passing inference is used, the standard training
procedure can be employed. One or more graphs are collected, and in each
iteration of the training loop, several vertices are sampled to form a
minibatch, the unknown quantity (for instance a label) is inferred, and errors
are backpropagated through the model according to Algorithm 1. Due to the fact
that HMill offers a principled way of dealing with missing data at different
levels of abstraction (Section 4.3), the only requirement on the data is that
HMill subtrees describing vertices (or edges) must follow the same schema. The
framework elegantly deals with other possible discrepancies in the data, such
as missing information on vertices or edges or empty neighborhoods. One
concern about working with any method for learning from graphs is limited
scalability. The number of vertices in a $k$-step neighborhood of a vertex
grows exponentially, and with increasing $k$ gradient computation quickly
becomes infeasible. This is still an open problem [Wu et al., 2019], however,
it can be partially remedied by sampling. One possible solution is to employ a
variant of importance sampling, which in HMill is easily done by replacing bag
data nodes representing neighborhoods with their weighted variants (Section
4.1). For instance, in problems from the computer security domain, positive
signs of malicious activity can be sometimes observed sparsely, however, they
convey the most useful information. In this case, the sampling procedure can
be designed to be biased towards malicious observations. The specific examples
are provided later in this chapter and in the next one. In our experience, the
performance gain achieved by increasing the number of steps $k$ in the
inference procedure quickly saturates. A pleasant side-effect of sampling from
bags (neighborhoods) is an introduction of stochasticity into learning, which
is known to be beneficial for optimization with the gradient descent class of
algorithms.

####  6.2.5 Differences to GNNs

We have shown how the HMill framework can be extended to perform message-
passing inference in graphs, using bag nodes to model neighborhoods and
product nodes to merge heterogeneous information. The resulting method retains
all favorable properties of the message-passing computation scheme while
offering several additional benefits. Due to the fact that universal function
approximators are used instead of analytical closed-form expressions for
message computation, the space of all inferable hypotheses becomes larger. At
the same time, the method does not require specifying a family of functions
for message computation, avoiding potential approximation error stemming from
selecting a wrong family. The HMill-based inference, as well as GNN s, assumes
that a subgraph induced by $k$-step neighborhood contains all relevant
information and that the further a node or edge is located in the graph, the
less important it is. The proposed inference method starts at the most distant
vertices and proceeds towards $v$, in each step using bag and product model
nodes to aggregate and merge distilled knowledge relevant for inferring the
class label or other unknown quantity. Upon closed inspection, it turns out
that in comparison to how GNN s operate, our method uses a slightly different
computational model, parameters are shared at different locations, and
different overall philosophy to the graph inference is adopted. Whereas GNN s
drew most inspiration from recurrent networks with each layer estimating graph
convolution operation, which leads to powerful, but hard-to-scale models, the
HMill-based approach regards the graph merely as a knowledge database for
querying information about vertices and their relationships. Hence, both
approaches are suitable for different circumstances. As the main goal of the
framework is to unlock new possibilities in the processing of raw data, the
scalability and modelling flexibility are two of the main concerns.
Specifically, in this thesis, we provide examples of use cases from the
computer security domain, where graphs are huge and node degrees come from a
wide range of values. In our opinion, bag nodes are the perfect candidate for
modelling neighborhoods, which in most graphs observed in real-world scenarios
follow a power law distribution, and described bag sampling schemes make
inference tractable. Moreover, the alternative node-centric view to inference
enables straightforward parallelization as inferring unknown values for
vertices can be considered independent. Also, because a single value can be
inferred at the time, sampling techniques at the level of batches are
possible. For instance, if we want to perform binary classification of
vertices in graphs, where the vast majority of vertices are positive, it is
still possible to construct balanced minibatches for training. Another benefit
is that a problem of missing data at different levels is solved in the HMill
framework in a principled way and any heterogeneity in the graph, such as
different vertex and edge types can easily be modelled with product nodes.
Finally, thanks to the modelling flexibility of HMill in general, complex
descriptions of vertices and edges can be integrated and the resulting model
is still one HMill model, which means that the inference procedure together
with the processing of representations is optimized jointly. We believe that
due to little requirements imposed on the input graphs in terms of size,
degree distributions and structure, the proposed approach is versatile enough
for utilization in many potential applications.

###  6.3 Behavior-based malware classification in graphs (use case)

An integral component in any antivirus software is a capability to recognize
and protect users from malicious programs. Nowadays, harmful software takes
many forms. Trojans allow unauthorized access to computers, worms spread
themselves over many machines, ransomware programs encrypt the whole hard
drive of the user and demand payment in exchange for decryption secret,
spyware records sensitive information like passwords, and finally, adware
programs generate revenue by displaying an advertisement on the screen. These
examples are just a few of many types of existing malware. Current
technologies designed for defense can be generally categorized into two
categories. Traditional signature-based mechanisms rely on extracting a
predefined ‘signature’ from the inspected file and comparing it to a known
database of threats. A signature could be based for instance on MD5 or SHA1
hash of the whole executable or its part in the most trivial case. The
advantages of the signature-based approach is low false positive rate, and,
due to its simplicity, high detection speed, crucial characteristics of any
antivirus software. On the other hand, once a new threat appears, anti-malware
solution providers have to generate its signature and distribute it to
clients. Oftentimes, this leads to blacklists becoming outdated and users
being vulnerable to the evolved malware for a period of time. The Cisco 2017
Annual CyberSecurity Report [Cisco Systems, 2017] states that 95% of all
analyzed malware files were not even a day old, which indicates that obsolete
signature databases are more the rule than the exception. Moreover, a great
deal of resources is required to keep blacklist databases updated all the
time. More sophisticated malware makes use of encrypted or polymorphic code
segments, which makes signature generation very hard, since many executables
with the same behavior but different signatures can be generated.

These issues can be addressed by designing a rule-based system or a machine-
learning system that takes into account the structure of the file or the
behavior of the program. In static analysis, dangerous attributes are sought
after within the structure of the executable, its code or its other
properties, whereas in dynamic analysis, the executable is run and observed in
a sandboxed environment. For this reason, this complementary approach to
signature-based is usually referred to as behavior-based. Once suspicious
behavior is detected, the antivirus may remove some important permissions of
the program or stop the execution altogether. Different aspects of the
program’s behavior can be observed, for example, changing register values,
installing rootkits, accessing specific files or services, locking the
executable from removal, disabling security controls, registering the
executable for autostart, attempting to detect and evade sandbox, and so
forth. Behavior-based techniques effectively tackle all of the signature-based
approach flaws, as malicious behavior is the defining feature of all malware
and is therefore difficult to conceal. Two programs with different signatures
but identical behavior are treated in the same way as well. Furthermore, this
approach detects new malware families, which might be programmed in a
completely different way, but exhibit the same behavior as known malware. The
downside of behavior-based methods is that they are usually slower. In some
cases, malware programs can resist detection and removal by spawning two
processes that restart each other, should one of them get killed. In other
cases, malware activities are exhibited not by a single process, but a group
of processes that may not all be instances of the same executable. Malware is
usually first installed with a dropper program, the malicious action is
performed by another program, and both can use additional modules. For
instance, Windows hooks play a major role in most of keyboard logger programs,
and droppers that install this malware have an extra dynamically linked
library or dll file [Williamson and Gorelik, 2014b]. As a result, malware
activities can be hard to detect when all the components are examined
independently. This motivates the use of a graph structure to encode the
overall behavior in the system for the detection.

####  6.3.1 Data description

For the experiment presented in this chapter, we once again partnered with
Avast researchers to try HMill-based graph inference on real-world data from
clients. Specifically, we focus on the Windows operating system, which is
installed on a bulk of clients. As described in patents [Williamson and
Gorelik, 2014b, Williamson and Gorelik, 2014a], the Avast client keeps a
regularly updated snapshot of the operating system for the detection purposes,
which is encoded as a graphical structure called Identity Protection (IDP)
graph. The graph consists of two types of vertices representing executable
files and processes, therefore
$T_{\mathcal{V}}=[\texttt{EXECUTABLE},\texttt{PROCESS}]$. Executable vertices
correspond to executable files in the system that have been run, and are
identified by the full file path. Therefore, multiple copies of the same file
will result in multiple vertices in the graph. Process vertices represent all
processes which either run in the system or spawned another process that is
still running. Executable files and their processes are represented
separately, because different processes may behave differently, for example
when foreign code is injected into a running process. Multiple types of graph
edges encode events that occurred in the system measured from the Windows
kernel and other sensors. Edge types (relationships) are summarized in Table
6.1. Note that the table contains only the most relevant relationships
published in the aforementioned patents, which is a proper subset of all $24$
types of edges. Table 6.1: Edge types (relationships) $T_{\mathcal{E}}$ used
in IDP graphs that were published in [Williamson and Gorelik, 2014b]. name
participants interpretation INSTALLS $p$ $e$ Process $p$ wrote an executable
$e$ to disk. SPAWNS $p_{1}$ $p_{2}$ Process $p_{1}$ started another process
$p_{2}$. INSTANCE OF $e$ $p$ Process $p$ is an instance of executable $e$.
REGISTERS $p$ $e$ Process $p$ stored some configuration data about executable
$e$, for instance a registry key causing reboot survival. CODE INJECT $p_{1}$
$p_{2}$ Process $p_{1}$ caused process $p_{2}$ to run foreign code. KILLS
$p_{1}$ $p_{2}$ Process $p_{1}$ attempted to kill another process $p_{2}$.
Moreover, each vertex in the graph is described by an XML document storing
additional multimodal information. So far, we have described only how JSON
documents are processed with the HMill framework. Even though XML documents
are processed analogically, in this case, documents can be directly converted
to JSON s, thus one can think of this as a preprocessing step and then apply
exactly the same procedure defined previously for JSON documents. An HMill
schema of the vertex description inferred from the training data can be found
in Figure C.6 in Appendix C. It includes a unique "id", a fully qualified
"name", a "size" of the file in bytes, a "sha256" hash, information about
"connections" made, "named_objects" used, or registry operations ("regops")
performed for executables, and finally, a set of ‘characteristics’.
Characteristics are binary predicates measured by means of static and dynamic
analysis of the given executable or process designed to discriminate between
malicious and clean software. For this experiment, we used a set of $164$ such
characteristics. Refer to Table C.1 in Appendix C for specific examples. Some
information is specified only for one type of vertices, for instance, it does
not make sense to define a size of a process node. The same holds for
characteristics, therefore, ternary logic with values
$\mathopen{}\mathclose{{}\left\\{\texttt{TRUE},\texttt{FALSE},\texttt{UNKNOWN}}\right\\}$
is used for their specification. In Appendix C, we provide Figure C.1
sketching how the shape of one weakly connected component of an IDP graph
looks like. Also, in Figures C.2 and C.4 we show specific instances of
subgraphs, the first one coming from an operating system without infection and
the other one from an infected system. For even more details, please refer to
the patents cited above.

####  6.3.2 Experimental settings

For experimental purposes, we obtained a labeled, anonymized dataset
consisting of around $4\cdot 10^{5}$ anonymized graphs collected from Avast
clients during February and March 2020. The goal was to classify a subset of
approximately $7,9\cdot 10^{6}$ executable nodes present in graphs into a
malicious or benign class based on only information encoded in the graph. We
did not attempt to classify all executable vertices in graphs since some of
them are signed by trusted companies and are therefore considered clean and
ignored by the antivirus client. In this experiment, we demonstrate that the
main strength of the HMill framework lies in its modelling flexibility. We
defined three types of increasingly more complex models operating with
different amounts of information, all depicted in Figure 6.6. Assume we want
to classify an executable vertex $v$. The first model (Figure 6.6(a)) $m_{1}$
processes only characteristics of $v$ encoded into a vector of length $2\cdot
164$ with one-hot encoding444We used two ‘bits’ per each characteristic. We
set the first bit if the characteristic is TRUE, the second one if it is
FALSE, and neither if it is UNKNOWN.. It is defined as a standard feedforward
neural network with two hidden layers of $300$ neurons with
$\operatorname{ReLu}$ activation function and $2$ output neurons, one for each
class. A similar model is deployed in Avast pipeline as well. In HMill
terminology, $m_{1}$ is a single array model node $\operatorname{am}(f)$,
where $f$ is the mapping implemented by the network. The second model $m_{2}$
(Figure 6.6(b)), defined as an HMill product model node, uses all information
available on $v$, which is summarized in the schema in Figure C.6. Thus, it
contains a subtree of depth one for modelling characteristics as well, with
the same structure as $m_{1}$.555In fact they differ in a minor detail that
$300$ neurons were used in the last layer instead of $2$ as in the case of
$m_{1}$, since the output it is used further in the model. In the figure, this
is indicated with a golden frame. Other than that, a product model node
$m_{2}$ contains a child for every key of the aforementioned schema, with the
exception of id, node_type666Modelling node types does not make sense in this
context since all modelled vertices are executables. and sha256. For sizes in
the size field, we used mapping $x\mapsto\log(1+x)$. The procedure for
constructing an HMill sample tree for $m_{2}$ from an XML document is the same
as used in the device identification problem from Section 5.2.2—we used
histograms of tri-grams for modelling string values, and one-hot encoding for
the representation of categorical variables. An HMill model tree was built so
that every neural network component in model nodes has the output dimension of
$50$ if not specified otherwise, in exactly the same way as in Section 5.2.2.
Last model $m_{3}$ implements a $1$-step message-passing inference as
described in this chapter, therefore, apart from data about $v$, it also
utilizes information stored in $1$-step neighborhood $\mathcal{N}(v)$ of the
vertex in the graph. This was implemented exactly as described in Section
6.2.3. The children (subtrees) of the topmost product model node of $m_{3}$
consist of three sets—one for modelling data stored in the vertex itself,
which is the same architecture as in $m_{2}$ (marked with a blue box in Figure
6.6(c)), and two bag model nodes for modelling incoming and outcoming edges.
For each sample, every edge (instance in one of the corresponding bag data
nodes) is represented as the one-hot encoded edge type and the description of
the neighbor, merged together in a product data node. The neighbor
description, in the figure depicted as a purple subtree, was implemented using
the same architecture as for the description of $v$, note however that model
parameters are not shared. Consequently, model $m_{3}$ can infer the label not
only from the information known about the executable vertex itself but also
from how the file (or its process instances) relates to other entities in the
whole system. Since both process and executable vertices may be present in
$\mathcal{N}(v)$, we include a one-hot encoded vertex type into all neighbor
representations. Again, all neural network component architectures were
defined using the simple rule described above for the case of $m_{2}$. We used
$\operatorname{ReLu}$ non-linearity everywhere, except of last layers of
instance models in bag model nodes, where we used the $\tanh$ activation
function since in our experience it performs better with subsequent
aggregations. The very last layer in the root of every model tree was followed
by a $\operatorname{softmax}$ layer to obtain class probabilities.

(a) $m_{1}$

(b) $m_{2}$

(c) $m_{3}$

Figure 6.6: Three different models used in experiments. Golden and blue
rectangles indicate that the same architecture (schema) was used for the
model. However, inner neural network components may have different numbers of
neurons in the last layer. For example, model $m_{1}$ (6.6(a)) has $2$ neurons
in the last layer of $f$, one for each class. (Sub)model in the golden subtree
in $m_{2}$ does not output $2$ values, as the output is transformed further in
the model.

In some cases, vertex $v$ has a high number of neighbors. This happens, for
example, when a ransomware program spawns a lot of processes to perform data
encryption. As most of the time all of these processes have the same
attributes, it is enough to inspect a small subset of them. Therefore, to save
time, we employed an importance sampling scheme and instance weighting
described in Section 4.1. We include into both bag data nodes describing
incoming and outcoming edges of vertex $v$ at maximum $K=100$ neighbors
connected with the edge of the same type. This way, all sparsely represented
relationships are placed in the bag every time, however neighbors connected to
$v$ with a common edge type are sampled. A $15\%$ proportion of all $4\cdot
10^{5}$ graphs in the dataset was left out for testing. We counted $8\cdot
10^{5}$ positive vertices out of $6\cdot 10^{6}$ vertices in the training
graphs and $1,5\cdot 10^{5}$ positive vertices out of $1.2\cdot 10^{6}$ in the
graphs held out for testing. All models were trained in the exact same way,
performing one epoch (iteration over the whole training set) using minibatches
of size $100$ and Adam optimizer [Kingma and Ba, 2014] with the default
parameters. Weights in each neural network layer were initialized with the
Glorot normal initialization [Glorot and Bengio, 2010], and all biases were
initialized to zero. Loss function was weighted binary cross entropy:
$\mathcal{L}(\theta)=-\frac{1}{n}\sum_{i=1}^{n}w_{1}y_{i}\log
p(1|v_{i};\theta)+w_{0}(1-y_{i})\log p(0|v_{i};\theta)$ (6.6) where $n$ is the
number of classified vertices in the batch, $v_{i}$ HMill representation of
$i$-th vertex in the batch,
$y_{i}\in\mathopen{}\mathclose{{}\left\\{0,1}\right\\}$ its binary label, and
$p(1|v_{i};\theta)$ conditional probability of $i$-th vertex being malicious
according to the model parametrized by $\theta$. Note, that
$p(0|v_{i};\theta)=1-p(1|v_{i};\theta)$ for each $i$. Weights $w_{0}$ and
$w_{1}$ are real numbers, specifically, we used $w_{1}=0.1$ and $w_{0}=0.9$ to
express that we prefer false-negative to false-positive errors.

####  6.3.3 Experimental results

Figure 6.7: Description of the performance of models $m_{1}$, $m_{2}$ and
$m_{3}$ with a Precision-Recall (PR) curve. A canonical way to compare the
performance of several binary classifiers in the presence of imbalanced
classes is to plot a PR curve or a Receiver Operating Characteristic (ROC)
curve. The standard scalar metrics computed from a confusion matrix require a
specific threshold value for their computation. PR and ROC curves enable us to
study values of such scalar metrics for multiple threshold values
simultaneously. PR curve captures a tradeoff between recall and precision, two
inversely related quantities. ROC curve describes how true positive rate
(recall) and false positive rate change for different values of the threshold.
To express the performance with a single number taking into account multiple
thresholds, we can measure the area under such curves. This gives rise to Area
under Precision-Recall curve (AUPRC) and Area under Receiver Operating
Characteristic curve (AUROC) metrics, which both range from zero to one. Even
though this is a convenient way to quantify the performance, the values of
AUPRC and AUROC may sometimes be misleading. Depending on the circumstances,
different parts of the curves may be of different significance. For instance,
in the cybersecurity domain, false positives are much less desirable than
false negatives. Usually, high enough threshold values are selected to keep
false positive rate sufficiently low, and precision sufficiently high. This
corresponds to rather small regions in both curve types, however, AUPRC and
AUROC are evaluated using the whole curve. For the same reason, we include an
ROC curve with a logarithmic $x$-axis (false positive rate). All curves drawn
in this work are made of one hundred points with $x$-values sampled uniformly
from the linear (or logarithmic) scale and $y$-values obtained with linear
interpolation. We trained all three aforementioned types of model $m_{1}$,
$m_{2}$, and $m_{3}$, and evaluated their performance. The resulting PR curve
is drawn in Figure 6.7, and ROC curve with linear and logarithmic scale of
$x$-axis is in Figure 6.8. Perhaps not surprisingly, the more information the
model processes, the better results it achieves. By inspecting a $1$-step
neighborhood of the vertex corresponding to an executable implemented in model
$m_{3}$, we were able to attain high true positive rate while keeping false
positive rate low. We could not achieve this with model $m_{2}$, which
processes only an XML document description of the vertex. This makes model
$m_{3}$ much more useful in practice than $m_{2}$, despite the little
difference in the AUROC metric. This experiment confirms that simply using
more data leads in most cases to better results. All three models $m_{1}$,
$m_{2}$ and $m_{3}$ were designed as ‘baseline’ and architectures and training
procedure can surely be improved. However, unutilized data sources are often
the low-hanging fruit. Thanks to HMill, modelling of such sources is
straightforward.

(a)

(b)

Figure 6.8: ROC curve corresponding to models $m_{1}$, $m_{2}$ and $m_{3}$
used in the experiment. We provide both variants with linear (6.8(a)) and
logarithmic (6.8(b)) scale of $x$-axis.

##  Chapter 7 Modelling interactions in heterogeneous networks

In this next-to-last chapter, we present one last and the most specific use
case, where the HMill framework comes in handy. In this problem, very little
is known about objects of interest themselves. Instead, interactions between
objects are provided in several bipartite graphs specifying binary relations
between objects. Thus, we move from a purely ‘feature-based’ approach
presented in Chapter 5 and also in Chapter 6 to learning merely from relations
between objects. The two approaches, each with its own advantages and flaws,
are to some extent orthogonal and can be combined as we discuss later. To
empirically test the proposed method, we investigate the task concerning an
extension of a blacklist of malicious computer domains based on their
interactions with other entities in the network. In the last part of the
chapter, we present the results of the experiments on the data provided by
researchers from Cisco Cognitive Intelligence. In the domain of network
security analysis, the focus is shifted from the level of individual client
machines to whole networks of clients and their protection by detecting and
neutralizing malicious activities happening in the network. Being able to
decide for any entity (IP address, domain, email) in the network whether it is
safe for communication or not is a crucial requirement for implementing
passive security measures in the first line of defense, for instance, simply
blacklisting known threats. However, keeping such blacklist updated is a
difficult task for the same reasons we listed earlier when discussing malware
signature databases—due to the quick growth of the number of existing threats,
it is nearly impossible to keep lists up to date manually. Therefore, methods
for the automatic extension of blacklists with uncovered threats are
researched.

###  7.1 Prior art

Prior art on the problem of blacklist extension can be broadly split into two
fundamental approaches. In classifier-based approach [Zhang et al., 2008,
Antonakakis et al., 2010, Bilge et al., 2011], blacklists are used to obtain
labels for training a classifier for determining the maliciousness of unseen
objects in the future. One of the most constraining assumptions to achieve
good performance is having high-quality training data. It is known that public
blacklists are oftentimes incomplete due to the delay in their updates or the
fact that some existing threats have not been discovered yet. If all entities
not present in the blacklist are labeled positive in the training phase, the
classifier will be inevitably shown samples with an incorrect label. The
second approach is graph-based [Coskun et al., 2010, Philips et al., 2012,
Carter et al., 2014, Oprea et al., 2015, Rezvani et al., 2015, Ahuja et al.,
2019] which assumes that malicious activities are localized in a graph,
forming communities with sharp boundaries [Collins et al., 2007, Yu et al.,
2010]. Methods following this approach build a graph representation of
relations between objects such that graph analysis algorithms can be utilized
to extend knowledge from the blacklist. High malicious score is first
initialized in vertices corresponding to objects found in the blacklist, which
are referred to as seeds or tips, and afterwards propagated through the graph
using formulas involving the maliciousness of the neighbors and weights
defined on edges. After some number of steps, vertices with highest malicious
scores that are not present in the blacklist are returned as candidates for
blacklist extension. As a result, no data apart from relations themselves are
collected. This is useful in situations when data is inaccessible or expensive
to obtain. On the other hand, to the best of our knowledge, algorithms for
propagating the maliciousness through graphs use fixed formulas for
propagation, cannot use complicated representations of relations, and have few
to none tunable parameters. For instance, the most popular Probabilistic
threat propagation (PTP) algorithm [Carter et al., 2013, Carter et al., 2014],
which we further elaborate on below, propagates malicious scores in the form
of scalar values throughout a graph with scalar weights on edges. Definition
of edge weights is the only parameter of the method together with a number of
steps performed. The small number of parameters limits the use of methods for
different problems, as the algorithms are oftentimes tailored to perform a
particular task on a specific dataset and may not be robust to changes in the
application domain. The graph-based approach to blacklist extension is similar
to the inference problem in probabilistic graphical models we described in the
previous chapter and in fact some methods are built on the same ideas [Carter
et al., 2013, Carter et al., 2014, Manadhata et al., 2014]. However, the
possibility of using GNN s for the task has not been studied for two main
reasons—their high computational complexity coupled with characteristically
large size of graphs observed in the domain, and a lack of public datasets.
Therefore, the cybersecurity domain requires different techniques, despite
many similarities in the task at hand. We end this section by reviewing a PTP
algorithm [Carter et al., 2013, Carter et al., 2014], which is a canonical
graph-based algorithm used in cybersecurity. Given an undirected graph, it
estimates the probability $P(v_{i})$ of a vertex $v_{i}$ being malicious from
its connection to other vertices and from a set of known malicious vertices,
which are assumed to come from the blacklist, using the following equation:
$P(v_{i})=\sum\limits_{v_{j}\in\mathcal{N}(v_{i})}w_{ij}P(v_{j}|v_{i}=0),$
(7.1) where $w_{ij}$ are positive edge weights assumed to be normalized to
one, $\sum_{j}w_{ij}=1$. Conditioning the threat to $v_{i}=0$ avoids unwanted
direct feedback of the vertex to itself. Since exact solving of (7.1) is
generally intractable for larger graphs, authors propose a message-passing
algorithm for finding an approximate solution. It initializes $P^{0}(v)$ to
one for vertices in the blacklist and zero for vertices outside of it, and
updates the solution in each iteration as
$P^{t}(v_{i})=\sum\limits_{v_{j}\in\mathcal{N}(v_{i})}w_{ij}(P^{t-1}(v_{j})-C^{t-1}(v_{i},v_{j})),$
(7.2) where $C^{t-1}(v_{i},v_{j})$ is the portion of $P^{t-1}(v_{j})$
propagated from vertex $v_{i}$ in the previous step. In each iteration, $P(v)$
for the blacklisted vertices is set to $1$ again to reinforce the signal. In
[Kazato et al., 2016], PTP was successfully used to infer malicious domains
from a bipartite DNS graph. To the best of our knowledge, the most scaled and
comprehensive use of PTP is in [Jusko, 2017], where it was used for malicious
domain discovery and successfully implemented in the pipeline of Cisco
Cognitive Intelligence products. Unlike in most prior art, the problem of
constructing the graph from binary relations and setting weights to edges is
discussed in [Jusko, 2017]. We consider methods from [Jusko, 2017] state of
the art and compare the proposed method to them.

###  7.2 Method description

Our method for blacklist detection is based on recent progress in both
cybersecurity and machine learning. All bipartite graphs representing
relations on the input are first transformed into unipartite graphs and
message-passing inference with HMill models is performed afterwards. Since
messages are learned, this leads to more general procedure and unlocks
possible applications outside of the domain as well, because the algorithm can
be optimized for a particular problem and circumstances. This part of the
thesis describes an extension of our work presented in [Mandlík, 2018], where
we learned messages for inference algorithm using only one bipartite graph on
the input. In this thesis, we expand the method to work with several bipartite
graphs with heterogeneous vertices and edges. Also, thanks to HMill, the
introduction of a rather complicated method seems natural. Since the resulting
model is implemented as an HMill model and raw data (features) can be
therefore specified in both vertices and edges, the method effectively bridges
the gap between the classifier-based and the graph-based approaches. We
elaborate on this matter later after the method is explained.

####  7.2.1 Input

Formally, the method expects $n$ undirected bipartite graphs
$\mathopen{}\mathclose{{}\left(\mathcal{G}_{1},\mathcal{G}_{2},\ldots,\mathcal{G}_{n}}\right)$,
where $\mathcal{G}_{i}$ is a tuple
$\mathopen{}\mathclose{{}\left(\mathcal{A}_{i}\cup\mathcal{B}_{i},\mathcal{E}_{i}}\right)$.
Here, all vertices are split into two disjunct sets $\mathcal{A}_{i}$ and
$\mathcal{B}_{i}$, such that for all $i=1,\ldots,n$ and all
$\mathopen{}\mathclose{{}\left\\{u,v}\right\\}\in\mathcal{E}_{i}$, it holds
$u\in\mathcal{A}_{i}$ and $v\in\mathcal{B}_{i}$. Each vertex in the graph
represents a specific object of some type in real-world. For instance, in the
cybersecurity domain, a vertex may correspond to a second-level domain (SLD),
a client, an IP address, an email address, and so forth. Moreover, we assume
that one vertex set in each of the graphs contains vertices of the same
constant type, that is, each vertex set $\mathcal{A}_{i}$ is a subset of some
superset $\mathcal{A}$ and $\operatorname{type}(v)$ is constant for each
$v\in\mathcal{A}$. It may happen that some vertex $v\in\mathcal{A}$ is present
in two different bipartite graphs, in other words,
$\mathcal{A}_{i}\cap\mathcal{A}_{j}\neq\varnothing$. For vertices in second
partites we only assume that for each $i$, $\operatorname{type}(v)$ is
constant across all vertices in $\mathcal{B}_{i}$, however, this type may vary
with different $i$. Each edge represents some relationship between two objects
corresponding to its incident vertices, and each graph $\mathcal{G}_{i}$ may
differ in a definition of the relationship in real-world. In Figure 7.1, there
is an example of two such bipartite graphs that can be used to describe the
behavior of entities in a network. In this example, left vertex sets of both
graphs contain vertices representing SLD s appearing in network traffic. The
right set in the first graph comprises clients connected to the network
(computers, smartphones, printers, and others), and each edge
$\mathopen{}\mathclose{{}\left\\{d,c}\right\\}$ encodes the fact, that a
client $c$ connected to a domain $d$ during the time when the network was
observed. In the second graph, the other set represents individual binary
files identified by their hashes, and each edge
$\mathopen{}\mathclose{{}\left\\{d,b}\right\\}$ indicates that a process
instance of $b$ has issued a request to a domain $d$.

(a) $\mathcal{G}_{1}$

(b) $\mathcal{G}_{2}$

Figure 7.1: Example of two bipartite graphs, a domain-client graph
$\mathcal{G}_{1}$ (7.1(a)) and a domain-binary graph $\mathcal{G}_{2}$
(7.1(b)), that can be used to describe the behavior of entities in a network.
Clients are identified by a four-character-long identifier and binary files by
their sha256 hash (shortened in the figure). See text for further details.
Additionally, we assume that a blacklist of malicious objects
$L\subseteq\mathcal{A}$ is available. In our example, this means that we know
some domains in graphs exhibit malicious activity. The task is to discover
additional domains that are not present in list $L$, based on the behavioral
information encoded in input bipartite graphs.

####  7.2.2 Graph transformation

In the next step, each graph
$\mathcal{G}_{i}=(\mathcal{A}_{i},\mathcal{B}_{i})$, is transformed to a
unipartite transformed graph $G_{i}=(\mathcal{A}_{i},E_{i})$ with vertices
equal to vertices in the left vertex set $\mathcal{A}_{i}$ and edges
constructed with the following rule: $\\{u,v\\}\in E_{i}\iff\exists
b\in\mathcal{B}_{i}\colon\\{\\{u,b\\},\\{v,b\\}\\}\subseteq\mathcal{E}_{i},$
(7.3) for all $u,v\in\mathcal{A}_{i}$. Thus, in the domain-client graph
example, the transformed graph contains an edge between two domains if and
only if there exists a client that communicated with both of them.
Analogically, in the domain-binary graph, each edge represents that at least
one binary file issued a request to both domains involved. The semantics of
this transformation were already discussed in [Liu et al., 2014, Jusko, 2017,
Manadhata et al., 2014]. For example, when a computer is infected over a
network, the malicious file is downloaded from some server. Once the malware
becomes active, the vital step for attackers is to establish a communication
channel to the infected machine. This usually means contacting another server,
which is called a command-and-control (C&C) server. As the infected client
connected to both servers, an edge connecting them together appears in the
transformed graph. The transformation of the domain-client and domain-binary
graphs from Figure 7.1 can be found in Figure 7.2.

(a) $\mathcal{G}_{1}\mapsto G_{1}$

(b) $\mathcal{G}_{2}\mapsto G_{2}$

Figure 7.2: Examples of graph transformation performed on both graphs from
Figure 7.1. Each client and binary file creates a clique in the transformed
graph of the same cardinality as its degree in the original bipartite graph.
Consequently, $b_{1}$ and $b_{5}$ in (7.2(b)) do not contribute at all. Each
edge in transformed graphs is labelled with all clients/binaries that caused
its appearance.

####  7.2.3 Message-passing phase

After performing graph transformation on all $n$ bipartite graphs, we are left
with $n$-tuple $\mathopen{}\mathclose{{}\left(G_{1},\ldots,G_{n}}\right)$,
where each $G_{i}$ is a unipartite graph of vertices of a single type, in our
case SLD s. Assume that we want to infer an unknown quantity for a vertex
$v\in\mathcal{A}$, for example, estimate the probability of $v$ being
malicious. In each of the $n$ transformed graphs, we perform $k$-step HMill-
based message-passing inference centered in the current vertex $v$ exactly in
the same way as described in the previous chapter. Recall that neighborhoods
are modelled with (weighted) bag model nodes and knowledge extracted from
heterogeneous multi-modal sources is merged together with product model nodes.
In this case, however, we are dealing with not one, but $n$ graphs at once.
This is again easily solved with the framework—as inference in each of the
transformed graphs is implemented by an HMill model, we create a product model
node as a new root of the model tree with each of the $n$ inferential models
as its child. Recall that when evaluated, each HMill model returns a vector.
We can therefore regard evaluation of the resulting product model node as
first distilling a vector describing interactions of $v$ in each of the
available relations and merging them afterwards to obtain a result. The
process of how a sample tree is obtained for our case where $n=2$ is sketched
in Figure 7.3. Given a vertex $v$ for prediction, all input bipartite graphs
are transformed to unipartite graphs. The neighborhood of $v$ in each graph is
represented in $n$ bag data nodes, which are merged together with a product
data node afterwards. Again, each instance in these bag data nodes represents
one edge $\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$ in the transformed
graph, which is depicted with different colors in the figure. Since $u$ is
also a part of this adjacency relation, we could describe the edge (instance)
in the bag node as a deeper subtree and perform multiple-step inference, as
explained in the last chapter. Also, information (features) on both vertices
and edges can be included, which we have already discussed as well. Figure
7.3: A process of obtaining an HMill sample tree for a specific vertex $v$.
Neighbors of $v$ in unipartite graphs are colored. We do not show the whole
transformed graph, but only its subgraph induced by vertices in the
neighborhood of $v$. This example shows the same graphs as examples in Figures
7.1 and 7.2. Vertex $v$ corresponds to albertlee.biz domain, aliased as
$d_{4}$ and $d_{3}$ in Figure 7.2. See text for details. In this section, we
show how to deal with circumstances when no additional information about
vertices and edges is known, and input to the method comprises only $n$
bipartite graphs representing binary relations between vertices. To describe
edge $\mathopen{}\mathclose{{}\left\\{u,v}\right\\}$ where
$u,v\in\mathcal{A}_{i}$ in transformed graph $G_{i}$, we can use any
information about $u$, $v$ or the edge itself. Recall, that when performing