# Pickup and Delivery Problem with Transfers

Santiago Hincapie Potes , Catalina Lesmes Ramirez Mathematical Engineering
Student<EMAIL_ADDRESS>Universidad EAFIT,
Medellín, Colombia.Mathematical Engineering Student<EMAIL_ADDRESS>Universidad EAFIT,
Medellín, Colombia.

###### Abstract

In this article we will be talking about Pickup and Delivery Problems with
Transfers (PDP-T), the main idea is to show first the Pickup and Delivery
Problem (PDP), which is the problem from which the PDP-T derives from, later
show a description of the problem and the mathematical formulation to get a
clear view of what we are going to solve, which leads us to present different
ways of solving this kind of problems.
Keywords:Delivery problem, local search, random search, heuristics, hybrid
methods

##  1 Introduction

In many logistics problems we require to pick some products in one place and
take them to another place, which is why we define the pickup and delivery
problem (PDP). In this problem a set of vehicles pick up and deliver a set of
items. The goal is to deliver the items at the lowest cost while obeying a set
of constraints, such as time windows and capacities. The PDP is a well-
studied, NP-hard problem, so approximation algorithms and heuristics have been
developed to address variants of the PDP. There are many techniques we can use
to solve the PDP problem, such as genetic algorithms, various metaheuristics,
taboo search heuristics and branch and cut algorithms.

To the PDP certain transfers can be added, which means that the vehicle which
is delivering the product transfers it to another vehicle before the delivery
is done. By adding those transfers we have the PDP with Transfers (PDP-T), in
which we consider transferring items between vehicles. We can convert the PDP
problem into a PDP-T problem by adding some variables and constraints.

This article is organized as follows: In Section 1 an introduction is
presented which gives a general notion of the problem, later in Section 2 we
give a brief literature review were some examples of solution methods for the
PDP-T are presented, after in Section 3 we give the problem description, were
all the variables, parameters and constraints of the mathematical model are
explained in detail, next in Section 4 the solution algorithms we give to the
problem are explained, later in Section 5 we present the results of the
computational experimental and finally in Sections 6 and 7 we give some
conclusions and future work suggestions.

##  2 Literature Review

As we saw before, there are different methods to solve the PDP-T.

As an example a branch and cut method using Benders Decomposition can be used.
In this method, the set of constraints is decomposed into pure integer and
mixed constraints, and then a branch-and-cut procedure is applied to the
resulting pure integer problem, by using real variables and constraints
related as cut generators. The key on the success of this method is that those
constraints defined by a logical sentence are not modeled using the big-M
technique, as usual in a branch-and-bound methodology behind the original PDP
formulation. This method may be applied only when the objective function is
either pure real or pure integer (Cortés et al. , , 2010).

Another method which solves the problem is a Very Large Neighborhood Search
with Transfers (VLNS-T) based on the Adaptive Very Large Neighborhood Search
(VLNS). The VLNS algorithm uses simulated annealing to randomly choose
neighboring schedules and iteratively improve the schedule. Neighboring
schedules are formed by removing random items and reinserting them with
heuristics. VLNS-T is based on the VLNS algorithm for the PDP without
transfers, a variant of simulated annealing in which the neighborhood of
states is ”very large”. In this case we remove random items from the schedule
and then reinsert them with multiple heuristics to find neighbors (Cat, ,
n.d.).

There are also some special cases of the PDP-T, as the Pickup and Delivery
Problem with Shuttle Routes (PDP-S) which relies on a structured network with
two categories of routes. Pickup routes visit a set of pickup points
independently of their delivery points and end at one delivery point. Shuttle
routes are direct trips between two delivery points. Requests can be
transported in one leg (pickup route) or two legs (pickup route plus shuttle
route) to their delivery point. The PDP-S applies to transportation systems
with a multitude of pickup points and a few, common delivery points (Masson et
al. , , 2010).

##  3 Problem Description and Mathematic Formulation

In the PDP, a heterogeneous vehicle fleet based at multiple terminals must
satisfy a set of transportation requests. Each request is defined by a pickup
point, a corresponding delivery point, and a demand to be transported between
these locations (Parragh et al. , , 2008a). Now, PDP-T allows the option for
passengers to transfer between vehicles, provided that the locations of the
transfer points are fixed and known.

###  3.1 Problem Formulation

Let $N$ be the set of transportation requests. For each transportation request
$i\in N$ a load of size $\bar{q_{i}}\in\mathbb{N}$ has to be transported from
a set of origins $N^{+}_{i}$ to a set of destinations $N_{i}^{-}$.
Each load is divided as follows $\displaystyle\bar{q_{i}}=\sum_{j\in
N_{i}^{+}}q_{j}=-\sum_{j\in N_{i}^{-}}q_{j}$, i.e., positive quantities for
pickups and negative quantities for deliveries.
Define $\displaystyle N^{+}:=\cup_{i\in N}N_{i}^{+}$ as the set of all origins
and $\displaystyle N^{-}:=\cup_{i\in N}N_{i}^{-}$ as the set of all
destinations. Let $V:=N^{+}\cup N^{-}$.
Furthermore, let $M$ be the set of vehicles. Each vehicle $k\in M$ has a
capacity $Q_{k}\in\mathbb{N}$ a start location $k^{+}$ and an end location
$k^{-}$.
Define $M^{+}:=\\{k^{+}|k\in M\\}$ as the set of start locations and
$M^{-}:=\\{k^{-}|k\in M\\}$ as the set of end locations. Let $W:=M^{+}\cup
M^{-}$.
For all $i,j\in V\cup W$ let $d_{ij}$ denote the travel distance, $t_{ij}$ the
travel time and $c_{ij}$ the travel cost. Note that the dwell time at origins
and destinations can be easily incorporated in the travel time and therefore
will not be considered explicitly.

To formulate the PDP as a mathematical program we introduce four types of
variables: for $i\in N$ and $k\in M$:

$z_{i}^{k}=\begin{cases}1&\text{if transportation request }i\text{ is assigned
to vehicle }k\\\ 0&\text{In other case}\end{cases}$

for $(i,j)\in(V\times V)\cup\\{(k^{+},j)|j\in V\\}\cup\\{(j,k^{-})|j\in V\\}$
and $k\in M$:

$x_{ij}^{k}=\begin{cases}1&\text{if vehicle }k\ \text{travels from location
}i\ \text{to location }j\\\ 0&\text{In other case}\end{cases}$

$D_{i}$ with $(i\in V\cup W)$, specifying the departure time at vertex $i$ and
$y_{i}$ with $(i\in V\cup W)$, specifying the load of the vehicle arriving at
vertex $i$ (Parragh et al. , , 2008a).
All this information can be summarized in Table 1:

Table 1: Parameters and variables found in (Parragh et al. , , 2008a) for the model Object | Meaning
---|---
$M$ | Set of vehicles
$C$ | Set of requests
$T$ | Set of transference point
$M^{+}$ | Set of origin depots for vehicles
$M^{-}$ | Set of destination depots for vehicles
$N^{+}$ | Set of origin nodes for requests
$N^{-}$ | Set of destination nodes for requests
$N$ | Set of nodes associated with requests
$V$ | Set of nodes
$q_{ij}$ | Size of request $i\in C$
$Q_{k}$ | Capacity of vehicle $k\in K$
$t_{ij}$ | Minimum ride time from node $i$ to node $j$
$c_{ij}$ | The travel cost
$d_{ij}$ | The travel distance
$z_{i}^{k}$ | bind transportation request and vehicles
$x_{ij}$ | bind routes and vehicles information
$D_{i}$ | specifying the departure time at specific vertex
$y_{i}$ | specifying the load of the vehicle arriving

Now the model is:

minimize $\displaystyle\sum_{i,j\in V\cup W}dij$ (1) subject to
$\displaystyle\sum_{k\in M}z_{i}^{k}=1$ for all $\displaystyle i\in N$ (2)
$\displaystyle\sum_{j\in V\cup W}x_{il}^{k}=z_{i}^{k}$ for all $\displaystyle
i\in N,l\in N^{+}_{i}\cup N_{i}^{-}k\in M$ (3) $\displaystyle\sum_{j\in
V\cup\\{k^{-}\\}}x^{k}_{k^{+}j}=1$ for all $\displaystyle k\in M$ (4)
$\displaystyle\sum_{j\in V\cup\\{k^{+}\\}}x^{k}_{k^{-}j}=1$ for all
$\displaystyle k\in M$ (5) $\displaystyle D_{k^{+}}=0$ for all $\displaystyle
k\in M$ (6) $\displaystyle D_{p}\leq D_{q}$ for all $\displaystyle i\in N,p\in
N_{i}^{+}q\in N_{i}^{-}$ (7) $\displaystyle x_{ij}^{k}=1\Rightarrow
D_{i}+t_{ij}\leq D_{j}$ for all $\displaystyle i,j\in V\cup W,k\in M$ (8)
$\displaystyle y_{k^{+}}=0$ for all $\displaystyle k\in M$ (9) $\displaystyle
y_{l}\leq\sum_{k\in M}Q_{k}z_{i}^{k}$ for all $\displaystyle i\in N,l\in
N_{i}^{+}\cup N_{i}^{-}$ (10) $\displaystyle D_{i}\geq 0$ for all
$\displaystyle i\in V\cup W$ (11) $\displaystyle y_{i}\geq 0$ for all
$\displaystyle i\in V\cup W$ (12)

Constraint (2) ensures that each transportation request is assigned to exactly
one vehicle. By constraint (3) a vehicle only enters or leaves a location $l$
if it is an origin or a destination of a transportation request assigned to
that vehicle. The next (4) and (5) make sure that each vehicle starts and ends
at the correct place. Also the constraints number (6), (7), (8) and (11)
together form the precedence constraints the others together form the capacity
constraints. Constraints (9), (10) and (12) together form the capacity
constraints.
This mathematical program, models the PDP. Now for a PDP-T model in the
literature it can bee seen that the transfer point is introduced, and the
procedure for the extended model is to iterative add constrains that involve
this transfer points. The process of inserting the transfers will be described
in Section 4.

##  4 Solution Algorithms

In this section we present solution algorithm for the PDP-T. In general, for
this part, we first solve the PDP using different methods (which would be
described in this section) and then we add the transship option to the PDP
solution, transship option always increases the computational time but makes
the solution better. In subsection 4.1 we present two constructive methods: a
greedy approach and multistart heuristic in Takoudjou et al. ,  (2012), then
in subsection 4.2 we present a random method GRASP which is applied to the
constructive methods mentioned before, later in subsection 4.3 we present
three local search methods: Variable neighborhood descent, Adaptive large
neighborhood search and simulated annealing, finally in subsection 4.4 we
present two hybrid methods which combine local search methods and genetic
algorithms.

###  4.1 Constructive Methods

####  4.1.1 Greedy approach

We implemented the greedy approach described in Coltin et al. ,  (2014), which
basically iterates through every item and vehicle and then inserts the best
pickup and delivery action into the schedule, it always chooses the option
that increases the least the cost, we repeat the process until no unassigned
items remain. Then we add the transfer option with another greedy idea,
inspired on Clarke saving algorithm, this method computes the savings of the
transfer, and it improves the solution, then we modify the solution by adding
this new note.

####  Clarke and Wright Algorithm

The Clarke and Wright savings algorithm is one of the most known heuristic for
VRP. It applies to problems for which the number of vehicles is not fixed (it
is a decision variable), and it works equally well for both directed and
undirected problems. When two routes $(0,...,i,0)$ and $(0,j,...,0)$ can
feasibly be merged into a single route $(0,...,i,j,...,0)$, a distance saving
$s_{ij}=c_{i0}+c_{0j}-c_{ij}$ is generated. The algorithm works as it follows:

Step 1. Savings computation

  * •

Compute the savings $s_{ij}=c_{i0}+c_{0j}-c_{ij}$ for $i,j=1,...,n$ and $i\neq
j$.

  * •

Create n vehicle routes $(0,i,0)$ for $i=1,...,n$.

  * •

Order the savings in a non increasing fashion.

Step 2. Route Extension (Sequential version)

  * •

Consider in turn each route $(0,i,...,j,0)$.

  * •

Determine the first saving $s_{ki}$ or $s_{jl}$ that can feasibly be used to
merge the current route with another route ending with $(k,0)$ or starting
with $(0,l)$.

  * •

Implement the merge and repeat this operation to the current route.

  * •

If not feasible merge exists, consider the next route and reapply the same
operations.

  * •

Stop when not route merge is feasible.

####  4.1.2 Hybrid Multistart Heuristic

The PDPT is a hard combinatorial optimization problem. The heuristics for this
problem must avoid being trapped by local optimum. To overcome local
optimality, a diversification procedure is needed. So we use the transshipment
heuristic from the hybrid heuristic for the PDP found in Takoudjou et al. ,
(2012) to obtain a solution to the PDP-T.
At each iteration, the transshipment heuristic is used to improve the PDP
solution and obtain a solution to the PDP-T. From the current PDP solution,
each request ($i$, $i+n$) $\in$ R; is removed from the solution. Then, ($i$,
$i+n$) is split into two different requests ($i$, $e_{t}$) and ($s_{t}$,
$i+n$), where $e_{t}$ and $s_{t}$ are the inbound/outbound doors of a
transshipment point $t$ $\in$ T. The best reinsertion cost of ($i$, $i+n$) in
the solution is computed. The best insertion cost to insert ($i$, $e_{t}$)
following by insertion of ($s_{t}$, $i+n$) in the solution is computed. The
cost to insert ($s_{t}$, $i+n$) following by the insertion of ($i$, $e_{t}$)
at their best position is computed. Between the three possibilities, the
insertion or the reinsertions offering the minimum cost is performed. As
described in Takoudjou et al. ,  (2012), this hybrid approach merges many
ideas. The solving principle is based on the following steps:

  * •

An initial PDP solution (Sol) is calculated by gradually inserting the
requests in routes associated to vehicles.
The initial solution is consisting of a single route containing only the
starting point and the ending point of each route. Requests are then
successively introduced into the tour of a vehicle offering the minimal
increase of the cost of transport.

  * •

Transshipment is then used to “destroy” and “repair” PDP solution to obtain a
better solution (PDPT solution).
At each iteration, the transshipment heuristic is used to improve the PDP
solution and obtain a solution to the PDPT. From the current PDP solution,
each request; is removed from the solution. Then, the request is split into
two different requests, the first go from the start point to the best transfer
point, the other go from the transfer point to the ending point.

###  4.2 Random Search Method

####  4.2.1 GRASP

GRASP is a multi-start metaheuristic for combinatorial optimization problems,
in which each iteration consists of two phases: construction and local search.
The construction phase builds a feasible solution, whose neighborhood is
investigated until a local minimum is found during the local search phase. The
best overall solution is kept as the result. Here we will only use one phase,
which is the construction of the solution Resende & Ribeiro,  (2010).

###  4.3 Local Search Methods

####  4.3.1 Variable Neighborhood Descent

The VND is used to explore the neighborhood of the current solution, which is
based on three operators defined below: ADR, RNR and SWR. With these three
operators, we can explore the solution space more intensively. An operator is
a move that transforms one solution to another with small modifications. In
best improvement operators the whole neighborhood is analyzed and the best
solution is kept. In first improvement local search the first better solution
found in the neighborhood is kept (Takoudjou et al. , , 2012).

##### SWR

(Swap requests between routes) In SWR neighborhood, two requests belonging to
two different routes are exchanged together provided that all PDP constraints
are satisfied.

##### RNR

(Remove and insert a request) In RNR neighborhood, provided that all PDP
constraints are satisfied, request belonging to one route is removed and
inserted in another route.

##### ADR

(Advance or delay a request) In ADR neighborhood, requests belonging to one
given route is advanced or delayed in the same route if PDP constraints are
all satisfied.

####  4.3.2 Adaptive Large Neighborhood Search

The general idea is to remove some requests from routes in the current
solution and then reinsert them elsewhere to arrive at a new solution. At each
iteration, one each of several removal and insertion heuristics is randomly
selected based on their designated weights. Letting the set of removal
heuristics be ${rh_{1},...,rh_{nr}}$ and their corresponding weights be
${rw_{1},...,rw_{nr}}$, $rh_{i}$ is chosen with probability
$\frac{rw_{i}}{\sum_{j=1}^{nr}rw_{j}}$, where $nr$ is the number of options.
Similarly for the insertion heuristics ${ih_{1},...,ih_{ni}}$ and their
corresponding weights ${iw_{1},...,iw_{ni}}$, the probability of $ih_{i}$
being $\frac{ih_{i}}{\sum_{j=1}^{n}ih_{j}}$ selected is, where $ni$ is the
number of options. After each iteration, if the solution realized is feasible
and better than the incumbent, the latter is updated. In any case, the process
is repeated until a predefined number of iterations $n_{II}^{max}$ is reached
(Masson et al. , , 2012).

####  4.3.3 Simulated Annealing

Simulated annealing is a metaheuristic that begins at some state, and chooses
a random “neighbor” of that state. With probability accept $(e,e_{0},t)$ the
new state is accepted as the current state, where $e$ is the “energy” of the
current state, $e_{0}$ is the energy of the new state, and $t$ is the
temperature, or the fraction of iterations of the algorithm currently
completed. If the new state is rejected we remain at the current state and
repeat with a new neighbor. The algorithm continues either for a fixed number
of iterations or until the energy crosses some threshold, when the best
solution that has been encountered thus far is returned (Coltin & Veloso, ,
2014).

###  4.4 Hybrid Methods

In this subsection we present two hybrid methods, both of this methods are
based on the use of a local search methods and genetic algorithms but they
possess different hybridizations. The first algorithm uses simulated annealing
for the population selection and for its mutation variance of the probability,
this is made with the purpose of achieves a better convergence, next to a
diverse population. The second hybrid method uses the taboo list to avoid
mating of chromosomes that have low Minkowski in order to keep diversity, also
use this information in order to change mutation probability. We show in
precise detail how the genetic algorithm was made, so in subsection 4.4.1 we
show the encoding of the problem, then we describe the initial population in
subsection 4.4.2, which we got in two different ways, we have a totally
randomized initial solution and a constructive solution, later we select the
population to be matched which is described in subsection 4.4.3, then in
subsection 4.4.4 we describe the crossover between the population, after that,
in subsection 4.4.5 we describe the mutation for both first and second
methods, then in subsection 4.4.6 we explain how the new population is chosen
for the first method and the second method, and finally in subsection 4.4.7 we
explain the stop criterion.

####  4.4.1 Encoding

Given the structure of the problem, the solutions can be represented as a set
of $N$ paths, one for each vehicle, considering that it is important to
consider the transfers as much as the possibility of a node recurrence, it was
opted to represent the solution as a directed multi-graph, in which the set of
vertexes represent some of the visited nodes and there exists a connexion
between two vertexes if one is a successor from the other, the direction goes
in the way of the route, additionally each one of the multi-sets is found
reserved for each route. To represent each multi-graph, a multidimensional
matrix was used, which positions $i,j,k$ possessed the locality distance $i$
to the locality $j$, if $j$ is the successor of $i$ in route $k$. For the
processes of the genetic algorithm to each element of the solution matrix a
derivate from the logistic function was applied (fixing the value from $0$ to
$0$), meaning that to each element a function as is shown below was applied:

$f(x)=\begin{cases}\frac{e^{-x}}{(1+e^{-x})^{2}}&\text{if }x>0\\\
0&\text{otherwise}\end{cases}$

this function was applied with the purpose of giving more importance to the
short routes, in this way the intention is to preserve this property. The
complexity in space of this encoding is $O(N|path|^{2})$.

####  4.4.2 Initial Population

For the computational experiments we performers with two district initial
population methods, a totally randomized approach and a set of district greedy
random solutions.

#### Totally Randomized Initial Solution

A set of 256 random multi-graphs were generated using the model of a
modification of the algorithm from Erdős-Rényi as described in (Godehardt, ,
1993). It was opted to generate random graphs in place of random matrices, due
to the fact that the solutions generated this way are usually more feasible
with a higher probability, additionally they use to be better.

#### Constructive Solution

We construct a set of 256 district random greedy solution using the algorithm
described in the section 4.2.

####  4.4.3 Selection

For the selection part we use Baker’s stochastic universal selection (SUS)
(Baker, , 1987), this algorithm uses a single random value to sample all the
solutions by choosing them at evenly spaced intervals. This gives weaker
members of the population (according to their fitness) a chance to be chosen
and thus reduces the unfair nature of fitness-proportional selection methods,
also the experimental work by Hancock (Hancock, , 1994) clearly demonstrates
the superiority of this approach.

#### Fitness Function

For a $K$ matrix its fitness function is the sum of the total distance of each
one of its routes, in case that the solution is not feasible (it does not
satisfy all the request) a $\infty$ value is assigned.

####  4.4.4 Crossover

To make the crossover the solutions which were previously selected were
matched and a new matrix was generated in which its entries were uniform
random numbers. The resultant son matrix possesses the gen of its first parent
in the position $i,j,k$ if the generated random matrix possesses a lower
number to $0.5$ in the corresponding position, in the other cases it possesses
the gen of its other parent.

####  4.4.5 Mutation

##### First method

For mutation operator we use a shift operator, that shifts the value of two
elements of the matrix depending on a parameter that we varied according to
the diversity in the population (measured in terms of the fitness variance)
and actual solution temperature (describe in the next step).

##### Second Method

For mutation operator we also use a shift operator but in this method we uses
a fixed probability $\alpha=0.05$

####  4.4.6 New population

##### First method

For the actualization of the population, we perform a probabilistic
replacement test, in this approach the solutions are accepted into the new
generation using a Boltzmann type calculation. At each generation the
temperature is decreased, thereby decreasing the probability of accepting
higher energy individuals.

##### Second method

For the population update, we compute the Minkowski distance between
population members and uses taboo list to avoid nears solutions, in order to
keep population diversity.

####  4.4.7 Stop criterion

For the stop criterion we compared the variance between the solutions of the
last $4$ generations, this with the purpose of ending the algorithm when it
reaches an stationary state.

##  5 Computational Experimentation

In this section we show the computational results given by the algorithms
described in section 4 that solve the PDP-T. In subsection 5.1, we present
tables with the instance size, the number of vehicles, optimal solutions and
the computational time, then we compare them to see which is the best
heuristic used.

All algorithms were coded on python 3.5.3 running on Intel Xeon E5-1620V4 ,
64-bit processor with 32 gigabyte RAM and Ubuntu 17.10. We used problem
instances from those in for the PDP, which are related to the well-known
Solomon instances. The datasets are available at
http://www.sintef.no/Projectweb/TOP/PDPTW/ . For each instance we introduce a
random number of transfer points, in a random node, to become the original PDP
instance in a PDP-T instance.

Each data set contains $X-Y$ hence, the lengths of the arcs in the network are
$L_{2}$ distances. For the locations of the vehicle depots, the origin and
final depots for each vehicle were randomly generated, scattered over the
region formed by the nodes. The number of vehicles are also random generated,
and finally we ignore the time windows restriction because are not consider in
the model.

###  5.1 Computational results

In this subsection we show the results of all the heuristics mentioned before
and compare them.

####  5.1.1 Constructive Methods

In this subsection we show the tables of the computational results for the
constructive methods, in Table 2 we present the results for the Greedy
approach and in Table 3 we show the results for the MULTISTART method.

Table 2: Greedy approach Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 2001.31 | 28.23
lc2_2_3 | 200 | 6 | 4920.43 | 50.48
lc2_4_5s | 400 | 12 | 7235.92 | 123.43
LC2_8_5 | 800 | 25 | 226748.20 | 240.12
Table 3: MULTISTART Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1825.83 | 40.81
lc2_2_3 | 200 | 6 | 4378.92 | 76.21
lc2_4_5s | 400 | 12 | 6891.06 | 173.79
LC2_8_5 | 800 | 25 | 201234.73 | 301.83

As we can see, the MULTISTART method gives better optimal solutions for the
PDP-T, but at the same time it requires more computational time.

####  5.1.2 Random Search Method

In this subsection the results for the GRASP which is the random search method
are shown, as we mentioned in Section 4, we add GRASP to the constructive
methods. In Table 4, the greedy approach with GRASP is shown and in Table 5 we
have the results for the MULTISTART with GRASP.

Table 4: Greedy approach (GRASP) Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1713.01 | 54.12
lc2_2_3 | 200 | 6 | 4001.12 | 80.15
lc2_4_5s | 400 | 12 | 6481.10 | 196.92
LC2_8_5 | 800 | 25 | 198917.87 | 320.16
Table 5: MULTISTART (GRASP) Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1521.31 | 58.63
lc2_2_3 | 200 | 6 | 3832.12 | 89.53
lc2_4_5s | 400 | 12 | 6312.18 | 200.12
LC2_8_5 | 800 | 25 | 190143.12 | 354.10

If we compare Tables 4 and 5 we can see that the optimal solutions are given
by the MULTISTART (GRASP), but it also takes more computational time than the
Greedy approach (GRASP). This results are coherent, since they follow the same
behavior as we mentioned in subsection 5.1.1, but we can also see that when we
add GRASP to our constructive methods the optimal solutions are better, even
thought they take more computational time.

####  5.1.3 Local Search Methods

In this subsection we show tables comparing all the results, plus one more,
the other is called ”MIX”, this table refers to a union of the VND and ALNS
methods, we wanted to see what happened if we merged two of them. So in Table
6 we present the ALNS, in Table 7 the VND, in Table 8 the SA and in table 9
the mix between the VND and the ALNS. Also in this subsection we also wanted
to show a 3D plot, showing the results of all the heuristics including the
initial solution, this plot is shown in Figure 1.

Table 6: ALNS Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1154.21 | 160.32
lc2_2_3 | 200 | 6 | 3213.31 | 233.21
lc2_4_5s | 400 | 12 | 5421.32 | 312.23
LC2_8_5 | 800 | 25 | 182313.32 | 600.12
Table 7: VND Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1032.12 | 187.21
lc2_2_3 | 200 | 6 | 3102.23 | 254.12
lc2_4_5s | 400 | 12 | 5250.12 | 352.21
LC2_8_5 | 800 | 25 | 173213.12 | 721.48
Table 8: SA Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1321.23 | 140.32
lc2_2_3 | 200 | 6 | 3543.34 | 221.34
lc2_4_5s | 400 | 12 | 5732.12 | 301.31
LC2_8_5 | 800 | 25 | 182313.32 | 572.12
Table 9: MIX Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 1023.62 | 350.14
lc2_2_3 | 200 | 6 | 3102.23 | 500.51
lc2_4_5s | 400 | 12 | 5291.17 | 698.92
LC2_8_5 | 800 | 25 | 177843.01 | 1425.92

size020406080100time050010001500200025003000objective020040060080010001200ALNSVNSSAMIXInitial
Solution

Figure 1: Comparison of the heuristic methods

In Figure 1, we can see in the plot the time, the objective function and the
size of the instance in contrast between each other, is a comparison of all
methods and the initial solution.

We can see from Tables 6 through 9 and Figure 1 that the best algorithm as
seen in computational time is the Simulated Annealing and that best algorithm
as seen in optimal solution is the mix between the ALNS and the VND. As until
now, our best results are given with the random search algorithms if we
compare them with the local search algorithms we can say that the best results
are given by the local search algorithms specifically with the mix of the ALNS
and the VNS. For this mix, we can see a very interesting result, even thought
it was the best result of all, we should expect it gave a much better result,
but the difference from the initial solution is not that wide.

####  5.1.4 Hybrid Methods

In this subsection we show the results given by the algorithms mentioned in
subsection 4.4, it is important to mention that the reaction of the initial
population and the associated individual operations such as encoding were
computed in parallel. So, in Tables 10 and 11 we show the results for the
genetic algorithms and simulated annealing with both the constructive and the
random initial solution, and in Tables 12 and 13 we show the results for the
genetic algorithms and taboo with both the constructive and the random initial
solutions.

Table 10: Genetic algorithm & SA, with random initial solution Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 720.60 | 789.12
lc2_2_3 | 200 | 6 | 2454.12 | 1325.21
lc2_4_5s | 400 | 12 | 4843.93 | 2012.32
LC2_8_5 | 800 | 25 | 14794.27 | 3024.15
Table 11: Genetic algorithm & SA, with Constructive solution Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 673.12 | 1254.65
lc2_2_3 | 200 | 6 | 2021.23 | 1743.54
lc2_4_5s | 400 | 12 | 4654.63 | 2433.33
LC2_8_5 | 800 | 25 | 13647.12 | 3502.12
Table 12: Genetic algorithm & Taboo, with random initial solution Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 690.43 | 823.43
lc2_2_3 | 200 | 6 | 2144.33 | 1513.43
lc2_4_5s | 400 | 12 | 4136.38 | 2243.43
LC2_8_5 | 800 | 25 | 13474.76 | 3324.69
Table 13: Genetic algorithm & Taboo, with Constructive solution Instance | Instance Size | Number of vehicles | Optimal Solution | Time
---|---|---|---|---
lc204 | 100 | 3 | 610.60 | 1242.43
lc2_2_3 | 200 | 6 | 1944.33 | 2091.35
lc2_4_5s | 400 | 12 | 4070.48 | 2695.49
LC2_8_5 | 800 | 25 | 12689.23 | 3820.94

When we compare the results given by Tables 10 through 13 we can see that the
constructive solution gives better results than the random initial solution
for both SA and Taboo, but they also take more computational time. In general
terms, the best solution as seen in optimal solution is given by the genetic
algorithm and Taboo with the constructive solution. The best solution as seen
in computational time is given by genetic algorithm and simulated annealing
with a random initial solution.

##  6 Conclusions

  * •

The best solution as seen in optimal solution is given by the genetic
algorithm and Taboo with the constructive solution in hybrid methods.

  * •

The best solution as seen in computational time is given by the greedy
approach in constructive methods.

  * •

An interesting matter with VND is that even though the solution is bad it has
a good computational time.

  * •

The mix between the ALNS and the VND is a kind of hybrid but it is slow and
gives bad solution, but the hybrids with the genetic algorithm give better
results due to the fact that they use the hybridization to cover its weakness.

  * •

We can see from all the computational times that when the instances are small,
the times between the heuristics are very similar but when the instances are
bigger, the times between the heuristics are very different.

##  7 Future Work

An interesting matter that it could be talked about in the future is about the
lower bound of this problem, because in none of the articles we have reviewed
the topic has risen up, and it’s important because we can not compare how much
closer we are getting to a better solution, we can see it gets better every
time but we can not see if we are closer to our goal. Also, it would be
interesting to run the algorithm several times rising the number of transfer
points with the purpose of seeing how it affects the computational time and
how much it diminishes the objective function.

## References

  * Cat, (n.d.)
  * Baker, (1987) Baker, James E. 1987.  Genetic Algorithms and their Applications (Hard).  Lawrence Erlbaum.
  * Bouros et al. , (2011) Bouros, Panagiotis, Sacharidis, Dimitris, Dalamagas, Theodore, & Sellis, Timos. 2011.  Dynamic Pickup and Delivery with Transfers.
  * Coltin & Veloso, (2014) Coltin, Brian, & Veloso, Manuela. 2014.  Online pickup and delivery planning with Transfers for mobile robots.  Proceedings - IEEE International Conference on Robotics and Automation, 5786–5791.
  * Coltin et al. , (2014) Coltin, Brian, Choset, Howie, & Smith, Stephen. 2014.  Multi-agent Pickup and Delivery Planning with Transfers.
  * Cortés et al. , (2010) Cortés, Cristián E., Matamala, Martín, & Contardo, Claudio. 2010\.  The pickup and delivery problem with transfers: Formulation and a branch-and-cut solution method.  European Journal of Operational Research, 200(3), 711–724.
  * Godehardt, (1993) Godehardt, Erhard A.J. 1993.  Probability Models for Random Multigraphs with Applications in Cluster Analysis.  Pages 93–108 of: Quo Vadis, Graph Theory? - A Source Book for Challenges and Directions.  Elsevier.
  * Hancock, (1994) Hancock, Peter J. B. 1994.  An empirical comparison of selection methods in evolutionary algorithms.  Pages 80–94 of: Evolutionary Computing.  Springer Berlin Heidelberg.
  * Li et al. , (2016) Li, Yuan, Chen, Haoxun, & Prins, Christian. 2016.  Adaptive large neighborhood search for the pickup and delivery problem with time windows, profits, and reserved requests.  European Journal of Operational Research, 252(1), 27–38.
  * Masson et al. , (2010) Masson, Renaud, Lehuede, Fabien, Peton, Olivier, & Ropke, Stefan. 2010.  The Pickup and Delivery Problem with Shuttle routes.  1–4.
  * Masson et al. , (2012) Masson, Renaud, Lehuédé, Fabien, & Péton, Olivier. 2012.  An Adaptive Large Neighborhood Search Heuristic for the Pickup and Delivery Problem with Time Windows.  Transportation Science, 47(3), 344–355.
  * Mladenović & Hansen, (1997) Mladenović, N., & Hansen, P. 1997.  Variable neighborhood search.  Computers & Operations Research, 24(11), 1097–1100.
  * Nanry & Wesley Barnes, (2000) Nanry, William P., & Wesley Barnes, J. 2000.  Solving the pickup and delivery problem with time windows using reactive tabu search.  Transportation Research Part B: Methodological, 34(2), 107–121.
  * Parragh et al. , (2008a) Parragh, Sophie N., Doerner, Karl F., & Hartl, Richard F. 2008a.  A survey on pickup and delivery problems.  Journal fur Betriebswirtschaft, 58(1), 21–51.
  * Parragh et al. , (2008b) Parragh, Sophie N., Doerner, Karl F., & Hartl, Richard F. 2008b.  A survey on pickup and delivery problems. Part I: Transportation between customers and depot.  Journal für Betriebswirtschaft, 58(2), 81–117.
  * Polat et al. , (2015) Polat, Olcay, Kalayci, Can B., Kulak, Osman, & Günther, Hans-Otto. 2015.  A perturbation based variable neighborhood search heuristic for solving the Vehicle Routing Problem with Simultaneous Pickup and Delivery with Time Limit.  European Journal of Operational Research, 242(2), 369–382.
  * Qu & Bard, (2012) Qu, Yuan, & Bard, Jonathan F. 2012.  A GRASP with adaptive large neighborhood search for pickup and delivery problems with transshipment.  Computers and Operations Research, 39(10), 2439–2456.
  * Resende & Ribeiro, (2010) Resende, Mauricio G.C., & Ribeiro, Celso C. 2010.  Greedy Randomized Adaptive Search Procedures: Advances, Hybridizations, and Applications.  Pages 283–319 of: Handbook of Metaheuristics.  Springer US.
  * Takoudjou et al. , (2012) Takoudjou, Rodrigue Tchapnga, Deschamps, Jean-christophe, & Dupas, Rémy. 2012\.  A hybrid multistart heuristic for the pickup and delivery problem with and without transshipment.  9th International Conference on Modeling, Optimization & SIMulation.


# MM-GNN: Mix-Moment Graph Neural Network towards Modeling Neighborhood
Feature Distribution

Wendong Bi  Institute of Computing Technology, Chinese Academy of
SciencesBeijingChina<EMAIL_ADDRESS>,  Lun Du  Microsoft Research
AsiaBeijingChina<EMAIL_ADDRESS>,  Qiang Fu  Microsoft Research
AsiaBeijingChina<EMAIL_ADDRESS>,  Yanlin Wang  Microsoft Research
AsiaBeijingChina<EMAIL_ADDRESS>,  Shi Han  Microsoft Research
AsiaBeijingChina<EMAIL_ADDRESS>and  Dongmei Zhang  Microsoft Research
AsiaBeijingChina<EMAIL_ADDRESS>
(2023)

###### Abstract.

Graph Neural Networks (GNNs) have shown expressive performance on graph
representation learning by aggregating information from neighbors. Recently,
some studies have discussed the importance of modeling neighborhood
distribution on the graph. However, most existing GNNs aggregate neighbors’
features through single statistic (e.g., mean, max, sum), which loses the
information related to neighbor’s feature distribution and therefore degrades
the model performance. In this paper, inspired by the method of moment in
statistical theory, we propose to model neighbor’s feature distribution with
multi-order moments. We design a novel GNN model, namely Mix-Moment Graph
Neural Network (MM-GNN), which includes a Multi-order Moment Embedding (MME)
module and an Element-wise Attention-based Moment Adaptor module. MM-GNN first
calculates the multi-order moments of the neighbors for each node as
signatures, and then use an Element-wise Attention-based Moment Adaptor to
assign larger weights to important moments for each node and update node
representations. We conduct extensive experiments on 15 real-world graphs
(including social networks, citation networks and web-page networks etc.) to
evaluate our model, and the results demonstrate the superiority of MM-GNN over
existing state-of-the-art models

Graph neural networks, Graph representation learning, Feature distribution,
Moment, Social networks

††journalyear: 2023††copyright: acmcopyright††conference: Proceedings of the
Sixteenth ACM International Conference on Web Search and Data Mining; February
27-March 3, 2023; Singapore, Singapore††booktitle: Proceedings of the
Sixteenth ACM International Conference on Web Search and Data Mining (WSDM
’23), February 27-March 3, 2023, Singapore, Singapore††price: 15.00††doi:
10.1145/3539597.3570457††isbn: 978-1-4503-9407-9/23/02††ccs: Computing
methodologies Neural networks††ccs: Information systems Social networks

##  1\. Introduction

Graph data is ubiquitous in the real world, such as social networks, citation
networks, etc. Graph Neural Networks (GNNs) have shown expressive power on
graph representation learning in recent years. Most existing GNNs follow the
neighborhood aggregation scheme, where each node aggregates features of
neighbors to update its own representation. Many GNNs (Kipf and Welling, 2016;
Veličković et al., 2017; Derr et al., 2020; Yang et al., 2016) with different
aggregators have been proposed to solve various downstream tasks, which can be
divided into graph-level tasks, edge-level tasks, and node-level tasks, such
as graph classification (Xu et al., 2018a; Bai et al., 2019), link prediction
(Zhang and Chen, 2018; Zhou et al., 2018; Guo et al., 2022), and node
classification (Dai et al., 2022; Lin et al., 2020, 2014; Fettal et al., 2022;
Jin et al., 2021; Zhao et al., 2021).

However, the aggregation schemes of existing GNNs usually use single
statistics (e.g., mean, max, sum). For example, GCN (Kipf and Welling, 2016)
uses mean aggregator, GraphSAGE (Hamilton et al., 2017) uses max aggregator,
and GIN uses sum aggregator. And some state-of-the-art GNNs (e.g. GCNII (Chen
et al., 2020), DAGNN (Liu et al., 2020)) still use single statistics, without
retaining the complete information of neighbor’s feature distribution.
Recently, researchers (Xu et al., 2018a; Du et al., 2021b; Ma et al., 2022;
Luo et al., 2022) begin to study the modeling of neighborhood distribution and
have discussed its importance for graph data from various perspectives.
However, the studies of neighbor’s feature distribution are limited to lower-
order characteristics without considering higher-order distribution
characteristics.

Considering that single statistic cannot represent complex distribution with
more than one degree of freedom, it’s intuitive that we should introduce more
higher-order statistics. To demonstrate the prominent effects of higher-order
distribution characteristics on distinguishing different nodes, we present an
example from Facebook Social Network (Traud et al., 2012) in Fig. 1.
Considering the graph composed of student/teacher nodes with class year as
attributes, Fig. 1 (b) shows the mean and variance of their neighbors’
features. And we cannot distinguish the label of node A and node B by the mean
of their neighbors’ features. However, we can succeed to distinguish the two
nodes if we consider the variance of their neighbors’ features. And it’s
intuitional because the teachers are more likely to know students from
different class years, thus the variance of teacher’s neighbors is larger. And
this motivates us to consider higher-order characteristics of neighbor’s
feature distribution for graph data, and we further demonstrate the critical
roles of higher-order characteristics by analysis on real data in Section 3.

Figure 1. An example in Facebook social network. (a) is an example of Facebook
college social network, where the student node A and the teacher node B have
different neighbors set. (b) shows the feature distributions of neighbors of A
and B by their mean and variance, where the mean information fails to
distinguish the two nodes.

A toy example of Facebook social network.

Inspired by the method of moments (Bowman and Shenton, 2004), which is used to
represent unknown probability distributions, we introduce methods of moments
to design a novel aggregator for GNNs. Specifically, we propose to use adapted
multi-order moments to represent the neighbor’s feature distribution and
further embed the distribution signals into the message passing mechanism of
GNNs. Then we design a novel Mix-Moment Graph Neural Network (MM-GNN) model,
which includes two major components: Multi-order Moments Embedding (MME) and
Attention-based Moment Adaptor (AMA). MM-GNN first computes multi-order
moments for neighbors of each node as their signatures. Based on the
observations that different order of moments may have various impact for each
node (Sec 3), an Element-wise Attention-based Moment Adaptor is designed to
select important moment signatures for specific nodes adaptively. As it is
impossible to compute infinite-orders of moments to represent the distribution
perfectly, we only compute finite-order moments in parallel at each graph
convolution year. Finally, we conduct extensive experiments to demonstrate the
superiority of MM-GNN over other state-of-the-art methods on 15 real-world
graphs, including social networks, citation networks, webpage networks and
image networks. The improved performance on all types of graphs demonstrate
the effectiveness of our method.

The main contributions of this paper are summarized as follows:

  1. (1)

We demonstrate that higher-order statistics are needed to model neighbor’s
feature distribution on graphs. And we introduce the method of moments to GNNs
from the perspective of modeling neighborhood feature distributions.

  2. (2)

We design a novel MM-GNN model with an aggregator based on the multi-order
moments and further design an Element-wise Attention-based Moment Adaptor to
adaptively adjust the weights of multi-order moments for each node.

  3. (3)

We conduct extensive experiments on 15 real-world graphs (i.e., social
networks, citation networks, webpage networks) to evaluate our method. The
results show that MM-GNN gains consistent improvements on each type of graphs.

##  2\. Preliminaries

In this section, we give the definitions of some important terminologies and
concepts appearing in this paper.

###  2.1. Graph Neural Networks (GNNs)

GNNs aim to learn representation for nodes on the graph. Let $G=(V,E)$ denotes
a graph, where $V=\\{v_{1},\cdots,v_{n}\\}$ is the node set, $n=|V|$ is the
number of nodes. Let $X\in\mathbb{R}^{n\times d}$ denote the feature matrix
and the $l$-th row of $X$ denoted as $x_{i}$ is the $d$-dimensional feature
vector of node $v_{i}$. $E$ is the edge set in $G$. Typically, GNN-based model
follows a neighborhood aggregation framework, where node representation is
updated by aggregating information of its first or multi-order neighboring
nodes. Let $h_{i}^{(l)}$ denotes the output vector of node $v_{i}$ at the
$l$-th hidden layer and let $h_{i}^{(0)}=x_{i}$. The update function of $l$-th
hidden layer is:

(1)
$h_{i}^{(l)}=\text{COMBINE}\left(h_{i}^{(l-1)},\text{AGGREGATE}\big{(}\\{h_{j}^{(l-1)}|v_{j}\in\mathcal{N}(v_{i})\\}\big{)}\right),$

where $\mathcal{N}(v_{i})$ is the set of neighbors of $v_{i}$. The AGGREGATE
function is aimed to gather information from neighbors and the goal of COMBINE
function is to fuse the information from neighbors and the central node. For
graph-level tasks, an additional READOUT function is required.

###  2.2. Method of Moments

####  2.2.1. Definitions about moments:

We give some basic definitions of statistical moments (Papoulis and Pillai,
2002; Bowman and Shenton, 2004) in this section. First, the definition of
moments is as follows:

###### Definition 2.1 (Moments in Mathematics).

For a random variable $Z$, its $k$-th order origin moment is denoted as
$\mathbb{E}(Z^{k})$ , where
$\mathbb{E}(Z^{k})=\int^{\infty}_{-\infty}z^{k}\cdot
f(z;\theta_{1},\cdots,\theta_{k})dz$. The $k$-th order central moment of $Z$
is
${\small\mathbb{E}\left(\big{(}Z-\mu_{z}\big{)}^{k}\right)}=\int^{\infty}_{-\infty}(z-\mu_{z})^{k}\cdot
f(z;\theta_{1},\cdots,\theta_{k})dz$. The $k$-th order standardized moment of
$Z$ is
${\small\mathbb{E}\left(\frac{\big{(}Z-\mu_{z}\big{)}^{k}}{\sigma_{z}^{k}}\right)}=\int^{\infty}_{-\infty}\frac{(z-\mu_{z})^{k}}{\sigma_{z}^{k}}\cdot
f(z;\theta_{1},\cdots,\theta_{k})dz$, where $\sigma_{z}$ is the standard
deviation of random variable Z.

Then, moments can be derived from the Moment-Generating Function (MGF), which
is an alternative specification of its probability distribution. The MGF
(Papoulis and Pillai, 2002) is defined as follows:

###### Definition 2.2 (Moment-Generating Function (MGF)).

Let $X$ be a random variable with cumulative distribution function $F_{X}$.
The Moment Generating Function (MGF) of $X$ (or $F_{X}$), denoted by
$M_{X}(t)$, is
$M_{X}(t)=\mathbb{E}(e^{tX})=\int_{-\infty}^{\infty}e^{tx}dF_{X}(x)$

An essential property is that MGF can uniquely determine the distribution. The
Taylor Expansion of MGF is as follows:

(2)
$M_{X}(t)=\mathbb{E}(e^{tX})=\sum_{k=0}^{\infty}\frac{t^{k}\mathbb{E}(X^{k})}{k!}=\sum_{k=0}^{\infty}\frac{t^{k}m_{k}}{k!}$

The $k$-th term of Eq. 2 is exactly the $k$-th order origin moment.

####  2.2.2. Relationship between moments and common statistics:

In order to intuitively explain the physical meaning of multi-order moments,
we give the relationship between moments and some common distribution
statistics of data. The common distribution statistics can be viewed as
special cases of moments. For a random variable $Z$, the expectation of $Z$ is
$\mathbb{E}(Z)$, which is the 1-st order origin moment. The variance of $Z$ is
$\mathbb{E}\left(\big{(}Z-\mu_{z}\big{)}^{2}\right)$, which is the 2-nd order
central moment. The skewness of $Z$ is
$\mathbb{E}\left(\frac{\big{(}Z-\mu_{z}\big{)}^{3}}{\sigma_{z}^{3}}\right)$,
which is the 3-rd order standardized moment.

##  3\. Findings: neighbor distribution benefits node classification

In this section, We verify the importance of neighbor’s feature distribution
through comprehensive data analysis. We first select three common statistics
(i.e., Mean, Variance, Skewness, see Sec. 2.2.2) to represent distribution of
neighbor’s feature. Then we calculate the selected statistics on the graphs of
Facebook social network (Traud et al., 2012), where the feature of nodes on
the graphs have practical meanings (e.g., gender, major, year) and easy to
analyze the physical meanings of their statistics. More details of the
datasets are presented in Sec. 6.1. And we further calculate the selected
statistics on the neighbors’ feature for all nodes and then take an average on
all feature dimensions. Finally, we choose two metrics (i.e. Fisher
Discrimination Index, Mutual Information) to measure the correlation between
statistics of neighbor distribution and the node label, and one statistic is
more important for node classification if it has stronger correlations with
the node label.

###  3.1. Data Analysis on Correlation between Neighbor Distribution and Label

(a) Fisher Discrimination Index. (b) Mutual Information.

Figure 2. Correlation analysis between the node label and neighbor’s feature
distribution represented by statistics (i.e., Mean, Variance, Skewness) on the
Facebook Social Networks.

####  3.1.1. Analysis by Fisher Discrimination Index

Inspired by Fisher’s Linear Discrimination (Fisher, 1936), we also use the
Fisher discrimination index (Ratio of inter-class variance to intra-class
variance) to measure the node classification capability distinguished by
different statistics of neighbors’ features distribution (e.g., Mean,
Variance):

$\small
S_{Fisher}=\frac{\sigma^{2}_{between}}{\sigma^{2}_{within}}=\sum_{\mathcal{C}_{i},\mathcal{C}_{j}\in\mathcal{C}}\frac{\left(\mu_{i}-\mu_{j}\right)^{2}}{\sigma^{2}_{i}+\sigma^{2}_{j}},$

where $S_{Fisher}$ denotes the Fisher Discrimination Index,
$\sigma^{2}_{between}$ and $\sigma^{2}_{within}$ denote the inter-class and
intra-class variation. $\mathcal{C}$ is the set of classes, $\mathcal{C}_{i}$
and $\mathcal{C}_{j}$ are two different categories, $\mu_{i}$ and
$\sigma^{2}_{i}$ are respectively the mean and variance of samples in the
class $\mathcal{C}_{i}$. As shown in Fig. 2 (a). We observe that, for
different feature dimensions, different statistic of neighbor distribution
contributes variously to distinguish node labels. For example, we observe that
the skewness statistic has the best discriminative ability for the feature
dimension ”year”, while the variance statistic has best discriminative ability
for the feature dimension ”major”.

####  3.1.2. Analysis by Mutual Information

Apart from Fisher Discrimination Index, we also analyze the Mutual Information
between distribution statistics and the node label to measure theirs
correlations. As shown in Fig. 2 (b), the results of mutual information also
demonstrate that different statistics have different strength of correlation
with the class label. In summary, all these results show that different order
of moments contribute quite differently to node classification. However, any
single statistic (e.g., mean, max, sum) cannot reflect complete neighborhood
feature distribution, we thus need to consider multi-order statistics such as
variance, skewness and other higher-order distribution statistics.

Figure 3. An overview of MM-GNN including Multi-order Moment Embedding (MME)
and Attention-based Moment Adaptor (AMA) modules. $H^{(l)}_{N(v_{i})}$ and
$P^{(l)}_{N(v_{i})}$ denotes the features and distributions (represented by
moments) of neighbors at the $l$-th layer. $h_{i}^{(l)}$ is the node
representation of node $v_{i}$ at the $l$-th layer. $M_{K}^{(l+1)}$ denotes
the $K$-th order moment of neighbors at the $(l+1)$-th layer. $Q_{i}$ and
$K_{i}^{k}$ denotes query and key vectors at the $l$-th layer.
$W_{query}^{(l+1)}$, $W_{key}^{(l+1)}$ and $W_{query}^{(l+1)}$ are learnable
transformation matrix at the $(l+1)$-th layer. $a_{k}^{(l+1)}$ is the learned
attention matrix for $k$-th order moment at the $(l+1)$-th layer.

Figure of model structure

##  4\. Mix-Moment Graph Neural Network

In this section, we introduce our Mix-Moment Graph Neural Network (MM-GNN)
model. An overview of MM-GNN is given in the Fig. 3. Compared with existing
GNNs, our MM-GNN introduces the Method of Moments to better represent
neighborhood feature distribution when updating node representations.

###  4.1. Modeling Neighbors’ Feature Distribution via Method of Moments

Based on the foundations of method of moments given in Sec. 2.2, we introduce
multi-order moments to represent the feature distribution of neighboring nodes
on the graph. And based on the analysis in Sec. 3, we observe that, for the
node classification task, the feature distribution of neighbors for the nodes
belonging to the same class follow the similar distributions, and vice versa.
Unlike other parametric methods, the outstanding power of estimating
distribution with method of moments is that we do not need to assume an
explicit form of the distribution function. In this paper, we propose to use
finite-order moments to represent the neighbors’ feature distribution for a
graph, and we further present the analysis on the theoretical and empirical
deviation of representing the neighbors’ feature distribution for graph data
with finite-order moments in Sec 5.2.

###  4.2. Multi-order Moments Embedding (MME)

As aforementioned, we model neighbors’ feature distribution by introducing the
Method of Moments to GNN. The Method of Moments illustrated in Section 2.2 is
a statistical method that use sample moments to estimate the distribution
function. Therefore, we use multi-order moments to estimate the feature
distribution of neighboring nodes. Moreover, moments of different nodes can be
computed in parallel at the aggregation step, and we select both origin moment
and central moment to represent neighborhood feature distribution. Here we
present the form of $k$-th order origin moment for $(l+1)$-th layer of MM-GNN:

(3)
$M_{k}^{(l+1)}(v_{i})=\mathbb{E}\big{(}(h_{j}^{(l)})^{k}\big{)}=\left(\frac{1}{|\mathcal{N}(v_{i})|}\sum_{j\in\mathcal{N}(v_{i})}((h_{j}^{(l)})^{k}\right)^{\frac{1}{k}}\cdot
W_{k}^{M}$

where $W_{k}^{M}\in\mathbb{R}^{D_{in}\times D_{hidden}}$ is the transformation
matrix for $k$-th order moment signature, $h_{j}^{(l)}\in\mathbb{R}^{D_{in}}$
is the node representation of $v_{j}$ output by the $l$-th layer, and
$h_{j}^{0}$ is the raw node feature $x_{j}$. Based on the Eq. 3, we can
observe that the $1$-st origin moment is exactly the expectation, which is
equivalent to the MEAN aggregation strategy. Thus previous GNNs with MEAN
aggregator are special cases of MM-GNN using $1$-st order origin moment.

Similarly, we can give the form of $k$-th order central moment for $(l+1)$-th
layer of GNNs:

(4)
$M_{k}^{(l+1)}(v_{i})=\mathbb{E}\big{(}(h_{j}^{(l)}-\mu_{i}^{(l)})^{k}\big{)}=\left(\frac{1}{|\mathcal{N}(v_{i})|}\sum_{j\in\mathcal{N}(v_{i})}(x_{j}-\mu_{i}^{(l)})^{k}\right)^{\frac{1}{k}}\cdot
W_{k}^{M}$

, where

(5)
$\mu_{i}^{(l)}=\frac{1}{|\mathcal{N}(v_{i})|}\sum_{j\in\mathcal{N}(v_{i})}h_{j}^{(l)}$

To keep the same order of magnitude, we normalize moments of different orders
by computing the $\frac{1}{k}$ power of $k$-th order moment vectors in the
model. Besides, we set a hyper-parameter for flexible model configuration of
using the central moment or origin moment.

###  4.3. Attention-based Moment Adaptor (AMA)

With the Multi-order Moments Embedding (MME) module, we get signatures
representing moments of different orders for each node. Then, we need to fuse
these signatures into one final representation for each node. As illustrated
in Section 3, different statistics (corresponding to moments) have different
effects on node classification, we design an Element-wise Attention-based
Moment Adaptor (AMA) that can adaptively select the important order of moments
for different nodes and feature dimensions. For each MM-GNN layer, after we
get the multi-order moment signatures with the MME module, we use the
transformed input node representation (the output of the previous layer) as
the query vector for each node, and then we use different moment-signatures
(output of the MME module of current layer) as key vectors. Then we compute
attentions as weights of different moments for each individual nodes. The
attention of $k$-th order moment for $v_{i}$ at $(l+1)$-th layer, denoted by
$a_{k}^{(l+1)}(v_{i})\in\mathbb{R}^{D_{hidden}}$, is:

(6) $a_{k}^{(l+1)}(v_{i})=\sigma\left([h_{i}^{(l)}\cdot
W_{query}^{(l+1)}||M_{k}^{(l+1)}(v_{i})\cdot W_{key}^{(l+1)}]\cdot
W_{a}^{(l+1)}\right)$

where $W_{a}^{(l)}\in\mathbb{R}^{2D_{hidden}\times D_{hidden}}$ is the
learnable attention matrix in $l$-th layer,
$W_{query}^{(l)},W_{key}^{(l)}\in\mathbb{R}^{D_{hidden}\times D_{hidden}}$ are
the transformation matrix for the query and key vectors in $l$-th layer. Then
different dimensions of $k$-th order moment signature for node $v_{i}$ use
different attention weights and therefore implement an element-wise attention
mechanism (attention of each order moment is a $D_{hidden}$-dimensional vector
rather than a scalar). Then we can get the output representation of node
$v_{i}$ at $(l+1)$-th layer:

(7) $h_{i}^{(l+1)}=\sum_{k=1}^{K}a_{k}^{(l+1)}(v_{i})\odot
M_{k}^{(l+1)}(v_{i}),\ a_{k}^{(l+1)}(v_{i})\in\mathbb{R}^{1\times D_{hidden}}$

where $\odot$ is the element-wise product and $K$ is the hyper-parameter for
the largest ordinal of moments used in MM-GNN. Besides, we add residual
connections across layers to further boost the performance of MM-GNN by
aggregating neighborhood feature distribution of multi-hop neighbors. We
finally use the output of the last GNN layer as the node representation and
further compute Cross Entropy Loss for node classification task and back
propagate to optimize the model parameters.

Table 1. Information of Facebook social networks. We present the number of nodes, number of edges for each social graph in this table. Each node has 6-dimensional attributes and belongs to 6 possible classes. Dataset | Northeastern | Caltech | UF | Hamilton | Howard | Simmons | GWU | UGA | Tulane
---|---|---|---|---|---|---|---|---|---
#Nodes | 13882 | 769 | 35123 | 2314 | 1909 | 1518 | 12193 | 24389 | 7752
#Edges | 763868 | 33312 | 2931320 | 192788 | 409700 | 65976 | 939056 | 2348114 | 567836
Table 2. Information of citation/webpage/image networks Dataset | #Nodes | #Edges | #Features | #Classes
---|---|---|---|---
Cora | 2,708 | 10,556 | 1,433 | 7
CiteSeer | 3,327 | 9,104 | 3,703 | 6
PubMed | 19,717 | 88,648 | 500 | 3
Chameleon | 2277 | 36101 | 2325 | 5
Squirrel | 5201 | 217073 | 2089 | 5
Flickr | 89,250 | 899,756 | 500 | 7

##  5\. Deep Analysis of MM-GNN

In this section, we first theoretically analyze the limitations of existing
GNNs using single statistic, which cannot retain the full information of
neighbor distribution. Then we analyze the deviation of estimating neighbors’
feature distribution with finite-order moments. Finally, we analyze the time
complexity of MM-GNN.

###  5.1. Limitations of Existing GNNs

In this sections, we prove the existing GNNs using single statistic (e.g. GCN,
GraphSAGE) have severe limitations when representing neighbor’s feature
distribution, which reduces the generalization ability of the model.
Complexity measure (Neyshabur et al., 2017) is the current mainstream method
to measure the generalization ability of the model. And we choose Consistency
of Representations (Natekar and Sharma, 2020) (champion of the NIPS2020
Competition on generalization measure). And the complexity measure of the
model, denoted as $\Gamma$, is:

(8) $\Gamma=\frac{1}{|\mathcal{C}|}\sum_{i=1}^{|\mathcal{C}|}\max_{i\neq
j}\frac{\mathcal{S}_{i}+\mathcal{S}_{j}}{\mathcal{M}_{i,j}},$

where $\mathcal{C}$ is the set of classes,
$\mathcal{C}_{i},\mathcal{C}_{j}\in\mathcal{C}$ are two different classes,
$\mathcal{S}_{i}=\left(\mathbb{E}_{v_{k}\sim\mathcal{C}_{i}}(|h_{k}-\mu_{\mathcal{C}_{i}}|^{p})\right)^{\frac{1}{p}}$
is the intra-class variance of class $\mathcal{C}_{i}$,
$\mathcal{M}_{i,j}=||\mu_{\mathcal{C}_{i}}-\mu_{\mathcal{C}_{j}}||_{p}$ is the
inter-class variance between $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$, $h_{k}$
is the representations of nodes $v_{k}$ learned by the model and
$\mu_{\mathcal{C}_{i}}=\mathbb{E}_{v_{k}\sim\mathcal{C}_{i}}(h_{k})$. And
higher complexity measure means lower generalization ability. Then we
demonstrate that GNNs using single statistic have weak generalization ability.
Take GCN (Kipf and Welling, 2016) as an example, we give the following
theorem:

###### Theorem 5.1.

Given a graph $G(V,E)$, we denote the nodes belonging to class
$\mathcal{C}_{i}$ as $\\{v_{i}|v_{i}\in\mathcal{C}_{i}\\}$, and assume that
the neighbor’s feature distribution of nodes in $\mathcal{C}_{i}$ follows the
same i.i.d Gaussian distribution $\mathrm{N}(\bm{\mu}_{i},\bm{\Sigma}_{i})$.
Given any two classes $\mathcal{C}_{i}\neq\mathcal{C}_{j}$, if
$\bm{\mu}_{i}=\bm{\mu}_{j}$ and $\bm{\Sigma}_{i}\neq\bm{\Sigma}_{j}$ , then
the complexity measure $\Gamma$ (Eq.8) of GCN with mean aggregation reaches
$+\infty$, and the GCN thus lose the generalization ability.

###### Proof.

The update function of GCN with mean aggregator is:

$h^{(l+1)}_{k}=\sum_{j\in\mathcal{N}(v_{k})}\frac{1}{d_{k}}\cdot
h_{j}^{(l)}\cdot
W^{(l)}=\left(\frac{1}{d_{k}}\cdot\sum_{j\in\mathcal{N}(v_{k})}h_{j}^{(l)}\right)\cdot
W^{(l)}$

To simplify the analysis, we omit the activation function. And we can get the
expectation of $h^{(l+1)}_{k}$ for nodes in $\mathcal{C}_{i}$ and
$\mathcal{C}_{j}$:

(9) $\mathbb{E}_{v_{k}\sim\mathcal{C}_{i}}(h_{k}^{(l+1)})=\bm{\mu}_{i}\cdot
W^{(l)}=\mathbb{E}_{v_{k}\sim\mathcal{C}_{j}}(h_{k}^{(l+1)})=\bm{\mu}_{j}\cdot
W^{(l)},$

which means the expectation of the learned representations of nodes in
$\mathcal{C}_{i}$ equals to that of nodes in $\mathcal{C}_{j}$. Then according
to Eq. 9:

$\mathcal{M}_{i,j}=||\mathbb{E}_{v_{k}\sim\mathcal{C}_{i}}(h_{k})-\mathbb{E}_{v_{k}\sim\mathcal{C}_{j}}(h_{k})||_{p}=\bm{0}$

Considering that $\mathcal{S}_{i}+\mathcal{S}_{i}>0$, the complexity measure
$\Gamma$ in Eq. 8 equals to $+\infty$. Then GCN loses generalization ability,
and this rule also applies to any other GNNs with mean aggregator. ∎

And we can intuitively generalize Theorem 5.1 to all GNNs using single
statistics. For example, GraphSAGE (Hamilton et al., 2017) use the max/mean
aggregator, GIN (Xu et al., 2018a) uses the sum aggregator. And most state-of-
the-art GNNs (e.g., GCNII (Chen et al., 2020), DAGNN (Liu et al., 2020)) still
use single statistic to represent the neighbor’s feature distribution , and
therefore have the same limitations when neighbor’s feature distribution has
more than one degree of freedom.

###  5.2. Deviation Analysis on Finite-order Moments.

We give the theoretical guarantee of estimating neighborhood feature
distributions on the graph with finite-order moments by a theorem and give its
proof as follows:

###### Theorem 5.2.

Given a graph $G$, we denote the deviation of using top-$K$ order moments to
represent the neighbor’s feature distribution as $e_{K}$, and $e_{K}$ wil not
larger than an upper bound: $e_{K}\leq\frac{(\epsilon\cdot c)^{k+1}}{(k+1)!}$,
where $c\in\mathbb{R}^{+}$ is a constant and $\epsilon\in\mathbb{R}$ is a
small number around 0. And larger $K$ leads to smaller estimation deviation
(convergence to 0).

###### Proof.

According to the Sec. 2.2, the $k$-th order moment can be derived from the
$k$-th term of the MGF’s Taylor Expansion (see eq. 2). Therefore, the top-$K$
order moments can be viewed as an estimation of the $K$-th order Taylor
expansion for MGF of neighbor’s feature distribution. In real cases, the
feature of samples is usually bounded. We assume that for $\forall x_{i,j}\in
X$, we have $|x_{i,j}|\leq c$, where $c\in\mathbb{R}^{+}$ is a constant.
Besides, sine the Eq. 2 is the Taylor expansion around $t=0$, then $t$ is a
small real number $t\leq\epsilon$, where $\epsilon\in\mathbb{R}$ is a small
number around 0. Then the deviation of estimating neighbor’s feature
distribution by $k$-th order Taylor expansion of MGF can be approximated by
its Lagrange Remainder:

(10) $R_{k}(x)=\frac{t^{k+1}\mathbb{E}(X^{k+1})}{(k+1)!}\leq\frac{(t\cdot
c)^{k+1}}{(k+1)!}\leq\frac{(\epsilon\cdot c)^{k+1}}{(k+1)!}$

Besides, in real scenarios, we usually normalize the data, so the upper bound
of the estimation error $\frac{(t\cdot c)^{k+1}}{(k+1)!}$ is further reduced.
And it’s obvious that the upper bound of estimation error converges to zero as
the moment ordinal $k$ increases. ∎

###  5.3. Time Complexity Analysis

Here we analyze the time complexity of training MM-GNN. Let $D_{in}$ denote
the dimension of input node feature, and $D_{hidden}$ denote the hidden
dimension of node representation at each layer. $|V|$ and $|E|$ denotes the
number of nodes on the graph. For each layer of MM-GNN, the aggregation and
transformation of graph convolution cost $O\left((|V|+|E|)\cdot
D_{in}D_{hidden}\right)$. Let $k$ denote the max ordinal of used moments, and
we used $k\leq 3$ in our experiments and thus can be seen as a small constant
integer. The calculation of query vector costs $O(|V|D_{hidden}^{2})$, and the
calculation of all key vectors costs $O\left(k\cdot(|V|+|E|)\cdot
D_{hidden}^{2}\right)$. The Attention-based Moment Adaptor costs
$O(k|V|D_{hidden}^{2})$. Considering that $D_{hidden}\leq D_{in}$ and $k$ is a
constant integer, the overall time complexity of MM-GNN layer is:

$\begin{aligned} T&=O\left((k+1)|V|D_{hidden}^{2}+(k+1)\cdot(|V|+|E|)\cdot
D_{in}D_{hidden}\right)\\\ &=O\left((|V|+|E|)\cdot
D_{in}D_{hidden}\right)\end{aligned},$

which equals to the complexity of other mainstream GNNs.

##  6\. Experiments

###  6.1. Datasets

The proposed MM-GNN is evaluated on 15 real-world graphs, including 9 social
graphs from Facebook social networks, 3 citation networks (Cora, CiteSeer,
PubMed), 2 webpage networks (Chameleon, Squirrel) and one image network
(Flickr). The detailed information of theses datasets are presented in Table 1
and Table 2.

Social network datasets:  Facebook100 (Traud et al., 2012) provides 100
college graphs, each graph describes the social relationship in a university.
Each graph has categorical node attributes with practical meaning (e.g.,
gender, major, class year.). Moreover, nodes in each dataset belong to six
different classes (a student/teacher status flag).

Citation network datasets:  We choose three public citation network datasets
(Cora, PubMed, CiteSeer) to evaluate our model. The nodes on such graphs
represent papers and edges represent citation relationships between papers.

Webpage network datasets:  Two webpage network datasets (Chameleon, Squirrel)
are used for model evaluation. Nodes of these datasets represent articles and
edges are mutual links between them.

Image network datasets:  We also choose one image network (Flickr) for
evaluation. The nodes represent images and edges means two images share the
same metadata (e.g., location, gallery, tags).

Table 3. Performance (accuracy) on node classification task ($60\%$ training set ratio) on Facebook social network datasets. Model | Northeastern | Caltech | UF | Hamilton | Howard | Simmons | GWU | UGA | Tulane
---|---|---|---|---|---|---|---|---|---
GCN | $90.58\pm 0.28$ | $88.47\pm 1.91$ | $83.94\pm 0.61$ | $92.26\pm 0.35$ | $91.21\pm 0.69$ | $88.74\pm 0.61$ | $86.66\pm 0.48$ | $85.73\pm 0.51$ | $87.93\pm 0.97$
GAT | $88.69\pm 0.30$ | $81.17\pm 2.15$ | $81.68\pm 0.59$ | $91.43\pm 1.25$ | $89.35\pm 0.44$ | $89.14\pm 0.57$ | $84.27\pm 0.88$ | $82.95\pm 0.75$ | $84.45\pm 1.40$
GraphSAGE | $91.77\pm 0.86$ | $92.14\pm 1.38$ | $86.21\pm 0.56$ | $94.75\pm 0.35$ | $93.13\pm 0.93$ | $92.31\pm 0.92$ | $89.43\pm 0.78$ | $88.11\pm 0.39$ | $89.87\pm 0.47$
GIN | $89.62\pm 1.10$ | $78.35\pm 4.24$ | $85.40\pm 0.71$ | $85.47\pm 6.20$ | $89.37\pm 2.05$ | $88.37\pm 1.18$ | $84.73\pm 2.80$ | $87.68\pm 0.40$ | $80.71\pm 4.78$
APPNP | $91.11\pm 0.27$ | $90.76\pm 2.38$ | $83.07\pm 0.54$ | $93.99\pm 0.47$ | $91.88\pm 1.22$ | $90.31\pm 1.20$ | $90.43\pm 0.38$ | $86.31\pm 0.41$ | $88.52\pm 0.44$
JKNet | $92.32\pm 0.39$ | $92.34\pm 1.09$ | $85.43\pm 0.71$ | $94.72\pm 0.76$ | $93.11\pm 0.78$ | $91.81\pm 0.71$ | $89.25\pm 0.65$ | $87.68\pm 0.40$ | $89.04\pm 0.43$
DAGNN | $92.43\pm 0.30$ | $91.15\pm 2.67$ | $86.71\pm 0.51$ | $94.75\pm 0.75$ | $93.44\pm 0.68$ | $92.49\pm 0.96$ | $89.83\pm 0.42$ | $88.20\pm 0.29$ | $89.29\pm 0.50$
MM-GNN | $93.16\pm 0.32$ | $93.28\pm 1.94$ | $88.20\pm 0.45$ | $95.46\pm 0.53$ | $94.08\pm 0.95$ | $93.61\pm 0.45$ | $90.78\pm 0.45$ | $89.44\pm 0.30$ | $90.22\pm 0.33$
Table 4. Performance (accuracy) on node classification task on citation/webpage/image network datasets Dataset | Cora | CiteSeer | PubMed | Chameleon | Squirrel | Flickr
---|---|---|---|---|---|---
GCN | $81.13\pm 0.26$ | $70.08\pm 0.43$ | $78.62\pm 0.63$ | $37.68\pm 3.06$ | $26.39\pm 0.88$ | $50.31\pm 0.30$
GAT | $81.36\pm 0.96$ | $70.93\pm 0.70$ | $78.19\pm 0.78$ | $44.34\pm 1.42$ | $29.82\pm 0.98$ | $50.59\pm 0.26$
GraphSAGE | $80.43\pm 0.39$ | $69.56\pm 0.29$ | $77.47\pm 0.76$ | $47.06\pm 1.88$ | $35.62\pm 1.21$ | $50.21\pm 0.31$
GIN | $76.77\pm 0.88$ | $67.67\pm 0.34$ | $77.03\pm 0.53$ | $32.18\pm 1.98$ | $25.08\pm 1.36$ | $46.65\pm 0.33$
APPNP | $82.57\pm 0.10$ | $70.53\pm 0.25$ | $78.33\pm 0.39$ | $40.44\pm 2.02$ | $29.20\pm 1.45$ | $49.21\pm 0.37$
JKNet | $79.43\pm 0.53$ | $70.36\pm 0.33$ | $77.80\pm 0.82$ | $42.43\pm 1.76$ | $35.52\pm 1.06$ | $49.83\pm 0.58$
DAGNN | $83.93\pm 0.63$ | $72.56\pm 1.76$ | $80.50\pm 0.75$ | $47.29\pm 0.62$ | $36.23\pm 1.58$ | $50.78\pm 0.62$
GPR-GNN | $81.85\pm 0.38$ | $70.64\pm 0.65$ | $79.78\pm 0.43$ | $60.83\pm 1.37$ | $49.21\pm 1.15$ | $50.01\pm 0.53$
MM-GNN | $84.21\pm 0.56$ | $73.03\pm 0.58$ | $\underline{80.26\pm 0.69}$ | $63.32\pm 1.31$ | $51.38\pm 1.73$ | $\bm{51.73\pm 0.35}$

###  6.2. Baselines

We compare our proposed MM-GNN model against eight widely-used GNN models,
including four baseline GNNs ( GCN (Kipf and Welling, 2016), GAT (Veličković
et al., 2017), GraphSAGE (Hamilton et al., 2017), GIN (Xu et al., 2018a)) and
four state-of-the-art GNNs (APPNP (Klicpera et al., 2018), JKNet (Xu et al.,
2018b), DAGNN (Liu et al., 2020), GPRGNN (Chien et al., 2020)). Considering
that the two webpage networks (Chameleon, Squirrel) are heterophily graph
datasets unlike other datasets used in this paper, we also choose GPR-GNN
(Chien et al., 2020), which designed for graphs with heterophily as one of the
baseline models.

###  6.3. Experimental Setup

We evaluate MM-GNN on the semi-supervised node classification task compared
with state-of-the-art methods. For each dataset, we all run our model and
other baseline models for 10 times and compute the mean and standard deviation
of the results. For Facebook social networks, we randomly generate 3 different
data splits with an average train/val/test split ratio of 60%/20%/20% in the
Table 3. And we further generate data splits with lower training ratio
(40%,20%,10%) and conduct experiments to validate the robustness of our model
in the Table 6. For citation networks (Cora, Citeseer, and Pubmed), we use the
public split recommended by (Kipf and Welling, 2016) with fixed 20 nodes per
class for training, 500 nodes for validation, and 1000 nodes for testing. For
webpage networks (Chameleon, Squirrel) (Pei et al., 2020) and image network
(Flickr (Zeng et al., 2019)), we all use the public splits recommended in the
original papers. The max moment ordinal $K$ for MM-GNN is searched in {1, 2,
3} in our experiments. To keep fairness, we search the shared hyper-parameters
for all models in the same searching spaces. Specifically, we train
200/400/600 epochs for each model. we search the hidden units in {16, 32, 64,
128, 256} and search the number of GNN layers in {2, 6, 10, 16}. We search
learning rate from $1e-4$ to $1e-1$ and the weight decay from $1e-4$ to
$1e-2$. For each experiment, we use the Adam optimizer (Kingma and Ba, 2014)
to train each GNN model on Nvidia Tesla V100 GPU.

###  6.4. Main Results

We compare the performance of MM-GNN to the state-of-the-art methods in Table
3 and Table 4. Compared with all baselines, MM-GNN generally achieves the best
performance on all datasets. For social networks, we observe from Table 3 that
our MM-GNN model gains significant improvement on all datasets. Besides, MM-
GNN also gain significant improvements on other public real-world datasets,
including citation networks, webpage networks and image networks. And the
results demonstrate that MM-GNN gains superior performance on the graphs of
various types.

###  6.5. Ablation Study

To estimate the effectiveness of each component in MM-GNN, we conduct
comprehensive ablation studies by removing certain component at a time. The
results are shown in the Table 5.

Table 5. Performance (accuracy) on node classification task compared with single-moment GNN models. M-1, M-2, M-3 represent the single-moment GNN models, where M-$k$ denotes the model whose aggregator only uses the $k$-th order moment of neighbors. Ensemble denotes the mean ensemble of the three single-moment models.  Model | M-1 | M-2 | M-3 | Ensemble | MM-GNN
---|---|---|---|---|---
Fusion mode | — | — | — | Mean | MLP | Attention
Howard | 93.15 | 93.61 | 93.83 | 93.70 | 94.00 | 94.37
Simmons | 90.90 | 92.36 | 92.45 | 91.85 | 92.99 | 93.66
GWU | 88.57 | 89.79 | 90.44 | 89.69 | 90.48 | 91.11
Cora | 82.31 | 81.70 | 80.07 | 81.98 | 82.58 | 84.21
CiteSeer | 71.31 | 70.70 | 71.53 | 71.28 | 72.33 | 73.03
PubMed | 79.23 | 78.01 | 78.27 | 78.93 | 79.35 | 80.26
Chameleon | 48.17 | 48.09 | 48.24 | 49.76 | 50.05 | 51.93
Squirrel | 35.58 | 36.29 | 35.97 | 36.09 | 36.35 | 38.92
Flickr | 50.37 | 51.09 | 50.88 | 51.15 | 51.37 | 51.73

####  6.5.1. Effects of MME

To validate the effects of Multi-order Moments Embedding (MME) module, we
design several model variants as shown in Table 5. The variants in the first
three columns use only one moment (e.g. 1st-order moment) for neighborhood
aggregation. The Ensemble-Mean in the table is another variant that ensemble
all single-moment models with the average logits of all single-moment
variants. And we observe that each single-moment model has different
performance on different datasets. And the simple ensemble of them usually
have better performance. However, our full MM-GNN (the rightmost column)
always has the best performance on all datasets against all variants.

####  6.5.2. Effects of AMA

To validate the effects of the Attention-based Moment Adaptor (AMA) module, we
design a variant of MM-GNN by removing this module and use a simple MLP to
fuse the output of MME module. MM-GNN with Attention-based Moment Adaptor
gains the state-of-the-art results on all datasets compared with the simple
MLP fusion method, which demonstrates the superiority of the Attention-based
Moment Adaptor module.

Table 6. Performance (accuracy) on node classification task with different training set ratio on Facebook social networks.  Dataset | Northeastern | UF | GWU | UGA | Howard
---|---|---|---|---|---
Training | $40\%$ | $20\%$ | $10\%$ | $40\%$ | $20\%$ | $10\%$ | $40\%$ | $20\%$ | $10\%$ | $40\%$ | $20\%$ | $10\%$ | $40\%$ | $20\%$ | $10\%$
GCN | 90.34 | 90.38 | 90.35 | 83.97 | 83.76 | 83.62 | 85.76 | 85.65 | 85.67 | 85.5 | 85.46 | 85.42 | 91.17 | 90.95 | 90.99
GAT | 88.57 | 88.63 | 88.69 | 79.79 | 79.89 | 79.85 | 83.95 | 84.07 | 83.98 | 82.87 | 82.92 | 82.79 | 89.09 | 88.92 | 88.93
GraphSAGE | 91.91 | 92.22 | 92.18 | 86.26 | 86.34 | 86.51 | 89.31 | 89.45 | 89.15 | 88.21 | 88.09 | 87.77 | 93.21 | 93.29 | 93.17
GIN | 89.45 | 89.59 | 89.57 | 85.64 | 85.44 | 85.57 | 84.27 | 83.98 | 84.07 | 87.57 | 87.44 | 87.34 | 89.40 | 89.24 | 89.31
APPNP | 89.76 | 89.94 | 89.85 | 81.89 | 81.92 | 81.87 | 87.14 | 87.11 | 87.12 | 84.79 | 84.75 | 84.67 | 89.20 | 89.51 | 89.27
JKNet | 92.23 | 92.35 | 92.28 | 85.51 | 85.67 | 85.49 | 89.12 | 89.06 | 89.02 | 87.91 | 87.82 | 87.81 | 93.15 | 93.14 | 93.17
DAGNN | 92.49 | 92.36 | 92.30 | 86.86 | 86.46 | 85.70 | 89.73 | 89.31 | 89.11 | 88.35 | 87.79 | 87.62 | 93.31 | 93.33 | 93.27
MM-GNN | 93.08 | 92.96 | 92.73 | 88.34 | 87.87 | 87.38 | 90.63 | 90.50 | 90.14 | 89.52 | 89.38 | 88.91 | 94.03 | 93.88 | 93.76

(a) Social network (b) Citation network

(c) Webpage/Image network

Figure 4. Average attention value of different moments. The value in each
block means the averaged attentions of all nodes and all feature dimensions.
And the deeper color means larger attention values learned by MM-GNN.

###  6.6. Robustness Analysis

In this section, we evaluate the robustness of our MM-GNN by changing the
training set ratio and the key hyper-parameters.

####  6.6.1. Analysis on effects of training set ratio:

We conducted extensive experiments with different training set ratios (10%,
20%, 40%) on the Facebook social networks. And due to space limitations, we
only show results on five of these datasets. As shown in Table 6, our MM-GNN
model has the best performance on each split of the five graphs even when the
proportion of training set is small.

(a) Simmons, GWU, Cora (b) Chameleon, Squirrel, Flickr

Figure 5. Hyper-parameter analysis on max ordinal of moments $K$ used in MM-
GNN

####  6.6.2. Analysis on the hyper-parameter $K$:

As shown in the Fig. 5, the performance of MM-GNN is improved with larger $K$
(max ordinal of moments). And when $K$ is no less than 3, the performance of
our model is convergent, which matches the data complexity.

###  6.7. Visualizations

The Attention-based Moment Adaptor (AMA) improves the classification
performance and also provides interpretability for MM-GNN, which plays an
critical role in selecting important orders of moments for each node on the
graph. To further explore the interpretability of the learned attention, we
analyze their effects through attention visualization. We output the attention
value learned by the first layer of MM-GNN for each node and compute the
averaged attention of all nodes on all hidden dimensions. As shown in the Fig.
4, the average attention values for each moment on different graph dataset
present different patterns. MM-GNN output larger weights for higher-order
moments on the social networks, while output smaller weights for higher-order
moments on citation networks. Moreover, the magnitude of learned attentions
are aligned with the performance of the single-moment model shown in Table 5,
which further demonstrates the interpretability of MM-GNN.

##  7\. Related work

In recent years, Graph Neural Networks (GNNs) have been widely studied on
graph representation learning for its expressive power (Du et al., 2021a; Yao
et al., 2022; Wang et al., 2019; Zhu et al., 2020; Monti et al., 2018; Wang et
al., 2022; Bi et al., 2022b). Spectral GNNs (Defferrard et al., 2016; Boscaini
et al., 2015; Xu et al., 2019) were first proposed based on graph signal
processing theory. Then Kipf et al. proposed GCN (Kipf and Welling, 2016),
which simplifies the convolution operation with low pass filtering on graphs
and bridges the gap between spectral GNNs and spatial GNNs. GAT further
(Veličković et al., 2017) introduced an attention mechanism to GNNs to learn
the importance of different neighbors, instead of using equal weights for
neighbors. Afterwards, some researchers turned their attention to considering
long-range dependencies (Dwivedi et al., 2021; Ying et al., 2021; Abu-El-Haija
et al., 2019). DAGNN (Liu et al., 2020) is a deeper GCN model with decoupled
transformation and aggregation strategy .

However, most existing GNNs update node representations by aggregating feature
from neighbors with single statistic (e.g., mean, max, sum). Though such
lower-order statistics can also reflect certain characteristics of neighbor
distribution, they lose the full information of the neighbor’s feature
distribution. Even for the simple Gaussian distribution which has two degree
of freedom (mean, variance), and the information of 2nd-order statistic
variance cannot be retained by mean/max/sum (1st-order statistic). Some graph
filter based methods (Dehmamy et al., 2019; Corso et al., 2020) proposed to
model graph signals by integrating multiple types of aggregators, without
considering the relationship between different signals and neighborhood
distribution. Recently, some works (Du et al., 2021b; Ma et al., 2022; Luo et
al., 2022; Bi et al., 2022a) begin to study the modeling of neighbor
distribution on the graph, but they do not consider the information loss
caused by single lower-order statistics.

##  8\. Conclusion

In this paper, we embrace the critical roles of neighbor’s feature
distribution for modeling graph data. And we introduce the multi-order moments
to model neighbor’s feature distribution on graphs, fixing the severe
limitations of existing GNNs using only single statistic for neighborhood
aggregation. Based on the method of moments, we propose a novel Mix-Moment GNN
model with adaptive attention mechanism, which includes two components, the
Multi-order Moments Embedding (MME) and the Attention-based Moment Adaptor
(AMA). MM-GNN use multi-order moments to represent neighbor’s feature
distribution and select important moments for nodes adaptively. The consistent
improvement of MM-GNN over other state-of-the-art GNNs on 15 real-world graphs
demonstrate the superiority of our method.

## References

  * (1)
  * Abu-El-Haija et al. (2019) Sami Abu-El-Haija, Bryan Perozzi, Amol Kapoor, Nazanin Alipourfard, Kristina Lerman, Hrayr Harutyunyan, Greg Ver Steeg, and Aram Galstyan. 2019.  Mixhop: Higher-order graph convolutional architectures via sparsified neighborhood mixing. In _international conference on machine learning_. PMLR, 21–29.
  * Bai et al. (2019) Yunsheng Bai, Hao Ding, Song Bian, Ting Chen, Yizhou Sun, and Wei Wang. 2019\.  Simgnn: A neural network approach to fast graph similarity computation. In _Proceedings of the Twelfth ACM International Conference on Web Search and Data Mining_. 384–392.
  * Bi et al. (2022a) Wendong Bi, Lun Du, Qiang Fu, Yanlin Wang, Shi Han, and Dongmei Zhang. 2022a.  Make Heterophily Graphs Better Fit GNN: A Graph Rewiring Approach.  _arXiv preprint arXiv:2209.08264_ (2022).
  * Bi et al. (2022b) Wendong Bi, Bingbing Xu, Xiaoqian Sun, Zidong Wang, Huawei Shen, and Xueqi Cheng. 2022b.  Company-as-Tribe: Company Financial Risk Assessment on Tribe-Style Graph with Hierarchical Graph Neural Networks. In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_. 2712–2720.
  * Boscaini et al. (2015) Davide Boscaini, Jonathan Masci, Simone Melzi, Michael M Bronstein, Umberto Castellani, and Pierre Vandergheynst. 2015.  Learning class-specific descriptors for deformable shapes using localized spectral convolutional networks. In _Computer Graphics Forum_ , Vol. 34. Wiley Online Library, 13–23.
  * Bowman and Shenton (2004) KO Bowman and LR Shenton. 2004.  Estimation: Method of moments.  _Encyclopedia of statistical sciences_ 3 (2004).
  * Chen et al. (2020) Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. 2020.  Simple and deep graph convolutional networks. In _International Conference on Machine Learning_. PMLR, 1725–1735.
  * Chien et al. (2020) Eli Chien, Jianhao Peng, Pan Li, and Olgica Milenkovic. 2020\.  Adaptive Universal Generalized PageRank Graph Neural Network. In _International Conference on Learning Representations_.
  * Corso et al. (2020) Gabriele Corso, Luca Cavalleri, Dominique Beaini, Pietro Liò, and Petar Veličković. 2020\.  Principal neighbourhood aggregation for graph nets.  _Advances in Neural Information Processing Systems_ 33 (2020), 13260–13271.
  * Dai et al. (2022) Enyan Dai, Wei Jin, Hui Liu, and Suhang Wang. 2022\.  Towards robust graph neural networks for noisy graphs with sparse labels. In _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_. 181–191.
  * Defferrard et al. (2016) Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. 2016.  Convolutional neural networks on graphs with fast localized spectral filtering.  _Advances in neural information processing systems_ 29 (2016).
  * Dehmamy et al. (2019) Nima Dehmamy, Albert-László Barabási, and Rose Yu. 2019.  Understanding the representation power of graph neural networks in learning graph topology.  _Advances in Neural Information Processing Systems_ 32 (2019).
  * Derr et al. (2020) Tyler Derr, Yao Ma, Wenqi Fan, Xiaorui Liu, Charu Aggarwal, and Jiliang Tang. 2020\.  Epidemic graph convolutional network. In _Proceedings of the 13th International Conference on Web Search and Data Mining_. 160–168.
  * Du et al. (2021a) Lun Du, Fei Gao, Xu Chen, Ran Jia, Junshan Wang, Jiang Zhang, Shi Han, and Dongmei Zhang. 2021a.  TabularNet: A Neural Network Architecture for Understanding Semantic Structures of Tabular Data. In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_. 322–331.
  * Du et al. (2021b) Lun Du, Xiaozhou Shi, Qiang Fu, Hengyu Liu, Shi Han, and Dongmei Zhang. 2021b.  GBK-GNN: Gated Bi-Kernel Graph Neural Networks for Modeling Both Homophily and Heterophily.  _arXiv preprint arXiv:2110.15777_ (2021).
  * Dwivedi et al. (2021) Vijay Prakash Dwivedi, Anh Tuan Luu, Thomas Laurent, Yoshua Bengio, and Xavier Bresson. 2021\.  Graph neural networks with learnable structural and positional representations.  _arXiv preprint arXiv:2110.07875_ (2021).
  * Fettal et al. (2022) Chakib Fettal, Lazhar Labiod, and Mohamed Nadif. 2022\.  Efficient Graph Convolution for Joint Node Representation Learning and Clustering. In _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_. 289–297.
  * Fisher (1936) Ronald A Fisher. 1936\.  The use of multiple measurements in taxonomic problems.  _Annals of eugenics_ 7, 2 (1936), 179–188.
  * Guo et al. (2022) Zhihao Guo, Feng Wang, Kaixuan Yao, Jiye Liang, and Zhiqiang Wang. 2022.  Multi-Scale Variational Graph AutoEncoder for Link Prediction. In _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_. 334–342.
  * Hamilton et al. (2017) William L Hamilton, Rex Ying, and Jure Leskovec. 2017\.  Inductive representation learning on large graphs. In _Proceedings of the 31st International Conference on Neural Information Processing Systems_. 1025–1035.
  * Jin et al. (2021) Wei Jin, Tyler Derr, Yiqi Wang, Yao Ma, Zitao Liu, and Jiliang Tang. 2021.  Node similarity preserving graph convolutional networks. In _Proceedings of the 14th ACM international conference on web search and data mining_. 148–156.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Kipf and Welling (2016) Thomas N Kipf and Max Welling. 2016.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ (2016).
  * Klicpera et al. (2018) Johannes Klicpera, Aleksandar Bojchevski, and Stephan Günnemann. 2018.  Predict then propagate: Graph neural networks meet personalized pagerank.  _arXiv preprint arXiv:1810.05997_ (2018).
  * Lin et al. (2020) Wenqing Lin, Feng He, Faqiang Zhang, Xu Cheng, and Hongyun Cai. 2020.  Initialization for network embedding: A graph partition approach. In _Proceedings of the 13th International Conference on Web Search and Data Mining_. 367–374.
  * Lin et al. (2014) Wenqing Lin, Xiaokui Xiao, and Gabriel Ghinita. 2014\.  Large-scale frequent subgraph mining in mapreduce. In _2014 IEEE 30th International Conference on Data Engineering_. IEEE, 844–855.
  * Liu et al. (2020) Meng Liu, Hongyang Gao, and Shuiwang Ji. 2020.  Towards Deeper Graph Neural Networks. In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. ACM.
  * Luo et al. (2022) Zihan Luo, Jianxun Lian, Hong Huang, Hai Jin, and Xing Xie. 2022.  Ada-GNN: Adapting to Local Patterns for Improving Graph Neural Networks. In _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_. 638–647.
  * Ma et al. (2022) Xiaojun Ma, Qin Chen, Yuanyi Ren, Guojie Song, and Liang Wang. 2022.  Meta-Weight Graph Neural Network: Push the Limits Beyond Global Homophily. In _Proceedings of the ACM Web Conference 2022_. 1270–1280.
  * Monti et al. (2018) Federico Monti, Karl Otness, and Michael M Bronstein. 2018\.  Motifnet: a motif-based graph convolutional network for directed graphs. In _2018 IEEE Data Science Workshop (DSW)_. IEEE, 225–228.
  * Natekar and Sharma (2020) Parth Natekar and Manik Sharma. 2020.  Representation based complexity measures for predicting generalization in deep learning.  _arXiv preprint arXiv:2012.02775_ (2020).
  * Neyshabur et al. (2017) Behnam Neyshabur, Srinadh Bhojanapalli, David McAllester, and Nati Srebro. 2017.  Exploring generalization in deep learning.  _Advances in neural information processing systems_ 30 (2017).
  * Papoulis and Pillai (2002) Athanasios Papoulis and S Unnikrishna Pillai. 2002.  _Probability, random variables, and stochastic processes_.  Tata McGraw-Hill Education.
  * Pei et al. (2020) Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. 2020.  Geom-gcn: Geometric graph convolutional networks.  _arXiv preprint arXiv:2002.05287_ (2020).
  * Traud et al. (2012) Amanda L Traud, Peter J Mucha, and Mason A Porter. 2012\.  Social structure of Facebook networks.  _Physica A: Statistical Mechanics and its Applications_ 391, 16 (2012), 4165–4180.
  * Veličković et al. (2017) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. 2017.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ (2017).
  * Wang et al. (2022) Tao Wang, Di Jin, Rui Wang, Dongxiao He, and Yuxiao Huang. 2022.  Powerful graph convolutional networks with adaptive propagation mechanism for homophily and heterophily. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 36. 4210–4218.
  * Wang et al. (2019) Yun Wang, Lun Du, Guojie Song, Xiaojun Ma, Lichen Jin, Wei Lin, and Fei Sun. 2019.  Tag2Gauss: Learning Tag Representations via Gaussian Distribution in Tagged Networks.. In _IJCAI_. 3799–3805.
  * Xu et al. (2019) Bingbing Xu, Huawei Shen, Qi Cao, Yunqi Qiu, and Xueqi Cheng. 2019.  Graph wavelet neural network.  _arXiv preprint arXiv:1904.07785_ (2019).
  * Xu et al. (2018a) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. 2018a.  How powerful are graph neural networks?  _arXiv preprint arXiv:1810.00826_ (2018).
  * Xu et al. (2018b) Keyulu Xu, Chengtao Li, Yonglong Tian, Tomohiro Sonobe, Ken-ichi Kawarabayashi, and Stefanie Jegelka. 2018b.  Representation learning on graphs with jumping knowledge networks. In _International Conference on Machine Learning_. PMLR, 5453–5462.
  * Yang et al. (2016) Zhilin Yang, William Cohen, and Ruslan Salakhudinov. 2016\.  Revisiting semi-supervised learning with graph embeddings. In _International conference on machine learning_. PMLR, 40–48.
  * Yao et al. (2022) Di Yao, Haonan Hu, Lun Du, Gao Cong, Shi Han, and Jingping Bi. 2022.  TrajGAT: A Graph-based Long-term Dependency Modeling Approach for Trajectory Similarity Computation. In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_. 2275–2285.
  * Ying et al. (2021) Chengxuan Ying, Tianle Cai, Shengjie Luo, Shuxin Zheng, Guolin Ke, Di He, Yanming Shen, and Tie-Yan Liu. 2021\.  Do Transformers Really Perform Badly for Graph Representation?  _Advances in Neural Information Processing Systems_ 34 (2021).
  * Zeng et al. (2019) Hanqing Zeng, Hongkuan Zhou, Ajitesh Srivastava, Rajgopal Kannan, and Viktor Prasanna. 2019\.  Graphsaint: Graph sampling based inductive learning method.  _arXiv preprint arXiv:1907.04931_ (2019).
  * Zhang and Chen (2018) Muhan Zhang and Yixin Chen. 2018.  Link prediction based on graph neural networks.  _Advances in Neural Information Processing Systems_ 31 (2018), 5165–5175.
  * Zhao et al. (2021) Tianxiang Zhao, Xiang Zhang, and Suhang Wang. 2021\.  Graphsmote: Imbalanced node classification on graphs with graph neural networks. In _Proceedings of the 14th ACM international conference on web search and data mining_. 833–841.
  * Zhou et al. (2018) Lekui Zhou, Yang Yang, Xiang Ren, Fei Wu, and Yueting Zhuang. 2018.  Dynamic network embedding by modeling triadic closure process. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 32.
  * Zhu et al. (2020) Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra. 2020\.  Beyond homophily in graph neural networks: Current limitations and effective designs.  _Advances in Neural Information Processing Systems_ 33 (2020), 7793–7804.