# Quantum Extreme Learning of molecular potential energy surfaces and force
fields

Gabriele Lo Monaco  Università degli Studi di Palermo, Dipartimento di Fisica
e Chimica - Emilio Segrè, via Archirafi 36, I-90123 Palermo, Italy  Marco
Bertini  Università degli Studi di Palermo, Dipartimento di Fisica e Chimica -
Emilio Segrè, via Archirafi 36, I-90123 Palermo, Italy  Salvatore Lorenzo
Università degli Studi di Palermo, Dipartimento di Fisica e Chimica - Emilio
Segrè, via Archirafi 36, I-90123 Palermo, Italy  G. Massimo Palma  Università
degli Studi di Palermo, Dipartimento di Fisica e Chimica - Emilio Segrè, via
Archirafi 36, I-90123 Palermo, Italy  NEST, Istituto Nanoscienze-CNR, Piazza
S. Silvestro 12, 56127 Pisa, Italy

###### Abstract

Quantum machine learning algorithms are expected to play a pivotal role in
quantum chemistry simulations in the immediate future. One such key
application is the training of a quantum neural network to learn the potential
energy surface and force field of molecular systems. We address this task by
using the quantum extreme learning machine paradigm. This particular
supervised learning routine allows for resource-efficient training, consisting
of a simple linear regression performed on a classical computer. We have
tested a setup that can be used to study molecules of any dimension and is
optimized for immediate use on NISQ devices with a limited number of native
gates. We have applied this setup to three case studies: lithium hydride,
water, and formamide, carrying out both noiseless simulations and actual
implementation on IBM quantum hardware. Compared to other supervised learning
routines, the proposed setup requires minimal quantum resources, making it
feasible for direct implementation on quantum platforms, while still achieving
a high level of predictive accuracy compared to simulations. Our encouraging
results pave the way towards the future application to more complex molecules,
being the proposed setup scalable.

##  1 Introduction

The simulation of complex phenomena involving large biochemical systems, such
as protein folding or enzyme activation, requires accurate predictions of
energies and inter-atomic forces for different molecular conformations.
Standard techniques such as molecular dynamics (MD) require _on-the-fly_
computing the energy and inter-atomic forces of the target biochemical system
at each step of the simulation. In most MD-like schemes, the Potential Energy
Surface (PES) and Force Fields (FF) are determined on empirical grounds or on
databases. In contrast, a higher degree of accuracy demands an ab initio
calculation of the PES, typically using density functional theory (DFT) [1].
The use of DFT on the fly leads to a much longer simulation running time and,
due to its huge computational cost, it turns out to be a viable route only for
short-time simulations (picoseconds) of small molecules.

A precise knowledge of the functional relation between the molecular geometry
and the corresponding energy on the PES, predicted by the DFT, would
incredibly speed up ab-initio MD simulations. As Behler and Parrinello
observed in [2], this functional relation may be well approximated by neural
networks (NNs) after appropriate training on the pairs of molecular
conformations and their DFT energies (or forces). The generation of the PES
and FFs using various machine learning (ML) paradigms has been further
investigated in [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. Given the
intrinsic quantum nature of the PES, it is natural to believe that quantum
algorithms may indeed help [17, 18]. Through the implementation of a
supervised learning setup where the role of the classical NN is played by a
parameterized quantum circuit (PQC), the authors of [19, 20] provided the
first evidence that quantum machine learning (QML) routines may offer improved
and more accurate PES predictions.

Since PQCs can be implemented on NISQ devices, they have been extensively
utilized in recent times for different purposes. However, they are not the
unique viable QML paradigm for supervised learning. An attractive alternative
is provided by quantum extreme learning machines (QELMs). This setup may be
advantageous for various reasons. First, it requires reduced training efforts,
as in the QELM paradigm no parameters in the quantum platform need to be tuned
and the training consists only of a simple classical linear regression.
Second, quantum extreme learning machines can be implemented on different
physical platforms [21, 22, 23, 24, 25, 26, 27, 28] and not necessarily on
gate-based quantum devices. We propose to use QELM to speed up ab initio
molecular dynamics computations. We show that a QELM trained on small sets of
geometry-energy pairs is able to efficiently predict the energy of a set of
test geometries with great accuracy.
The QELM setup that we propose significantly reduces the quantum resources
required for the training process compared to other supervised learning
paradigms such as VQE (as discussed in [19]). By quantum resources, we refer
to either the depth of the quantum circuits employed, as reported in Table 1
and Table 2 and the number of runs on the quantum platform. In fact, in QELM
the number of runs only depends on the size of the training set and the chosen
shot statistics; the optimization stage only consists of simple linear
regression on the output of the quantum circuit, and no parameter on the
quantum device needs to be optimized. In VQE, instead, the number of runs also
depends on the number of iterations necessary for the optimization to converge
and on the number of parameters in the quantum circuit to be optimized.

Section 2 contains a self-consistent review of the QELM paradigm. The details
about our particular realization of QELM adapted to PES and FF generations are
all collected in section 3.1 and section 3.2. Finally, in section 3.3, we
collect all the results about the simulations and quantum-hardware realization
of the QELM training in the case of lithium hydride, water, and formamide.

##  2 Basics of QELM

###  2.1 QELM as supervised training of a POVM outcome

The QELM paradigm with classical input data is a supervised learning model.
Its task is to _learn_ a function
$f:\mathbb{R}^{\mathcal{X}}\rightarrow\mathbb{R}^{\mathcal{Y}}$ that well
approximates the actual functional relation between inputs and outputs of some
training set $\\{\boldsymbol{x}^{\rm tr}_{i},\boldsymbol{y}^{\rm
tr}_{i}\\}_{i=1}^{M_{\rm tr}}$. The quality of the training is assessed on a
separate test set of previously unseen pairs $\\{\boldsymbol{x}^{\rm
test}_{i},\boldsymbol{y}^{\rm test}_{i}\\}_{i=1}^{M_{\rm test}}$.

In the case at hand, the training set consists of the coordinates
parameterizing the conformations of a given molecular system and the
corresponding energy and forces, computed classically with standard quantum
chemistry routines via Hartree-Fock (HF) approximation or DFT. While a
possible accuracy threshold may be the chemical accuracy $1.6\cdot
10^{-3}\,{\rm Hartree}$, it should be emphasized that the chemical accuracy
plays a role only in the comparison with experimental results, where energy
differences related to chemical processes are measured, and not in the
evaluation of the absolute energy of a molecular geometry obtained through
calculations. Even if chemical accuracy is a good target precision to keep in
mind, the goodness of a learning training should be actually assessed on the
typical fluctuations of the variable $\boldsymbol{y}$ in the test dataset. If
the error is well below such typical fluctuations, we have a good chance to
accurately predict the real difference between the energy or the force of two
given sampled geometries.

In a QELM, the classical $\boldsymbol{x}$ data are encoded in the quantum
state of some _input qubits_ , initially prepared in a state $\rho_{0}$. The
encoding is performed by applying a parameterized unitary transformation:

$\rho_{\boldsymbol{x}}\,=\,U(\boldsymbol{x})\rho_{0}\,U^{\dagger}(\boldsymbol{x})$
(1)

After such encoding, a POVM measurement
$\boldsymbol{\mu}=\\{\mu_{a}\\}_{a=1}^{\Sigma}$ is performed on the input
qubits.

Each outcome $a$ is measured with probability

$p_{a}(\boldsymbol{x})=\Tr(\mu_{a}\,\rho_{\boldsymbol{x}})\,.$ (2)

In the QELM protocol, the function $f$ is approximated as:

$f(\boldsymbol{x})\,=\,W\boldsymbol{p}(\boldsymbol{x})\,,$ (3)

where $W$ is a $\mathcal{Y}\times\Sigma$ matrix. The optimization of the
unknown matrix $W$ in order to fit the training data can be performed using
standard techniques and it consists of a simple linear regression.

$W=\boldsymbol{Y}_{\rm tr}\,\boldsymbol{P}_{\rm tr}^{+}$ (4)

where $\boldsymbol{P}_{\rm tr}$ is a $\Sigma\otimes M_{\rm tr}$ matrix whose
columns are the probability vectors $\boldsymbol{p}(\boldsymbol{x}_{i}^{\rm
tr})$ of the training set. The superscript $+$ denotes the Moore-Penrose
pseudoinverse and $\boldsymbol{Y}_{\tr}$ is a $\mathcal{Y}\otimes M_{\rm tr}$
matrix whose columns are the output $\boldsymbol{y}_{\rm tr}$ of the training
set. The matrix $W$ is thus the one minimizing the root mean square error
(RMSE) on the training set $\epsilon_{\rm tr}=||\boldsymbol{Y}_{\rm
tr}-W\boldsymbol{P}_{\rm tr}||_{F}/\sqrt{M_{\rm tr}}$ where $||\cdot||_{F}$ is
the Frobenius norm. What is relevant is the root mean square error on the test
set:

$\epsilon_{\rm QELM}=\frac{1}{\sqrt{M_{\rm test}}}||\boldsymbol{Y}_{\rm
test}-W\boldsymbol{P}_{\rm test}||_{F}\,.$ (5)

This is an estimate of the ability of the trained QELM to provide accurate
predictions on a new set of test input data. As already mentioned, the root
mean square error alone is not a good quantifier and it should be compared
with the typical fluctuations in the output values $\boldsymbol{y}$ if we
sample two random elements from the dataset. This scale is measured by the
variance

${\rm var}(\boldsymbol{Y})=\frac{1}{M_{\rm test}}||\boldsymbol{Y}_{\rm
test}-\overline{\boldsymbol{Y}}_{\rm test}||^{2}_{F}\,,$ (6)

where $\overline{\boldsymbol{Y}}_{\rm test}$ is the $\cal{Y}\times M_{\rm
test}$ matrix whose column are all equal to the mean value
$\overline{\boldsymbol{y}}_{\rm test}$. If $\epsilon_{\rm QELM}$ is much
smaller of this rooted variance, this means that we will able to predict the
variation of the variable $\boldsymbol{y}$ with great precision.

Let us point out that it is always possible to rewrite (3) in the following
way:

$f(\boldsymbol{x})=W\boldsymbol{p}(\boldsymbol{x})=\Tr\left(\tilde{\boldsymbol{\mu}}\rho_{x}\right)$
(7)

where $\tilde{\boldsymbol{\mu}}=W\boldsymbol{\mu}$ is a vector of
$\mathcal{Y}$ operators acting on the Hilbert space of the input qubits. This
means that we are actually approximating the function $f$ as a vector of
expectation values. On real quantum platforms, the probabilities (2) are
estimated performing a finite number of shots (_finite statistics_). This
implies that a statistical error is always present and how to minimize it is
an open problem of crucial relevance [29, 30, 31, 32]. Furthermore, the choice
of the embedding unitaries $U(\boldsymbol{x})$ and the choice of POVM
$\boldsymbol{\mu}$ strongly affects the performances of the QELM both at
infinite and finite statistics.

###  2.2 Implementing informationally complete POVMs with a reservoir

As we have pointed out, QELM is equivalent to estimating the expectation value
of some observable. Looking at (7), this is possible only if such observable
can be decomposed as a linear combination of the effect of the POVM and this
is possible for any operator only if the the POVM is informationally complete
(IC). From an experimental point of view, implementing an IC-POVM is in
principle a hard task. To overcome this issue, one may let the input qubits
interact with some ancillae (the _reservoir_) prepared in a fixed state
$\rho_{R}$. Then, an easy-to-implement measurement $\boldsymbol{\nu}$ is
performed. For instance, some qubits can be measured in the computational
basis. The string of probabilities obtained with this procedure can be written
as:

$p_{a}(\boldsymbol{x})=\Tr\left(\nu_{a}\,\mathcal{U}\rho_{\boldsymbol{x}}\otimes\rho_{R}\,\mathcal{U}^{\dagger}\right)\,=\,\Tr\left(\Phi^{\dagger}(\nu_{a})\rho_{x}\right)\,,$
(8)

where we defined the quantum channel
$\Phi(\rho)=\mathcal{U}\rho\otimes\rho_{R}\mathcal{U}^{\dagger}$ and
$\Phi^{\dagger}$ denotes its dual. Comparing this expression with (2), it is
evident that the role of the interaction of the input qubits with the
reservoir is to generate an effective IC-POVM
$\tilde{\boldsymbol{\mu}}=\Phi^{\dagger}(\boldsymbol{\nu})$ on the input
qubits. The reconstruction performances of QELM are affected by the choice of
$\boldsymbol{\nu}$ and by the scrambling properties of the interaction
$\mathcal{U}$ between the input qubits and the ancillae [33]. It is common in
the literature to assume the measurement $\boldsymbol{\nu}$ to act only on the
reservoir. However, it is evident from our discussion that this is not
necessary. In the following, we will assume the full set of input qubits and
ancillae to be measured in the computation basis. In this way, the number of
ancillae required to generate an effective IC-POVM is minimized.

###  2.3 Expressivity and QELM

Different choices of embedding of the classical input may influence the
capacity of the QELM to reconstruct the functional relation between the
elements of the training set. In order to understand how this choice may
affect the performances, it is useful to observe that the QELM paradigm shows
some similarities with the quantum kernel method, a feature that seems to be
common to all quantum supervised learning routines [33, 34]. In other words,
we embed the classical data $\boldsymbol{x}$ into an element of a higher
dimensional space, $\rho_{\boldsymbol{x}}\in\mathcal{D}$. The input Hilbert
space is endowed with an inner product defining the kernel of the model:

$\kappa(\boldsymbol{x},\boldsymbol{y})=\Tr\left(\rho_{\boldsymbol{x}}\rho_{\boldsymbol{y}}\right)\,.$
(9)

Different kernels $\kappa(\boldsymbol{x},\boldsymbol{y})$ allow to approximate
different kinds of functions. For instance, a popular choice is the _rotation
encoding_ :

$U(\boldsymbol{x})=\prod_{i=1}^{d_{\rm
in}}\exp\left(-\frac{i}{2}x_{i}\sigma^{(i)}_{x}\right)\,,$ (10)

where $\sigma^{(i)}_{x}$ is the Pauli matrix $x$ acting on the $i$-th input
qubit. The corresponding kernel is the _translation invariant squared cosine
kernel_ [34]:

$\kappa(\boldsymbol{x},\boldsymbol{y})=\prod_{i=1}^{d_{\rm
in}}|\cos(x_{i}-y_{i})|^{2}$ (11)

The rotation embedding can be considered as a particular instance of the more
general _Fourier embedding_ :

$U(\boldsymbol{x})=e^{-x_{\mathcal{X}}G}W^{(\mathcal{X})}\dots
W^{(2)}e^{-x_{1}G}W^{(1)}$ (12)

where $W^{(1)},W^{(2)},\dots$ are random unitary matrices acting on the input
qubits, and $G$ is a Hermitian operator, for instance, a Pauli string. The
name of this embedding comes from the peculiar form of the associated kernel:

$\kappa(\boldsymbol{x},\boldsymbol{y})=\sum_{\boldsymbol{n},\boldsymbol{m}}\alpha_{\boldsymbol{n},\boldsymbol{m}}e^{i(\boldsymbol{n}\cdot\boldsymbol{x}-\boldsymbol{m}\cdot\boldsymbol{y})}$
(13)

where $c_{\boldsymbol{m},\boldsymbol{n}}$ are complex coefficients and
$\boldsymbol{n},\boldsymbol{m}\in\mathbb{\mathbb{R}}^{\mathcal{X}}$. The
Fourier frequencies $\boldsymbol{m}$ are proportional to the gaps
$\lambda_{i}-\lambda_{j}$, where $\lambda_{i,j}$ are eigenvalues of $G$ [35].
The number of distinct Fourier frequencies is bounded from above by
$2^{2N-1}-1$, where $N$ is the number of encoding qubits. As the number of
frequencies increases, the model becomes more and more expressive. For this
reason, studying the kernel associated with the embedding is a useful way to
assess expressiveness [33, 34, 36, 37]. As shown in [35], the Fourier encoding
models are universal approximators for appropriate choices of generator $G$.

Figure 1: The QELM setup used for PES prediction with $4$ input qubits. The
mixing unitary $\tilde{W}$ is made of three blocks, each of the same form as
$W^{(1)}$, with random rotation angles drawn from the same ensemble. The two-
qubit entangling gates in each layer are ECR gates.

Figure 2: Parameterization of the geometry of a single molecule of lithium
hydride.

Figure 3: Parameterization of the geometry of a single molecule of water.

Figure 4: Parameterization of the geometry of a single molecule of formamide.
The two ${\rm NO}$ bond lengths are taken equal and parameterized by the same
coordinate $r_{4}.$

##  3 PES predictions with QELM: simulations

In this section we train a QELM to learn the PES and the FF for three
different molecules of increasing complexity: lithium hydride (${\rm LiH}$),
water (${\rm H}_{2}{\rm O}$) and formamide (${\rm HCONH}_{2}$). The datasets
used during the training are generated in the same way for all the
investigated molecules and their generation is discussed in greater detail in
section 3.1. The QELM model used can be adapted to molecules with an arbitrary
number of atoms and is presented in section 3.2.

###  3.1 Datasets

The conformations of any molecule are represented in terms of generalized
coordinates, either bond lengths or bond angles, that are usually used in the
Z-matrix description of any species. The lithium hydride conformations are
described in terms of a single coordinate, the bond length (Figure 4). The
conformations of a water molecule are parameterized by three coordinates: the
${\rm HO}$ bond lengths $r_{1},r_{2}$ and the oxygen and the ${\rm HOH}$ bond
angle $\phi$ (Figure 4). The parameterization of formamide deserves a few more
words. Such organic molecule is a particular instance of amides, the smallest
molecules presenting a peptide bond. They have particular biological
relevance: proteins are polymers whose units are connected by amide groups.
The distinctive strength of the peptide (double) bonds forces the formamide to
be a planar molecule in its ground state configuration. For such reason, we
parameterized a formamide molecule assuming the geometry to be planar and the
two ${\rm NH}$ bond lengths are assumed to be the same. The remaining
conformations are parameterized in terms of $8$ generalized coordinates
(instead of $12$) and their choice is shown in Figure 4.

The training and test configurations are then obtained by sampling random
tuples of generalized coordinates in given intervals. In the three cases under
consideration, these intervals are centered around the ground state
configuration in order to focus on the geometries that are typically explored
during a molecular dynamics simulations. In fact, for our three molecular
species, the PES are characterized by a deep potential well which it’s hard to
escape. Given that the PES of ${\rm LiH}$ is actually a one-dimensional curve,
in this case, we explore a wider part of the configuration landscape sampling
the bond lengths from a larger interval than ${\rm H}_{2}{\rm O}$ or ${\rm
HCONH}_{2}$. The electronic structure calculations performed to generate the
used dataset were conducted using the PySCF software implemented in Pennylane
package (${\rm LiH}$) or Gaussian 16 software [38] (${\rm H}_{2}{\rm O}$ and
${\rm HCONH}_{2}$). Since the main focus of this work is to assess the
agreement of QELM with generic PES data obtained from traditional quantum
chemistry methods, we are not restricted to a specific theoretical framework.
For this reason, the level of theory chosen to investigate the different
chemical species changes as a consequence of both the complexity of the
molecule and the necessity of a properly sized dataset.

In the case of lithium hydride and water, calculations were performed using
DFT with M06 exchange-correlation functional [39] and STO-3G (cc-pVDZ) basis
set. For ${\rm LiH}$, the generated dataset is generated performing $170$
single point calculation on as many geometries, obtained by choosing the bond
length values randomly in the range 0.9Å $\leq r\leq$ 4.5Å. The QELM is
trained over $50$ geometries while $120$ configurations are left for test. For
${\rm H}_{2}{\rm O}$, the dataset is generated performing $900$ single point
calculations on as many geometries. Each geometry is obtained by choosing each
value for the three internal coordinates at random, in a range centered on
their respective equilibrium values, i.e. 0.964 $\pm$ 0.2Å for distances, and
102.792∘ $\pm$ 13.0∘ for the ${\rm HOH}$ angle. The QELM is trained over $300$
elements of the whole dataset, while $600$ configurations are left for test.
In the formamide case, the dataset consists of approximately 6500 structures
(even if much less are actually used for the training). Each of the geometries
that constitute the dataset was generated by randomly sampling values for the
aforementioned bonds and angles within specific intervals. In particular,
these intervals are centered on the value of the respective internal
equilibrium coordinates, with ranges of: i) $\pm 0.10$ Å for bonds involving a
hydrogen atom; ii) $\pm 0.15$ Å for the remaining bonds; and iii) $\pm
8.00^{\circ}$ for all angles. For each structure, a calculation was performed
at the HF/SVP level of theory to obtain the electronic energy and atomic force
values.

${\rm LiH}$ |  # shots |  RMSE(E) |  RMSE($F_{1}$) |  # qubits |  depth
---|---|---|---|---|---
statevector | $\infty$ | $6.9\cdot 10^{-5}$ | $4.2\cdot 10^{-4}$ | $5$ | 27
QASM simulator | $4\cdot 10^{4}$ | $9.3\cdot 10^{-3}$ | $1.8\cdot 10^{-2}$ | $4$ | 26
IBM_BRISBANE | $4\cdot 10^{4}$ | $3.2\cdot 10^{-2}$ | $5.0\cdot 10^{-2}$ | $4$ | 26
${\rm H}_{2}{\rm O}$ |  # shots |  RMSE(E) |  RMSE($F_{1}$) |  # qubits |  depth
statevector | $\infty$ | $2.0\cdot 10^{-6}$ | $9.3\cdot 10^{-6}$ | $7$ | 41
QASM simulator | $2\cdot 10^{4}$ | $3.6\cdot 10^{-3}$ | $1.2\cdot 10^{-2}$ | $5$ | 39
IBM_BRISBANE | $2\cdot 10^{4}$ | $5.2\cdot 10^{-3}$ | $2.2\cdot 10^{-2}$ | $5$ | 39
${\rm HCONH}_{2}$ |  # shots |  RMSE(E) |  RMSE($F_{1}$) |  # qubits |  depth
statevector | $\infty$ | $3.4\cdot 10^{-4}$ | $9.7\cdot 10^{-4}$ | $9$ | 73
quasm simulator | $2\cdot 10^{4}$ | $7.0\cdot 10^{-3}$ | $2.1\cdot 10^{-2}$ | $7$ | 71
IBM_BRISBANE | $2\cdot 10^{4}$ | $1.3\cdot 10^{-2}$ | $4.4\cdot 10^{-2}$ | $7$ | 71
Table 1:  Summary of the performances of the QELM trained to learn PESs and
FFs of a molecule of ${\rm LiH}$, ${\rm H}_{2}{\rm O}$ and ${\rm HCONH}_{2}$.
The RMSE is measured in ${\rm Ha}$ for the energy and ${\rm Ha}/$Å for the
force.

###  3.2 QELM setup

The gate parameters in quantum hardware are usually rotation angles. For this
reason, we need to perform a rescaling of the coordinate, with particular
attention to radial coordinates. Each bond length is rescaled with respect to
a reference one $\overline{r}$. The reference distances are chosen to be
$\overline{r}_{{\rm LiH}}=6$Å and $\overline{r}_{{\rm H}_{2}{\rm
O}}=\overline{r}_{{\rm HCONH}_{2}}=2$Å. The bond angles are instead rescaled
by a factor $2$. The input data of water are thus of the following form:

$\boldsymbol{x}=\left\\{\frac{r_{1}}{\overline{r}}\pi\,,\,\frac{r_{2}}{\overline{r}}\pi\,,\,\frac{\phi}{2}\right\\}\,.$
(14)

with analogous expressions for lithium hydride and formamide. The string
$\boldsymbol{x}$ is fed in the state of $N$ qubits using a Fourier encoding
(12) with:

$G=\frac{1}{2}\sum_{i=1}^{N}\sigma_{z}^{(i)}$ (15)

This amounts to act on each input qubit with an $R_{z}$ gate. Observe that the
number of input qubits of the encoding is not fixed and can be increased at
wish in order to improve performances. The random unitaries
$W^{(1)},W^{(2)},\dots$ in the Fourier encoding are made of one layer of the
form depicted in Figure 1, where $\theta_{i}$ are random rotation angles
uniformly drawn in the interval $(0,\pi/2)$. All the input qubits are measured
in the computational basis, with no addition of ancillae. Let us also stress
that the Fourier encoding as presented in (12) can be adapted to any molecule
or system of molecules with an arbitrary number of generalized coordinates.

The kernel associated with the choice of generator (15) is a Fourier kernel
(13) with integer frequencies $-N,-N+1,\dots,N-1,N$. Even if the number of
frequencies scales only linearly with the number of qubits, it must be
stressed that the choice of embedding and gates in the setup is driven by the
wish to minimize the depth of the quantum circuit once it is implemented on
real quantum hardware. In fact, as it will be clear in the following, the
accuracy of the PES predictions grows as the number of encoding qubits
increases, reflecting the wider spectrum of Fourier frequencies of the model
(see Figure 5).

In the following, we will present the actual experimental implementation of
the proposed setup on the IBM_BRISBANE device, a quantum processor with $127$
superconducting transmon qubits. Such a device admits thee native single-qubit
gates ($X$, $\sqrt{X}$ and ${\rm R_{z}}(\phi)$ for any angle) and a single
two-qubit gate, the _Echoed cross-resonance_ gate (ECR):

${\rm ECR}_{12}={\rm R_{xx}}(\pi/2){\rm X}_{2}$ (16)

The ECR gate is equivalent to the CNOT gate up to single-qubit rotations. Any
other gates must be decomposed as a combination of native gates and this may
cause a huge increase in the circuit depth. For instance:

$\displaystyle{\rm R_{x}}(\phi)$ $\displaystyle={\rm
R_{z}}\left(\frac{\pi}{2}\right)\sqrt{X}\,{\rm
R_{z}}(\pi+\phi)\,\sqrt{X}\,{\rm R_{z}}\left(\frac{5\pi}{2}\right)$ (17)
$\displaystyle{\rm R_{y}}(\phi)$ $\displaystyle={\rm
R_{z}}(-\pi)\,\sqrt{X}\,{\rm R_{z}}(\pi-\phi)\,\sqrt{X}$ (18)

Since the ${\rm R_{y}}$ gates are slightly cheaper in terms of native gates,
they should be preferred in order to minimize the circuit depth. For this
reason, our circuit is built out of ${\rm R_{z,y}}$ and ECR gates only.

###  3.3 Results

For each atomic species that we considered, we performed three different kinds
of analysis: first, we simulated the performances of the QELM in the ideal
case of infinite statistics and absence of noise. We refer to this case as
_statevector_ simulation. Then we reduce the level of approximation simulating
the QELM training with finite statistics but in noiseless conditions. To
address this task, we use the _QASM simulator_ of the Qiskit package [40].
Finally, we implement the QELM setup on IBM_BRISBANE, taking thus into account
also the actual noise present on NISQ devices. It must be stressed that we do
not perform any error mitigation on the quantum-hardware output data.

The statevector simulation is the ideal playground to test the achievable
precision of the QELM prediction when the parameters of the setup are varied.
The only freedom left in our configuration is in the number of encoding qubits
and the dimension of the training set. In Figure 5 we plot the behavior of the
statevector RMSE for the energy for all the molecules under investigation when
both $N$ and $M_{\rm train}$ are varied. As one can observe, as the number of
qubits increases, the performance of the QELM improves. This is consistent
with the idea that a larger reservoir can be used to approximate a wider range
of functions [33]. What is interesting is that the RMSE always reaches a
plateau above a certain dimension of the training set. Above this threshold
value, any further information does not significantly affect the accuracy of
the predictions on the test set. The reason relies on the linear regression
that is carried out in post-processing in the QELM paradigm. As explained in
Figure 2, the training consists of a single linear regression once the matrix

$W=\boldsymbol{Y}_{\rm tr}\boldsymbol{\,}P_{\rm tr}^{+}$ (19)

is determined. The reservoir dependence of the training is all contained in
the probability matrix $\boldsymbol{P}_{\rm tr}^{+}$. It is well known that
some information about the accuracy of the training can be extracted by
looking at the singular values of the probability matrix only [29, 41].
However, the number of singular values of $\boldsymbol{P}_{\rm tr}^{+}$ is
always less than $\min(2^{N},M_{\rm tr})$. This implies that, as long as the
number of qubits is kept fixed, the larger the value of $M_{\rm tr}$ is
compared to $2^{N}$, the less any new training element will impact
performance. Looking at Figure 5, the RMSE plateau is reached for:

$2^{N}<M_{\rm tr}\lesssim 2^{N+1}\,.$ (20)

In the first row of each sub-table of Table 1, we collect the results of
training for the three atomic species for a specific choice of $N$. The RMSE
of the energy is in all three cases well below the typical energy fluctuations
in the dataset, estimated by the square root of the variance, $\sqrt{{\rm
var}(\boldsymbol{E}_{\rm tr})}$, where $\boldsymbol{E}_{\rm tr}$ is the vector
whose elements are the training energies. In all the three cases, ${\rm LiH}$,
${\rm H}_{2}{\rm O}$ and ${\rm HCONH}_{2}$, the variance is of order
$2\cdot{10}^{-2}\,{\rm Ha}$. For the force associated to bond lengths instead,
$\sqrt{{\rm var}(\boldsymbol{F}_{r})}\approx 10^{-1}{\rm Ha}/$Å.

It is interesting to compare these results with the performances (in the
infinite-statistics, noiseless ideal case) of the VQE routine proposed in
[19], where an analogous study of ${\rm LiH}$ and ${\rm H}_{2}{\rm O}$ is
carried out Table 2.

Table 2: Performances (statevector simulation) of the VQE routine proposed in [19] to learn PES and FF of LiH and H2O, together with the classical machine learning routines used as benchmark. ${\rm LiH}$ |  RMSE(E) |  # qubits |  depth | # of parameters
---|---|---|---|---
VQE | $1.4\cdot 10^{-4}\,{\rm Ha}$ | $3$ | 682 | 73
DNN | $1.5\cdot 10^{-3}\,{\rm Ha}$ | $\times$ | $\times$ | 73
${\rm H}_{2}{\rm O}$ |  RMSE(E) |  # qubits |  depth | # of parameters
VQE | $1.8\cdot 10^{-4}\,{\rm Ha}$ | $3$ | 818 | 87
DNN | $2.2\cdot 10^{-4}\,{\rm Ha}$ | $\times$ | $\times$ | 87
BPNN | $2.5\cdot 10^{-5}\,{\rm Ha}$ | $\times$ | $\times$ | 1642

In both cases, the QELM outperforms the VQE, not only in terms of RMSE but
also in terms of needed resources. In fact, even if the VQE routine of [19]
uses a smaller circuit ($3$ qubits), it needs a circuit with a depth
approximately $200$ times bigger than in the QELM counterpart. Moreover, the
VQE circuit needs to be run many more times, in order to obtain the gradients
of the parameters (see section A) and in order to perform more training
epochs.

In Table 2, we also included for completeness the comparison of the study
conducted by the authors of [19] using a fully connected multi-layered deep
neural network (DNN) with hyperbolic tangent activation function (five layers
with $[7,4,5,2,1]$ units for LiH and 6 layers with $[7,4,6,2,2,1]$ units for
H2O) and the Behler-Parrinello-like neural network (BPNN) of [42]. The number
of parameters to be optimized is taken to be the same for VQE and DNN for each
molecular species while it is greater for BPNN. As it is evident, also the VQE
can outperform the classical networks. The QELM not only improves the VQE
performances, but it does so while keeping minimal quantum computational cost.

In realistic implementations, the probability matrices $\boldsymbol{P}_{\rm
tr}$ and $\boldsymbol{P}_{\rm test}$ are estimated with finite precision, due
to the fact that probabilities are estimated by performing a limited number of
shots on a quantum device. This implies that the accuracy of the QELM
predictions is in general lower than the ideal case. The results of the
training either in the noiseless simulation or the real-hardware
implementation on IBM_BRISBANE are again collected in table Table 1, where
they can be compared with their statevector counterpart. One can see that
because of the finite statistics, the RMSE increases by a factor that depends
on the specific molecule. However, the RMSE in the QASM simulations and
IBM_BRISBANE implementation are always of the same order of magnitude, meaning
that the noise does not affect significantly the performances and the
noiseless simulations are reliable. The dependence of the RMSE on the number
of training elements follows the same argument as above. Moreover, looking at
Figure 6 it is evident that another effect of the finite statistics is that
the RMSEs for different numbers of qubits are much closer to each other on the
plateau.

In order to have a visual perspective of the accuracy of our predictions at
finite statistics (either in QASM simulation or IBM_BRISBANE realization), we
plot in Figure 7, the value of the force computed with classical methods
against the respective QELM prediction (in the ${\rm LiH}$ case we rather plot
the value of the forces in the $r-F$ plane).

Figure 5:  Behaviour of the energy RMSE in statevector simulations for a
different number of the encoding qubits and different sizes of the training
set. Figure 6:  Behaviour of the energy RMSE in QASM simulations and
IBM_BRISBANE realizations for different numbers of the encoding qubits and
different sizes of training set. Figure 7:  Comparison between classical
forces and QELM predictions obtained with QASM simulations and IBM_BRISBANE
realizations. In the ${\rm LiH}$ case, we explicitly plot the force field
curve.

##  4 Conclusions

In this work, we have demonstrated a fits application of QELM to quantum
chemistry. The QELM is trained to learn the functional dependence of the
potential energy surface and force fields of a molecule on a set of
generalized coordinates that parameterize its possible conformations. The
predictions of the QELM can be then exploited to speed up ab-initio molecular
dynamics computations that require the knowledge of the PES and FF at HS or
DFT precision.

We have studied in detail three case studies, ${\rm LiH}$, ${\rm H}_{2}{\rm
O}$ and ${\rm HCONH}_{2}$ but the proposed setup is scalable, in the sense
that can be applied to molecules with any number of atoms and degrees of
freedom. The setup is also optimized for practical implementation on quantum
devices, with IBM superconducting transmon quantum computers in mind. However,
by changing the set of native gates, the setup can be easily adapted to other
platforms, e.g. trapped-ion. The precision of the prediction can be
arbitrarily increased by utilizing more encoding qubits. The number of qubits
in turn fixes the dimension of the training set that should be used.

The accuracy of the QELM predictions is first tested in noiseless simulations
at infinite statistics, finding an exceptional agreement. We then investigate
the performances at finite statistics, either with noiseless simulation and
also with practical implementations on the IBM_BRISBANE quantum device. In
both cases we find encouraging results, suggesting that the proposed method
can be effectively competitive in a near future. It must be stressed that QELM
is particularly efficient in terms of quantum resource needed (it does not
require any optimization of quantum circuit parameters via gradient descent)
and in terms of classical post-processing (we do not perform any error
mitigation).

We leave to future work the application to more complicated molecules in order
to understand the ultimate frontiers of our proposal and the study of
efficient methods to overcome the limitations deriving from finite statistics.

We would like to thank Luca Innocenti for invaluable discussions. We
acknowledge the use of IBM Quantum Credits for this work. The views expressed
are those of the authors and do not reflect the official policy or position of
IBM or the IBM Quantum team. GLM and GMP acknowledge funding from the European
Union - NextGenerationEU through the Italian Ministry of University and
Research under PNRR - M4C2-I1.3 Project PE-00000019 ”HEAL ITALIA” (CUP
B73C22001250006 ). MB thanks “SiciliAn MicronanOTecH Research And innovation
CEnter - SAMOTHRACE” (MUR, PNRR-M4C2, ECS0000002), spoke 3, Università degli
Studi di Palermo, for financial support. SL and GMP acknowledge support by MUR
under PRIN Project No. 2022FEXLYB Quantum Reservoir Computing (QuReCo) and by
the “National Centre for HPC, Big Data and Quantum Computing (HPC)” Project
CN00000013 HyQELM – SPOKE 10.

## References

## References

  * [1] Car R and Parrinello M 1985 Physical review letters 55 2471
  * [2] Behler J and Parrinello M 2007 Physical review letters 98 146401
  * [3] Behler J 2011 The Journal of chemical physics 134
  * [4] Rupp M, Tkatchenko A, Müller K R and Von Lilienfeld O A 2012 Physical review letters 108 058301
  * [5] Behler J 2016 The Journal of chemical physics 145
  * [6] Gastegger M, Schwiedrzik L, Bittermann M, Berzsenyi F and Marquetand P 2018 The Journal of chemical physics 148
  * [7] Imbalzano G, Anelli A, Giofré D, Klees S, Behler J and Ceriotti M 2018 The Journal of chemical physics 148
  * [8] Cheng B, Mazzola G, Pickard C J and Ceriotti M 2020 Nature 585 217–220
  * [9] Unke O T, Chmiela S, Sauceda H E, Gastegger M, Poltavsky I, Schütt K T, Tkatchenko A and Müller K R 2021 Chemical Reviews 121 10142–10186
  * [10] Krems R 2019 Physical Chemistry Chemical Physics 21 13392–13410
  * [11] Asnaashari K and Krems R V 2021 Machine Learning: Science and Technology 3 015005
  * [12] Jasinski A, Montaner J, Forrey R, Yang B, Stancil P, Balakrishnan N, Dai J, Vargas-Hernández R and Krems R 2020 Physical Review Research 2 032051
  * [13] Cui J and Krems R V 2016 Journal of Physics B: Atomic, Molecular and Optical Physics 49 224001
  * [14] Dai J and Krems R V 2023 Machine Learning: Science and Technology 4 045027
  * [15] Blank T B, Brown S D, Calhoun A W and Doren D J 1995 The Journal of chemical physics 103 4129–4137
  * [16] Guan S H, Shang C and Liu Z P 2021 Chinese Journal of Chemistry 39 3144–3154
  * [17] Ollitrault P J, Miessen A and Tavernelli I 2021 Accounts of Chemical Research 54 4229–4238
  * [18] Sajjan M, Li J, Selvarajan R, Sureshbabu S H, Kale S S, Gupta R, Singh V and Kais S 2022 Chemical Society Reviews 51 6475–6573
  * [19] Kiss O, Tacchino F, Vallecorsa S and Tavernelli I 2022 Machine Learning: Science and Technology 3 035004
  * [20] Dai J and Krems R V 2022 The Journal of Chemical Physics 156
  * [21] Fujii K and Nakajima K 2017 Physical Review Applied 8 024030
  * [22] Chen J and Nurdin H I 2019 Quantum Information Processing 18 1–36
  * [23] Nakajima K, Fujii K, Negoro M, Mitarai K and Kitagawa M 2019 Physical Review Applied 11 034021
  * [24] Martínez-Peña R, Nokkala J, Giorgi G L, Zambrini R and Soriano M C 2020 Cognitive Computation 1–12
  * [25] Chen J, Nurdin H I and Yamamoto N 2020 Physical Review Applied 14 024065
  * [26] Kutvonen A, Fujii K and Sagawa T 2020 Scientific reports 10 14687
  * [27] Martínez-Peña R, Giorgi G L, Nokkala J, Soriano M C and Zambrini R 2021 Physical Review Letters 127 100502
  * [28] Nokkala J, Martínez-Peña R, Giorgi G L, Parigi V, Soriano M C and Zambrini R 2021 Communications Physics 4 53
  * [29] Innocenti L, Lorenzo S, Palmisano I, Ferraro A, Paternostro M and Palma G M 2023 Communications Physics 6 118
  * [30] Innocenti L, Lorenzo S, Palmisano I, Albarelli F, Ferraro A, Paternostro M and Palma G M 2023 arXiv preprint arXiv:2301.13229
  * [31] Domingo L, Carlo G and Borondo F 2022 arXiv preprint arXiv:2205.10107
  * [32] Thanasilp S, Wang S, Cerezo M and Holmes Z 2022 arXiv preprint arXiv:2208.11060
  * [33] Xiong W, Facelli G, Sahebi M, Agnel O, Chotibut T, Thanasilp S and Holmes Z 2023 arXiv preprint arXiv:2312.15124
  * [34] Schuld M and Petruccione F 2021 Quantum Models as Kernel Methods (Cham: Springer International Publishing) pp 217–245 ISBN 978-3-030-83098-4 URL https://doi.org/10.1007/978-3-030-83098-4_6
  * [35] Schuld M, Sweke R and Meyer J J 2021 Physical Review A 103 032430
  * [36] Jäger J and Krems R V 2023 Nature Communications 14 576 ISSN 2041-1723
  * [37] Jerbi S, Fiderer L J, Poulsen Nautrup H, Kübler J M, Briegel H J and Dunjko V 2023 Nature Communications 14 517 ISSN 2041-1723
  * [38] Frisch M J, Trucks G W, Schlegel H B, Scuseria G E, Robb M A, Cheeseman J R, Scalmani G, Barone V, Petersson G A, Nakatsuji H, Li X, Caricato M, Marenich A V, Bloino J, Janesko B G, Gomperts R, Mennucci B, Hratchian H P, Ortiz J V, Izmaylov A F, Sonnenberg J L, Williams-Young D, Ding F, Lipparini F, Egidi F, Goings J, Peng B, Petrone A, Henderson T, Ranasinghe D, Zakrzewski V G, Gao J, Rega N, Zheng G, Liang W, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Throssell K, Montgomery Jr J A, Peralta J E, Ogliaro F, Bearpark M J, Heyd J J, Brothers E N, Kudin K N, Staroverov V N, Keith T A, Kobayashi R, Normand J, Raghavachari K, Rendell A P, Burant J C, Iyengar S S, Tomasi J, Cossi M, Millam J M, Klene M, Adamo C, Cammi R, Ochterski J W, Martin R L, Morokuma K, Farkas O, Foresman J B and Fox D J 2016 Gaussian˜16 Revision C.01 gaussian Inc. Wallingford CT
  * [39] Zhao Y and Truhlar D G 2008 Theoretical chemistry accounts 120 215–241
  * [40] Qiskit contributors 2023 Qiskit: An open-source framework for quantum computing
  * [41] Higham N J 2002 Accuracy and stability of numerical algorithms (SIAM)
  * [42] Singraber A, Behler J and Dellago C 2019 Journal of chemical theory and computation 15 1827–1840
  * [43] Mitarai K, Negoro M, Kitagawa M and Fujii K 2018 Physical Review A 98 032309
  * [44] Schuld M, Bergholm V, Gogolin C, Izaac J and Killoran N 2019 Physical Review A 99 032331

##  Appendix A Force field generation with quantum routines

All supervised QML routines consist in evaluating the expectation values of
some operator $\Tr(\mathcal{O}(\boldsymbol{\theta},\boldsymbol{x})\rho_{0})$
where $\rho_{0}$ is some reference state. $\boldsymbol{\theta}$ is a
collection of parameters that in VQE routines gets optimized to fit the
training data. The optimization consists of a gradient descend algorithm
performed on a classical computer, while the gradients of the parameters are
estimated on a quantum computer through the parameter-shift rule [43, 44]:

$\partial_{\theta}f(\theta)\,=\,\,\lambda\left(f\left(\theta+\frac{\pi}{4\lambda}\right)-f\left(\theta-\frac{\pi}{4\lambda}\right)\right)\,,$
(21)

where $f(\theta)=\Tr(e^{-i\theta G}\mathcal{O}e^{i\theta G}\rho)$ and $G$ is a
Hermitian operator whose spectrum has only two distinct eigenvalues
$\pm\lambda$. When $G$ has more than two distinct eigenvalues, it is possible
to construct generalized shift rules as in (21) where the sum runs over more
than two shifts. Using the parameter-shift rule, it is also possible to
compute the derivative with respect to the input data. In the context of PES
predictions, where the VQE is used to predict conformational energies, these
derivatives coincide with the forces. This knowledge can be used to train a
unique circuit able to provide either energies and forces. However, in order
to get all this information, the circuit needs to be run
$2\mathcal{\chi}\mathcal{S}+1$ times, where $\mathcal{\chi}$ is the number of
coordinates and $\mathcal{S}$ is the number of shifts appearing in the
parameter-shift rule.

The parameter-shift rule can be also applied in the QELM setup but this does
not lead to great advantages. In fact, in the QELM routines, a unique encoding
of the coordinates $\boldsymbol{x}$ is enough to predict either energies and
forces: the unique difference is that a bigger matrix $\boldsymbol{W}$ is
obtained at the classical stage of the algorithm. This greatly reduces the
running time of a QELM simulation for force field predictions.


# Compositional Fairness Constraints for Graph Embeddings

Avishek Joey Bose  William L. Hamilton

###### Abstract

Learning high-quality node embeddings is a key building block for machine
learning models that operate on graph data, such as social networks and
recommender systems. However, existing graph embedding techniques are unable
to cope with fairness constraints, e.g., ensuring that the learned
representations do not correlate with certain attributes, such as age or
gender. Here, we introduce an adversarial framework to enforce fairness
constraints on graph embeddings. Our approach is compositional—meaning that it
can flexibly accommodate different combinations of fairness constraints during
inference. For instance, in the context of social recommendations, our
framework would allow one user to request that their recommendations are
invariant to both their age and gender, while also allowing another user to
request invariance to just their age. Experiments on standard knowledge graph
and recommender system benchmarks highlight the utility of our proposed
framework.

Machine Learning, ICML

##  1 Introduction

Learning low-dimensional embeddings of the nodes in a graph is a fundamental
technique underlying state-of-the-art approaches to link prediction and
recommender systems (Hamilton et al., 2017b). However, in many
applications—especially those involving social graphs—it is desirable to
exercise control over the information contained within learned node
embeddings. For instance, we may want to ensure that recommendations are fair
or balanced with respect to certain attributes (e.g., that they do not depend
on a user’s race or gender) or we may want to ensure privacy by not exposing
certain attributes through learned node representations. In this work we
investigate the feasibility of enforcing such invariance constraints on
(social) graph embeddings.

Figure 1: Overview of our approach: Our goal is to generate graph embeddings
that are invariant to particular sensitive attributes (e.g., age or gender).
We train a set of “filters” to prevent adversarial discriminators from
classifying the sensitive information from the filtered embeddings. After
training, these filters can be composed together in different combinations,
allowing the flexible generation of embeddings that are invariant w.r.t. any
subset of the sensitive attributes.

While enforcing invariance constraints on general classification models
(Chouldechova, 2017; Gajane & Pechenizkiy, 2017; Kamishima et al., 2012) and
collaborative filtering algorithms (Yao & Huang, 2017) has received
considerable attention in recent years, these techniques have yet to be
considered within the context of graph embeddings—a setting that introduces
particular challenges due to the non-i.i.d. and non-Euclidean nature of
relational, graph data.

Moreover, in the case of social graphs and large-scale recommender systems, it
is often the case that there are many possible sensitive attributes that we
may want to enforce invariance constraints over. Previous work on enforcing
invariance (or “fairness”) in social applications has generally focused on
situations that involve one sensitive attribute (e.g., age in the context of
credit or loan decisions; Zemel et al. (2013)), but in the context of social
graph embeddings there can be an extremely large number of possible sensitive
attributes. In fact, in extreme settings we may even want to be fair with
respect to the existence of individual edges. For instance, a user on a social
networking platform might want that platform’s recommender system to ignore
the fact that they are friends with a certain other user, or that they engaged
with a particular piece of content.

Our contributions. We introduce an adversarial framework to enforce
compositional fairness constraints on graph embeddings for multiple sensitive
attributes. The key idea behind our approach is that we learn a set of
adversarial filters that remove information about particular sensitive
attributes. Crucially, each of these learned filters can be optionally applied
after training, so the model can flexibly generate embeddings that are
invariant with respect to different combinations of sensitive attributes. As
the space of possible combinations of sensitive attributes can be
combinatorially large, we demonstrate that our compositional strategy can
generate invariant embeddings even on unseen combinations at test time. Our
contribution is at the intersection of research on (social) graph embedding
and algorithmic fairness. We build upon the success of recent adversarial
approaches to fairness (Edwards & Storkey, 2015), disentanglement (Mathieu et
al., 2016), and transfer learning (Madras et al., 2018)–extending these
approaches to the domain of graph representation learning and introducing new
algorithmic techniques to accommodate compositional constraints during
inference.

##  2 Related Work

We now briefly highlight core related work on (social) graph embeddings and
algorithmic fairness, which our research builds upon.

###  2.1 Graph Embedding

At the core of our proposed methodology is the notion of learning low-
dimensional embeddings of graph-structured data, especially social data. Graph
embedding techniques have a long history in the social sciences, with
connections to early research on “sociograms” (small hand-constructed social
networks) and latent variable models of social interactions (Faust, 1988;
Majone, 1972). In more recent years, the task of embedding graph-structured
data has received increasing attention from the machine learning and data
mining communities (Cai et al., 2018; Hamilton et al., 2017b). Generally, the
goal of these works is to map graph nodes to low-dimensional vector
embeddings, such that the original graph can be reconstructed from these
embeddings. Traditional approaches to this problem include Laplacian eigenmaps
(Belkin & Niyogi, 2002) and matrix factorization techniques (Ng et al., 2001),
with recent years witnessing a surge in methods that rely on random-walk based
objectives (Grover & Leskovec, 2016; Perozzi et al., 2014), deep autoencoders
(Wang et al., 2016), and graph neural networks (Hamilton et al., 2017a; Kipf &
Welling, 2016).

Learned graph embeddings can be used for a wide variety of tasks, including
node classification, relation prediction, and clustering (Hamilton et al.,
2017b). Here, we focus on the relation prediction task, i.e., using the
learned representations to predict previously unobserved relationships between
the input nodes. The relation prediction task is exceptionally general—for
example, it generalizes basic recommender systems, knowledge base completion,
and even node classification (ibid.).

###  2.2 Algorithmic Fairness

Unlike previous research on graph embedding, in this work we focus on the
challenge of enforcing fairnes or invariance constraints on the learned
representations. Recent work on fairness in machine learning, including work
on fairness in collaborative filtering, involves making predictions that are
balanced or invariant with respect to certain sensitive variables (e.g., age
or gender) (Chouldechova, 2017; Gajane & Pechenizkiy, 2017; Kamishima et al.,
2012; Madras et al., 2018; Zemel et al., 2013; Yao & Huang, 2017). Formally,
in the standard “fair classification” setting we consider a data point
$\mathbf{x}\in\mathbb{R}^{n}$, its class label $y\in\mathcal{Y}$, and a binary
sensitive attribute $a\in\\{0,1\\}$ (e.g., indicating gender). The high-level
goal is then to train a model to predict $y$ from $\mathbf{x}$, while making
this prediction invariant or fair with respect to $a$ (Madras et al., 2018).
There are many specific definitions of fairness, such as whether fairness
refers to parity or satisfying certain preferences (see (Gajane & Pechenizkiy,
2017) for a detailed discussion). In the context of fair machine learning, our
core contribution is motivating, implementing, and evaluating an approach to
enforce fairness within the context of graph embeddings. There are a number of
complications introduced by this setting—for instance, rather than one
classification task, we instead have thousands or even millions of
interdependent edge relationships.

##  3 Preliminaries

We consider the general case of embedding a heterogeneous or multi-relational
(social) graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, which consists of a
set of directed edge triples $e=\langle u,r,v\rangle\in\mathcal{E}$, where
$u,v\in\mathcal{V}$ are nodes and $r\in\mathcal{R}$ is a relation type. We
further assume that each node is of a particular type,
$\mathcal{T}\subseteq\mathcal{V}$, and that relations may have constraints
regarding the types of nodes that they can connect.

Relation Prediction. The general relation prediction task on such a graph is
as follows. Let $\mathcal{E}_{\textrm{train}}\subset\mathcal{E}$ denote a set
of observed training edges and let $\bar{\mathcal{E}}=\\{\langle
v_{i},r,v_{j}\rangle:v_{i},v_{j}\in\mathcal{V},r\in R\\}\setminus\mathcal{E}$
denote the set of negative edges that are not present in the true graph
$\mathcal{G}$. Given $\mathcal{E}_{\textrm{train}}$, we aim to learn a scoring
function $s$ such that

$s(e)>s(e^{\prime}),\forall e\in\mathcal{E},e^{\prime}\in\bar{\mathcal{E}}.$
(1)

In other words, the learned scoring function should ideally score any true
edge higher than any negative edge.

Embedding-based Models. In the context of graph embeddings, we aim to solve
this relation prediction task by learning a function
$\textsc{enc}:\mathcal{V}\mapsto\mathbb{R}^{d}$ that maps each node
$v\in\mathcal{V}$ to an embedding $\mathbf{z}_{v}=\textsc{enc}(v)$. In this
case, the signature of the score function becomes
$s:\mathbb{R}^{d}\times\mathcal{R}\times\mathbb{R}^{d}\mapsto\mathbb{R}$,
i.e., it takes two node embeddings
$\mathbf{z}_{u},\mathbf{z}_{v}\in\mathbb{R}^{d}$ and a relation
$r\in\mathcal{R}$ and scores the likelihood that the edge $e=<u,r,v>$ exists
in the graph. Generally, the intuition in embedding-based approaches is that
the distance between two node embeddings should encode the likelihood that
there is an edge between the nodes. Following standard practice, we consider
the optimization of these scoring functions using contrastive learning methods
that make use of a corruption distribution such as noise contrastive
estimation (Dyer, 2014; Mnih & Teh, 2012) and similar variants (Bose et al.,
2018), where the loss over a batch of edges
$\mathcal{E}_{\textrm{batch}}\subseteq\mathcal{E}_{\textrm{train}}$ is given
by:

$\sum_{e\in\mathcal{E}_{\textrm{edge}}}L_{\textrm{edge}}(s(e),s(e^{-}_{1}),...,s(e^{-}_{m})),$
(2)

where $L_{\textrm{edge}}$ is a per-edge loss function and
$e^{-}_{1},...,e^{-}_{m}\in\bar{\mathcal{E}}$ are “negative samples”, i.e.,
randomly sampled edges that do not exist in the graph. Loss functions of this
form generally attempt to maximize the likelihood of true edges compared to
the negative samples.

Fairness. In order to incorporate the notion of fairness into the graph
embedding setup, we assume that for exactly one node type $\mathcal{T}^{*}$,
all nodes of this type, i.e. all $u\in\mathcal{T}^{*}$, have $K$ categorical
sensitive attributes, $a^{k}_{u}\in\mathcal{A}_{k},k=1...,K$, and for
simplicity, we assume that there are no other features or attributes
associated with the nodes and edges in the graph.111Though this assumption can
easily be relaxed. The challenge in enforcing fairness is thus to ensure that
the learned node embeddings, $\mathbf{z}_{u}$, are not biased or unfair with
respect to these sensitive attributes—a point which we formalize in the next
section.

##  4 Invariant Graph Embeddings

We first motivate and argue in favor of a particular form of “fairness” (or
rather invariance) within the context of graph embeddings. Following this, we
outline our compositional and adversarial approach for enforcing these
invariance constraints on graph embeddings.

###  4.1 Pragmatic Fairness as Invariance

In this paper we consider a simple, user-centric formulation of fairness
within the context of social graph embeddings. Using gender as an example of a
sensitive attribute and movie recommendation as an example relation prediction
task, our approach is guided by the following question: If one gives a user a
button that says “Please ignore my gender when recommending movies”, what does
a user expect from the system after this button is pressed? Here, we accept it
as non-controversial that the expectation from the user is that recommendation
does not depend in any way on their gender, i.e., that the recommendation
would be the same regardless of their gender. Formally, given a user $u$, this
expectation amounts to an assumption of independence,

$s(e)\perp a_{u}\qquad\forall v\in\mathcal{V},r\in\mathcal{R}$ (3)

between the recommendation—i.e., the score of the edge,
$s(e)=s(\langle\mathbf{z}_{u},r,\mathbf{z}_{v}\rangle)$—and the sensitive
attribute $a_{u}$.

One issue in directly enforcing Equation (3) is that there are many
(potentially millions) of possible edges that we might want to score for every
node $u\in\mathcal{T}^{*}$, making it intractable to enforce independence on
each of these decisions individually. However, if we assume that the score
function $s(\langle\mathbf{z}_{u},r,\mathbf{z}_{v}\rangle)$ depends on $u$
only through $u^{\prime}s$’s embedding, $\mathbf{z}_{u}$, then we can
guarantee the independence in Equation (3) for all edge predictions by
enforcing what we call representational invariance:

$\mathbf{z}_{u}\perp a_{u},\>\>\forall u\in\mathcal{V}.$ (4)

In other words, we require that the mutual information
$I(\mathbf{z}_{u},a_{u})$ is $0$.

Generalizing to the setting of multiple sensitive attributes, for a given set
of sensitive attributes $S\subseteq\\{1,...,K\\}$, we would require that

$I(\mathbf{z}_{u},a^{k}_{u})=0,k\in S,\forall u\in\mathcal{V},$ (5)

which amounts to the assumption of $S$ independent invariance constraints on
the $S$ distinct sensitive attributes.222Note that this does not necessarily
imply “subgroup fairness” on the joint distribution (Kearns et al., 2017).
Importantly, we assume that the set $S$ is not fixed (e.g., different users
might request different invariance constraints). In the language of
algorithmic fairness, the representational invariance of Equation (5) implies
that traditional demographic parity constraints are satisfied on the sensitive
attributes and recommendations.

###  4.2 Model Definition

In this work, we enforce representational invariance constraints on the node
embeddings (Equation 5) by introducing an adversarial loss and a technique to
“filter” the embeddings generated by the enc function. Note, again, that a
unique challenge here is that $S$—the set of sensitive attributes we want to
be invariant with respect to—is not fixed across nodes; i.e., we may want to
enforce invariance on different sets of sensitive attributes for different
nodes.

Note also that the framework presented in this section is quite general and
can function with arbitrary combinations of base node embedding functions enc
and edge-prediction losses $L_{\textrm{edge}}$ (see Equation 2). We discuss
three concrete instantiations of this framework in Section 5.

Compositional Encoder. The first key component of our model is generalizing
the enc embedding function to optionally “filter” out the information about
certain sensitive attributes. In particular, for every sensitive attribute
$k\in\\{1,...,K\\}$ we define a filter function
$f_{k}:\mathcal{R}^{d}\mapsto\mathcal{R}^{d}$ that is trained to remove the
information about the $k$th sensitive attribute. If we want the node embedding
to be invariant w.r.t. some set of sensitive attributes
$S\subseteq\\{1,...,K\\}$, we then generate its embedding by “composing” the
output of the $|S|$ filtered embeddings using a compositional encoder:

$\textsc{c-enc}(u,S)=\frac{1}{|S|}\sum_{k\in S}f_{k}(\textsc{enc}(u))$ (6)

To train $\textsc{c-enc}(u,S)$ we sample a binary mask to determine the set
$S$ at every iteration. In this work, we sample the binary mask as a sequence
of $k$ independent Bernoulli draws with a common fixed probability $p=0.5$;
however, other application-specific distributions (e.g., incorporating
dependencies between the attributes) could be employed. Sampling random binary
masks forces the model to produce invariant embeddings for different
combinations of sensitive attributes during training with the hope of
generalizing to unseen combinations during inference time — a phenomena that
we empirically validate in Section 5.2.

Adversarial Loss. To train the compositional encoder, we employ an adversarial
regularizer. For each sensitive attribute $k\in K$, we define a discriminator
$D_{k}:\mathbb{R}^{d}\times\mathcal{A}_{k}\mapsto[0,1]$, which attempts to
predict the $k$th sensitive attribute from the node embeddings. Assuming we
are given an edge-prediction loss function $L_{\textrm{edge}}$ (as in Equation
2), we can then define our new adversarially regularized per-edge loss as

$\displaystyle L(e)=$ $\displaystyle
L_{\textrm{edge}}(s(e),s(e^{-}_{1}),...,s(e^{-}_{m}))$
$\displaystyle+\lambda\sum_{k\in
S}\sum_{a^{k}\in\mathcal{A}_{k}}\log(D_{k}(\textsc{c-enc}(u,S),a^{k})),$ (7)

where $\lambda$ is a hyperparameter controlling the strength of the
adversarial regularization. To optimize this loss in a minibatch setting, we
alternate between two types of stochastic gradient descent updates: (1) $T$
minibatch updates minimizing $L(e)$ with respect to c-enc (with all the
$D_{k}$ fixed), and (2) $T^{\prime}$ minibatch updates minimizing $-L(e)$ with
respect to $D_{k},k=1...,K$ (with c-enc fixed).

Theoretical Considerations. For clarity and simplicity, we consider the case
of a single binary sensitive attribute, with the theoretical intuitions
naturally generalizing to the multi-attribute and multi-class settings.
Assuming a single binary sensitive attribute $a_{k}$, by simple application of
Proposition 2 in Goodfellow et al. (2014), we have:333Theorem 1 holds as a
consequence of Proposition 2 in Goodfellow et al. (2014) if we simply replace
the task of distinguishing real/fake data by classifying a binary sensitive
attribute.

######  Theorem 1.

If c-enc and $D_{k}$ have enough capacity, $T^{\prime}$ is large enough so
that $D_{k}$ is allowed to reach its optimum on $-L(e)$ (with c-enc fixed),
and c-enc is optimized according to $L(e)$ (with $D$ fixed), then
$I(\mathbf{z}_{u},a_{u})\rightarrow 0,\forall u\in\mathcal{T}^{*}$ as
$\lambda\rightarrow\infty$.

That is, if we increase the weight of the adversarial regularizer to infinity,
the equilibrium of the minimax game in Equation (4.2) occurs when there is
zero mutual information between the sensitive attribute and the embeddings. Of
course, as $\lambda\rightarrow\infty$ trivial solutions to this game exist
(e.g., c-enc simply outputting a constant value) and in practice setting
$\lambda<\infty$ leads to a tradeoff between performance on edge prediction
and representational invariance.

##  5 Experiments

We investigated the impact of enforcing invariance on graph embeddings using
three datasets: Freebase15k-237444www.microsoft.com/en-
us/download/details.aspx?id=52312,
MovieLens-1M555grouplens.org/datasets/movielens/1m/, and an edge-prediction
dataset derived from Reddit.666reddit.com The dataset statistics are given in
Table 5.1. Our experimental setup closely mirrors that of (Madras et al.,
2018) where we jointly train the main model with adversaries, but when testing
invariance, we train a new classifier (with the same capacity as the
discriminator) to predict the senstive attributes from the learned embeddings.

The goal of our experiments was to answer three questions:

  1. (Q1)

The cost of invariance. What is the tradeoff between enforcing invariance on
the representations and accuracy on the main edge prediction task?

  2. (Q2)

The impact of compositionality. How does the performance of a compositional
approach, which jointly enforces fairness over a set of sensitive attributes,
compare to a more traditional model that only enforces fairness on a single
attribute?

  3. (Q3)

Invariance on unseen combinations. In settings with many sensitive attributes,
is our approach able to enforce invariance even on combinations of sensitive
attributes that it never saw during training?

Throughout these experiments, we rely on two baselines: First, we compare
against baselines that do not include any invariance constraints, i.e., models
with $\lambda=0$. Second, we compare against a non-compositional adversarial
approach where we separately train $K$ distinct encoders and $K$ distinct
adversaries for each of the $K$ sensitive attributes in the data. This non-
compositional adversary is essentially an extension of Edwards & Storkey
(2015)’s approach to the graph embedding domain.

###  5.1 Setup and Datasets

Table 1: Statistics for the three datasets, including the total number of
nodes $(|\mathcal{V}|$) and number of nodes with sensitive attributes
$|\mathcal{T}^{*}|$, the number of sensitive attributes and their types and
the total number of edges in the graph.

Dataset | $|\mathcal{V}|$ | $|\mathcal{T}^{*}|$ |  #Sensitive |  |  |
---|---|---|---|---|---|---
| Edges |  Binary | Attributes |  |  |
|  Multiclass | Attributes? |  |  |  |
| Attributes? |  |  |  |  |
FB15k-237 | 14,940 | 14,940 | 3 | 168,618 | $\surd$ | $\times$
MovieLens1M | 9,940 | 6,040 | 3 | 1,000,209 | $\surd$ | $\surd$
Reddit Comments | 385,735 | 366,797 | 10 | 7,255,096 | $\surd$ | $\times$

Before describing our experimental results, we first outline the key
properties of the datasets we used, as well as the specific encoders and edge-
prediction loss functions used.

In all experiments, we used multi-layer perceptrons (MLPs) with leaky ReLU
activation functions (Xu et al., 2015) as the discriminators $D_{k}$ and
filters $f_{k}$. The Appendix contains details on the exact hyperparameters
(e.g., number of layers and sizes) used for all the different experiments, as
well as details on the training procedures (e.g., number of epochs and data
splits). Code to reproduce our results is available at:
https://github.com/joeybose/Flexible-Fairness-Constraints.

#### Freebase15k-237

Freebase 15k-237 is a standard benchmark used for knowledge base completion
(Toutanova et al., 2015). In this work, we use Freebase 15k-237 as a semi-
synthetic testbed to evaluate the impact of adversarial regularization. Taking
the entity attribute labels from Moon et al. (2017), we used the $3$-most
common attribute labels (e.g., /award/award_nominee) as “sensitive”
attributes. The goal in this dataset is to perform the standard knowledge base
completion task, while having the entity embeddings be invariant with respect
to these “sensitive” attribute labels. While synthetic, this dataset provides
a useful reference point due to its popularity in the graph embedding
literature.

For our encoder and edge-prediction loss function, we follow Ji et al.
(2015)’s TransD approach, since we found this approach gave significant
performance boosts compared to simpler models (e.g., TransE). In this model,
the encoding of a node/entity depends on the edge relation being predicted, as
well as on whether the entity is the head or tail in a relation (i.e., the
edge direction matters). In particular, the embedding of the head node (i.e.,
the source node) in an edge relation is given by:

$\textsc{enc}(u,\langle
u,r,v\rangle)=(\mathbf{r}_{p}\mathbf{u}_{p}^{\top}+\mathbf{I}^{d\times
d})\mathbf{u},$ (8)

where $\mathbf{u},\mathbf{u}_{p},\mathbf{r}_{p}\in\mathbb{R}^{d}$ are
trainable embedding parameters and $\mathbf{I}^{d\times d}$ is a
$d$-dimensional identity matrix. The encoding function for the tail node is
defined analogously. The score function for this approach is given by

$\displaystyle s(\langle u,r,v\rangle)=-\|\textsc{enc}(u,\langle
u,r,v\rangle)+$ $\displaystyle\mathbf{r}$
$\displaystyle-\textsc{enc}(v,\langle u,r,$ $\displaystyle v\rangle)\|_{2},$

where $\mathbf{r}\in\mathbb{R}^{d}$ is another trainable embedding parameter
(one per relation). Finally, we use a standard max-margin loss with a single
negative sample per positive edge:

$L_{\textrm{edge}}(s(e),s(e^{-}))=\max(0,1-s(e)+s(e)^{-}).$ (9)

#### Movielens-1M

Our second dataset is derived from the MovieLens-1M recommender system
benchmark (Harper & Konstan, 2016). This is a standard recommender system
benchmark, where the goal is to predict the rating that users assign movies.
However, unlike previous work, in our experiments we treat the user features
(age, gender, and occupation) as sensitive attributes (rather than as
additional feature information for the recommendation task). Following Berg et
al. (2017) we treat this recommendation task as an edge prediction problem
between users and movies, viewing the different possible ratings as different
edge relations.

For this dataset we use a simple “embedding-lookup” encoder, where each user
and movie is associated with a unique embedding vector in $\mathbb{R}^{d}$. As
a scoring function, we follow Berg et al. (2017) and use a log-likelihood
approach:

$s(\langle
u,r,m\rangle)=\mathbf{z}_{u}^{\top}\mathbf{Q}_{r}\mathbf{z}_{v}-\log(\sum_{r^{\prime}\in\mathcal{R}}\mathbf{z}_{u}^{\top}\mathbf{Q}_{r^{\prime}}\mathbf{z}_{v}),$

The relation matrices $\mathbf{Q}_{r}\in\mathbb{R}^{d\times d}$ are computed
as:

$\mathbf{Q}_{r}=a_{r,1}\mathbf{P}_{1}+a_{r,2}\mathbf{P}_{2},$

where $a_{r,1},a_{r,1}\in\mathcal{R}$ and
$\mathbf{P}_{1},\mathbf{P}_{2}\in\mathcal{R}^{d\times d}$ are trainable
parameters. In this case, the loss function is simply the negative of the log-
likelihood score.

#### Reddit

The final dataset we consider is based on the social media website Reddit—a
popular, discussion-based website where users can post and comment on content
in different topical communities, called “subreddits”. For this dataset, we
consider a traditional edge prediction task, where the goal is to predict
interactions between users and subreddit communities.

To construct the edge prediction task, we examined all comments from the month
of November in $2017$, and we placed an edge between a user and a community if
this user commented on that community at least once within this time period.
We then took the $10$-core of this graph to remove low-degree nodes, which
resulted in a graph with approximately $366$K users, $18$K communities, and
$7$M edges. Given this graph, the main task is to train an edge-prediction
model on $90\%$ of the user-subreddit edges and then predict missing edges in
a held-out test set of the remaining edges.

Reddit is a pseudonymous website with no public user attributes. Thus, to
define sensitive attributes, we treat certain subreddit nodes as sensitive
nodes, and the sensitive attributes for users are whether or not they have an
edge connecting to these sensitive nodes. In other words, the fairness
objective in this setting is to force the model to be invariant to whether or
not a user commented on a particular community. To select the “sensitive”
subreddit communities, we randomly sampled $10$ from the top-100 communities
by degree.777We excluded the top-5 highest-degree outlying communities. Note
that this setting represents the extreme case where we want the model to be
invariant with respect to the existence of particular edges in the input
graph.

As with MovieLens-1M, we use a simple “embedding-lookup” encoder. In this
case, there is only a single relation type—indicating whether a Reddit user
has commented on a “subreddit” community. Thus, we employ a simple dot-product
based scoring function, $s(\langle
u,r,v\rangle)=\mathbf{z}_{u}^{\top}\mathbf{z}_{v}$, and we use a max-margin
loss as in Equation (9).

###  5.2 Results

We now address the core experimental questions (Q1-Q3).

Table 2: Ability to predict sensitive attributes on the MovieLens data when using various embedding approaches. For gender attribute the score is AUC while for age and occupation attributes the score is micro averaged F1. The columns represent the different embedding approaches (e.g., with or without adversarial regularizatin) while the rows are the attribute being classified. MovieLens1M |  Baseline |  |  |  |  |  |
---|---|---|---|---|---|---|---
|  Gender |  |  |  |  |  |
|  Age |  |  |  |  |  |
|  Occupation |  |  |  |  |  |
|  Comp. |  |  |  |  |  |
|  Majority |  |  |  |  |  |
|  Random |  |  |  |  |  |
| Classifier |  |  |  |  |  |
Gender | 0.712 | 0.532 | 0.541 | 0.551 | 0.511 | 0.5 | 0.5
Age | 0.412 | 0.341 | 0.333 | 0.321 | 0.313 | 0.367 | 0.141
Occupation | 0.146 | 0.141 | 0.108 | 0.131 | 0.121 | 0.126 | 0.05
Table 3: Ability to predict sensitive attributes on the Freebase15k-237 data
when using various embedding approaches. AUC scores are reported, since all
the sensitive attributes are binary. The mean rank on the main edge-prediction
task is also reported.

FB15k-237 |  Baseline |  |
---|---|---|---
|  Non |  |
|  Comp. |  |
| Adversary |  |
Attribute 0 | 0.97 | 0.82 | 0.77
Attribute 1 | 0.99 | 0.81 | 0.79
Attribute 2 | 0.98 | 0.81 | 0.81
Mean Rank | 285 | 320 | 542
Figure 2: Performance on the edge prediction (i.e., recommendation) task on
MovieLens, using RMSE as in Berg et al. (2017). Figure 3: Performance on the
edge prediction (i.e., recommendation) task on the Reddit data. Evaluation is
using the AUC score, since there is only one edge/relation type. Figure 4:
Ability to predict sensitive attributes on the Reddit data when using various
embedding approaches. Bar plots correspond to the average AUC across the $10$
binary sensitive attributes.

#### Q1: The Cost of Invariance

In order to quantify the extent to which the learned embeddings are invariant
to the sensitive attributes (e.g., after adversarial training), we freeze the
trained compositional encoder c-enc and train an new MLP classifier to predict
each sensitive attribute from the filtered embeddings (i.e., we train one new
classifier per sensitive attribute). We also evaluate the performance of these
filtered embeddings on the original prediction tasks. In the best case, a
newly trained MLP classifier should have random accuracy when attempting to
predict the sensitive attributes from the filtered embeddings, but these
embeddings should still provide strong performance on the main edge prediction
task. Thus, for binary sensitive attributes, an ideal result is an AUC score
of $0.5$ when attempting to predict the sensitive attributes from the learned
embeddings.

Figure 5: Tradeoff of Gender AUC score for a compositional adversary versus
different $\lambda$

Overall, we found that on the more realistic social recommendation
datasets—i.e., the MovieLens-1M and Reddit datasets—our approach was able to
achieve a reasonable tradeoff, with the near-complete removal of the sensitive
information leading to a roughly 10% relative error increase on the edge
prediction tasks. In other words, on these two datasets the sensitive
attributes were nearly impossible to predict from the filtered embeddings,
while the accuracy on the main edge prediction task was roughly 10% worse than
a baseline approach that does not include the invariance constraints. Table
5.2 and Figure 2 summarize these results for the MovieLens data, where we can
see that the accuracy of classifying the sensitive attributes is on-par with a
majority-vote classifier (Table 5.2) while the RMSE degrades from $0.865$ to
$1.01$ with the compositional adversary. Figures 5 and 6 illustrate this
tradeoff and show how the RMSE for the edge prediction task and ability to
predict the sensitive attributes change as we vary the regularization
strength, $\lambda$. As expected, increasing $\lambda$ does indeed produce
more invariant embeddings but at the cost of higher RMSE values. Figures 3 and
4 similiarly summarize these results on Reddit.

Interestingly, we found that on the Freebase15k-237 dataset it was not
possible to completely remove the sensitive information without incurring a
significant cost on the original edge prediction task. This result is not
entirely surprising, since for this dataset the “sensitive” attributes were
synthetically constructed from entity type annotations, which are presumably
very relevant to the main edge/relation prediction task. However, it is an
interesting point of reference that demonstrates the potential limitations of
removing sensitive information from learned graph embeddings.

Figure 6: Tradeoff of RMSE for a compositional adversary versus different
$\lambda$

#### Q2: The Impact of Compositionality

In all our experiments, we observed that our compositional approach performed
favorably compared to an approach that individually enforced fairness on each
individual attribute. In fact, on the MovieLens-1M data (and the synthetic
Freebase15k-237 data), the compostionally trained adversary outperformed the
individually trained adversaries in terms of removing information about the
sensitive attributes (Table 5.2). In other words, training a model to jointly
remove information about the sensitive attributes using the compositional
encoder (Equation 6) removed more information about the sensitive attributes
than training separate adversarially regularized embedding models for each
sensitive attribute. This result is not entirely surprising, as it essentially
indicates that the different sensitive attributes (age, gender, and
occupation) are correlated in this dataset. Nonetheless, it is a positive
result indicating that the extra flexibility afforded by the compositional
approach does not necessarily lead to a decrease in performance. That said, on
the Reddit data we observed the opposite trend and found that the
compositional approach performed worse in terms of its ability to remove
information about the sensitive attributes (Figure 4) as well as a small drop
on the performance of the main edge prediction task (Figure 3).

#### Q3: Invariance on Unseen Combinations

One of the key benefits of the compositional encoder is that it can flexibly
generate embeddings that are invariant to any subset of
$S\subseteq\\{1,...,K\\}$ of the sensitive attributes in a domain. In other
words, at inference time, it is possible to generate $2^{K}$ distinct
embeddings for an individual node, depending on the exact set of invariance
constraints. However, given this combinatorially large output space, a natural
question is whether this approach performs well when generalizing to unseen
combinations of sensitive attributes.

We tested this phenomenon on the Reddit dataset, since it has the largest
number of sensitive attributes (10, compared to 3 sensitive attributes for the
other two datasets). During training we held out $10\%$ of the combinations of
sensitive attributes, and we then evaluated the model’s ability to enforce
invariance on this held-out set. As we can see in Figure 4, the performance
drop for the held-out combinations is very small ($0.025$), indicating that
our compositional approach is capable of effectively generalizing to unseen
combinations. The Appendix contains further results demonstrating how this
trends scales gracefully when we increase the number of sensitive attributes
from 10 to 50.

#### Quantifying Bias

In all of the above results, we used the ability to classify the sensitive
attributes as a proxy for bias being contained within the embeddings. While
this is a standard approach, e.g., see Edwards & Storkey (2015), and an
intuitive method for evaluating representational invariance—a natural question
is whether the adversarial regularization also decreases bias in the edge
prediction tasks. Ideally, after filtering the embeddings, we would have that
the edge predictions themselves are not biased according to the sensitive
attributes.

To quantify this issue, we computed a “prediction bias” score for the
MovieLens1M dataset: For each movie, we computed the absolute difference
between the average rating predicted for each possible value of a sensitive
attribute and we then averaged these scores over all movies. Thus, for
example, the bias score for gender corresponds to the average absolute
difference in predicted ratings for male vs. female users, across all movies.
From the perspective of fairness our adversary imposes a soft demographic
parity constraint on the main task. A reduction in prediction bias across the
different subgroups represents an empirical measure of achieving demographic
parity. Figure 7 highlights these results, which show that adversarial
regularization does indeed drastically reduce prediction bias. Interestingly,
using a compositional adversary works better than a single adversary for a
specific sensitive attribute which we hypothesize is due to correlation
between sensitive attributes.

Figure 7: Prediction Bias for different Sensitive Attributes under three
settings in MovieLens1M.

##  6 Discussion and Conclusion

Our work sheds light on how fairness can be enforced in graph representation
learning—a setting that is highly relevant to large-scale social
recommendation and networking platforms. We found that using our proposed
compositional adversary allows us to flexibly accomodate unseen combinations
of fairness constraints without explicitly training on them. Crucially, this
highlights how fairness could be deployed in a real-word, user-driven setting,
where it is necessary to optionally enforce a large number of possible
invariance constraints over learned graph representations.

In terms of limitations and directions for future work, one important
limitation is that we only consider one type of adversarial loss to enforce
fairness. While this adversarial loss is theoretically motivated and known to
perform well, there are other recent variations in the literature (e.g.,
Madras et al. (2018))—as well as related non-adversarial regularizers (e.g.,
Zemel et al. (2013)). Also, while we considered imposing fairness over sets of
attributes, we did not explicitly model subgroup-level fairness (Kearns et
al., 2017). Extending and testing our framework with these alternatives is a
natural direction for future work.

There are also important questions about how our framework translates to real-
world production systems. For instance, in this work we enforced fairness with
respect to randomly sampled sets of attributes, but in real-world
environments, these sets of attributes would come from user preferences, which
may themselves be biased; e.g., it might be more common for female users to
request fairness than male users potentially leading to new kinds of
demographic inequalities. Understanding how these preference biases could
impact our framework is an important direction for future inquiry.

### Acknowledgements

The authors would like to thank the anonymous ICML reviewers for their helpful
comments. In addition we would like to thank Koustuv Sinha, Riashat Islam and
Andre Cianflone for helpful feedback on earlier drafts of this work. This
research was funded in part by an academic grant from Microsoft Research, as
well as a Canada CIFAR Chair in AI, held by Prof. Hamilton.

## References

  * Belkin & Niyogi (2002) Belkin, M. and Niyogi, P.  Laplacian eigenmaps and spectral techniques for embedding and clustering.  In _NIPS_ , 2002.
  * Berg et al. (2017) Berg, R. v. d., Kipf, T. N., and Welling, M.  Graph convolutional matrix completion.  _arXiv preprint arXiv:1706.02263_ , 2017.
  * Bose et al. (2018) Bose, A. J., Ling, H., and Cao, Y.  Adversarial contrastive estimation.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Long Papers)_ , 2018.  URL https://arxiv.org/abs/1805.03642.
  * Cai et al. (2018) Cai, H., Zheng, V. W., and Chang, K.  A comprehensive survey of graph embedding: problems, techniques and applications.  _IEEE Transactions on Knowledge and Data Engineering_ , 2018.
  * Chouldechova (2017) Chouldechova, A.  Fair prediction with disparate impact: A study of bias in recidivism prediction instruments.  _Big data_ , 5(2):153–163, 2017.
  * Dyer (2014) Dyer, C.  Notes on noise contrastive estimation and negative sampling.  _arXiv preprint arXiv:1410.8251_ , 2014.
  * Edwards & Storkey (2015) Edwards, H. and Storkey, A.  Censoring representations with an adversary.  In _ICLR_ , 2015.
  * Faust (1988) Faust, K.  Comparison of methods for positional analysis: Structural and general equivalences.  _Soc. Net._ , 10(4):313–341, 1988.
  * Gajane & Pechenizkiy (2017) Gajane, P. and Pechenizkiy, M.  On formalizing fairness in prediction with machine learning.  _arXiv preprint arXiv:1710.03184_ , 2017.
  * Goodfellow et al. (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y.  Generative adversarial nets.  In _NIPS_ , pp. 2672–2680, 2014.
  * Grover & Leskovec (2016) Grover, A. and Leskovec, J.  node2vec: Scalable feature learning for networks.  In _KDD_ , 2016.
  * Hamilton et al. (2017a) Hamilton, W., Ying, R., and Leskovec, J.  Inductive representation learning on large graphs.  In _NIPS_ , 2017a.
  * Hamilton et al. (2017b) Hamilton, W., Ying, R., and Leskovec, J.  Representation learning on graphs: Methods and applications.  _IEEE Data Engineering Bulletin_ , 2017b.
  * Harper & Konstan (2016) Harper, F. M. and Konstan, J. A.  The movielens datasets: History and context.  _Acm transactions on interactive intelligent systems (tiis)_ , 5(4):19, 2016.
  * Ji et al. (2015) Ji, G., He, S., Xu, L., Liu, K., and Zhao, J.  Knowledge graph embedding via dynamic mapping matrix.  In _ACL_ , 2015.
  * Kamishima et al. (2012) Kamishima, T., Akaho, S., Asoh, H., and Sakuma, J.  Fairness-aware classifier with prejudice remover regularizer.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , pp. 35–50. Springer, 2012.
  * Kearns et al. (2017) Kearns, M., Neel, S., Roth, A., and Wu, Z.  Preventing fairness gerrymandering: Auditing and learning for subgroup fairness.  _arXiv:1711.05144_ , 2017.
  * Kipf & Welling (2016) Kipf, T. and Welling, M.  Variational graph auto-encoders.  In _NIPS Workshop on Bayesian Deep Learning_ , 2016.
  * Madras et al. (2018) Madras, D., Creager, E., Pitassi, T., and Zemel, R.  Learning adversarially fair and transferable representations.  2018\.
  * Majone (1972) Majone, G.  Social space and social distance: some remarks on metric methods in data analysis.  _Quality and Quantity_ , 6(2):239–252, 1972.
  * Mathieu et al. (2016) Mathieu, M. F., Zhao, J. J., Zhao, J., Ramesh, A., Sprechmann, P., and LeCun, Y.  Disentangling factors of variation in deep representation using adversarial training.  In _Advances in Neural Information Processing Systems_ , pp. 5040–5048, 2016.
  * Mnih & Teh (2012) Mnih, A. and Teh, Y. W.  A fast and simple algorithm for training neural probabilistic language models.  2012\.
  * Moon et al. (2017) Moon, C., Jones, P., and Samatova, N. F.  Learning entity type embeddings for knowledge graph completion.  In _CIKM_ , 2017.
  * Ng et al. (2001) Ng, A., Jordan, M., Weiss, Y., et al.  On spectral clustering: Analysis and an algorithm.  In _NIPS_ , 2001.
  * Perozzi et al. (2014) Perozzi, B., Al-Rfou, R., and Skiena, S.  Deepwalk: Online learning of social representations.  In _KDD_ , 2014.
  * Toutanova et al. (2015) Toutanova, K., Chen, D., Pantel, P., Poon, H., Choudhury, P., and Gamon, M.  Representing text for joint embedding of text and knowledge bases.  In _EMNLP_ , 2015.
  * Wang et al. (2016) Wang, D., Cui, P., and Zhu, W.  Structural deep network embedding.  In _KDD_ , 2016.
  * Xu et al. (2015) Xu, B., Wang, N., Chen, T., and Li, M.  Empirical evaluation of rectified activations in convolutional network.  _Deep Learning Workshop, ICML 2015_ , 2015.
  * Yao & Huang (2017) Yao, S. and Huang, B.  New fairness metrics for recommendation that embrace differences.  In _Workshop on Fairness, Accountability, and Transparency in Machine Learning (FATML)_ , June 2017.
  * Zemel et al. (2013) Zemel, R., Wu, Y., Swersky, K., Pitassi, T., and Dwork, C.  Learning fair representations.  In _International Conference on Machine Learning_ , pp. 325–333, 2013.

##  Appendix A Implementation Details

We implement each discriminator and adversarial filter as multi-layer
perceptrons (MLPs) with a leaky ReLU non-linearity between layers, and we use
the Adam optimizer with default parameters. Unless otherwise specified we use
$\lambda=1000$ for all experiments and datasets. For fair comparison our
discriminator during training time and subsequent sensitive attribute
classifier share the same architecture and capacity. Finally, for every step
performed by the main encoding model the Discriminator is updated $5$ times.
We found that this was necessary to provide a sufficient supervisory signal to
the main encoding model.

##  Appendix B FB15k-237 Details

To generate negative triplets we randomly sample either a head or tail entity
during training, with a ratio of $20$ negatives for each positive triplet. The
TransD model is trained for $100$ epochs with an embedding dimension of $20$,
selected using cross-validation, while the sensitive attribute classifers are
trained for $50$ epochs. The discriminators, sensitive attribute classifier
and adversarial filters are modelled as MLP’s with $4$,$4$ and $2$ layers
respectively. Lastly, we use the training, validation and testing splits
provided in the datasets.

##  Appendix C MovieLens1M

As with FB15k-237 we use model the discriminators and sensitive attribute
classifiers are modelled as MLP’s but $9$ layers with dropout with $p=0.3$
between layers while the adversarial filter remains unchanged from FB15k-237.
We found that regularization was crucial to the performance of main model and
we use BatchNorm after the embedding lookup in the main model which has an
embedding dimensionality of $30$. As only user nodes contain sensitive
attributes our discriminators do not compute losses using movie nodes.
Finally, to train our sensitive attribute classifier we construct a $90\%$
split of all users while the remaining user nodes are used for test. The same
ratio of train/test is used for the actual dataset which constains
users,movies and corresponding ratings for said movies. Finally, we train the
main model and sensitive attribute classifiers for $200$ epochs.

##  Appendix D Reddit

Like FB15k-237 we generate negative triplets by either sampling head or tail
entities which are either users or subreddits but unnlike FB15k-237 we keep
the ratio of negatives and positives the same. We also inherit the same
architectures for discriminator, sensitive attribute classifier and attribute
filters used in MovieLens1M. The main model however uses an embedding
dimensionality of $50$. Similar to MovieLens1M only user nodes contain
sensitive attributes and as such the discriminator and sensitive attribute
classifier does not compute losses with respect to subreddit nodes. Also, our
training set comprises of a $90\%$ split of all edges while the the remaining
$10\%$ is used as a test set. To test compositional generalizability we held
out $10\%$ of user nodes. Lastly, we train the main model for $50$ epochs and
the sensitive attribute classifier for $100$ epochs.

##  Appendix E Additional Results on Reddit

To the test degree of which invariance is affected by the number of sensitive
attributes we report additional results on the Reddit dataset. Specifically,
we report results for the Held out set with 20, 30, 40, and 50 sensitive
attributes. Overall, these results show no statistically significant
degradation in terms of invariance performance or task accuracy.

Table 4: Average AUC values across top-k sensitive attributes for Reddit. The
results are reported on a Held Out test of different combinations of
attributes.
Reddit |  Held out AUC
---|---
20 Sensitive Attributes | 0.569
30 Sensitive Attributes | 0.569
40 Sensitive Attributes | 0.556
50 Sensitive Attributes | 0.519