# A frequency and sensitivity tunable microresonator array for high-speed
quantum processor readout

J. D. Whittaker  D-Wave Systems, Inc., Burnaby, British Columbia, V5G 4M9,
Canada  L. J. Swenson  D-Wave Systems, Inc., Burnaby, British Columbia, V5G
4M9, Canada  M. H. Volkmann  D-Wave Systems, Inc., Burnaby, British Columbia,
V5G 4M9, Canada  P. Spear  D-Wave Systems, Inc., Burnaby, British Columbia,
V5G 4M9, Canada  F. Altomare  D-Wave Systems, Inc., Burnaby, British Columbia,
V5G 4M9, Canada  A. J. Berkley  D-Wave Systems, Inc., Burnaby, British
Columbia, V5G 4M9, Canada  B. Bumble  Jet Propulsion Laboratory, California
Institute of Technology, Pasadena, California 91109, USA  P. Bunyk  D-Wave
Systems, Inc., Burnaby, British Columbia, V5G 4M9, Canada  P. K. Day  Jet
Propulsion Laboratory, California Institute of Technology, Pasadena,
California 91109, USA  B. H. Eom  Jet Propulsion Laboratory, California
Institute of Technology, Pasadena, California 91109, USA  R. Harris  D-Wave
Systems, Inc., Burnaby, British Columbia, V5G 4M9, Canada  J. P. Hilton
D-Wave Systems, Inc., Burnaby, British Columbia, V5G 4M9, Canada  E. Hoskinson
D-Wave Systems, Inc., Burnaby, British Columbia, V5G 4M9, Canada  M. W.
Johnson  D-Wave Systems, Inc., Burnaby, British Columbia, V5G 4M9, Canada  A.
Kleinsasser  Jet Propulsion Laboratory, California Institute of Technology,
Pasadena, California 91109, USA  E. Ladizinsky  D-Wave Systems, Inc., Burnaby,
British Columbia, V5G 4M9, Canada  T. Lanting  D-Wave Systems, Inc., Burnaby,
British Columbia, V5G 4M9, Canada  T. Oh  D-Wave Systems, Inc., Burnaby,
British Columbia, V5G 4M9, Canada  I. Perminov  D-Wave Systems, Inc., Burnaby,
British Columbia, V5G 4M9, Canada  E. Tolkacheva  D-Wave Systems, Inc.,
Burnaby, British Columbia, V5G 4M9, Canada  J. Yao  D-Wave Systems, Inc.,
Burnaby, British Columbia, V5G 4M9, Canada<EMAIL_ADDRESS>
###### Abstract

Superconducting microresonators have been successfully utilized as detection
elements for a wide variety of applications. With multiplexing factors
exceeding 1,000 detectors per transmission line, they are the most scalable
low-temperature detector technology demonstrated to date. For high-throughput
applications, fewer detectors can be coupled to a single wire but utilize a
larger per-detector bandwidth. For all existing designs, fluctuations in
fabrication tolerances result in a non-uniform shift in resonance frequency
and sensitivity, which ultimately limits the efficiency of bandwidth
utilization. Here we present the design, implementation, and initial
characterization of a superconducting microresonator readout integrating two
tunable inductances per detector. We demonstrate that these tuning elements
provide independent control of both the detector frequency and sensitivity,
allowing us to maximize the transmission line bandwidth utilization. Finally
we discuss the integration of these detectors in a multilayer fabrication
stack for high-speed readout of the D-Wave quantum processor, highlighting the
use of control and routing circuitry composed of single-flux-quantum loops to
minimize the number of control wires at the lowest temperature stage.

##  I Introduction

Since their first demonstration in 2003,day:2003  superconducting
microresonators have become increasingly popular for low-temperature detector
applications. They are an attractive option due to their intrinsic frequency-
domain multiplexing capability, the wide availability of high-speed digital
electronics, and breakthroughs in material physics allowing very high
sensitivities to be achieved.leduc:2010 ; gao:2011 ; cecil:2012 ;
devisser:2011 ; barends:2008 ; vissers:2012 ; catalano:2015  Arrays of these
devices have been successfully demonstrated at sub-mm,monfardini:2010 ;
monfardini:2011 ; maloney:2010  far-infrared,swenson:2012 ; mckenney:2012
infrared/optical,mazin:2010 , and x-rayquaranta:2013  wavelengths, for phonon
imagingswenson:2010 ; cruciani:2012 ; moore:2012  and in quantum computing
circuity.jerger:2011 ; chen:2012 .

Despite these achievements, superconducting microresonator arrays typically
suffer from inefficient bandwidth utilization due to unavoidable fabrication
imperfections. As an example, the operating frequency $f$ of an LC resonator
utilizing a parallel plate capacitor suitable for integration with a multi-
layer fabrication stack will be sensitive to perturbations of the dielectric
thickness $d$ according to the expression $\delta f/f=\delta d/2d$. A
realistic $\pm 10\%$ tolerance in layer thickness for a modern fabrication
facility results in up to $\pm$300 MHz shift for a nominal 6 GHz
microresonator. Global variations affecting all devices uniformly can usually
be compensated by a simple shift in the readout center frequency, however,
local variations are more problematic. These cause each microresonator to
exhibit a resonance frequency $f_{0}$ which deviates randomly from its
expected value - often by more than a linewidth. The array designer is thus
left with a choice of either sparsely populating the available electronics
bandwidth with resonances such that stochastic resonance collisions are
avoided, or utilizing a high multiplexing factor and accepting a decreased
array yield due to overlapping resonances.

Figure 1: Schematic for a single FASTR detector. As described in the text, the
microresonator parameters are adjusted with local magnetic flux biases to the
DC-SQUID loops which act as tunable nonlinear inductors. Tuning these biases
allows in-situ compensation for the unavoidable device variations inherent in
fabrication. While both DC-SQUIDs are tunable, only a single loop, labeled
SENSE, is coupled to the measurement system. For the D-Wave quantum processor,
this is the last stage of a data shift register made from quantum flux
parametrons (QFPs).

One way to make an array with efficient bandwidth utilization is to use
superconducting microresonators whose frequency and sensitivity can be tuned
in-situ, allowing compensation for local fabrication variations. Tuning of
these parameters allows an array of detectors to exhibit homogeneous
sensitivity and uniform frequency spacing. While superconducting
microresonators with tunable frequencies have been previously
demonstrated,mates:2008 ; whittaker:2014 ; vissers:2015 ; palacios-laloy:2008
this single adjustment knob is insufficient in the context of an efficiently
packed array. The detector sensitivity is equally important when the array is
subject to a signal with a large dynamic range, as more sensitive detectors
will exhibit a larger frequency shift than less sensitive detectors resulting
in detector frequency collisions. Adding a second tuning knob that also
provides sensitivity tuning allows maximum packing efficiency to be achieved
for a fixed electronic bandwidth. A schematic of a Frequency And Sensitivity
Tunable Resonator (FASTR) detector is shown in Fig. 1.

The tuning parameters are realized using two DC-SQUID loops in series with a
fixed geometric inductance $L_{g}$. Each of these DC-SQUIDs are nonlinear
inductances that can be adjusted using an external flux bias. This inductive
branch is in parallel with a shunt capacitor $C_{s}$. External flux can be
applied locally to each DC-SQUID loop either with individual analog wires,
requiring $2N$ wires for an array of $N$ detectors, or with an array of flux
DACs.bunyk:2014  This latter option can take advantage of an addressing
network which requires only $O(\sqrt[3]{2n_{res}})$ wires, making control of
even very large arrays straightforward. Flux sensitivity is achieved by
coupling one of the two loops to the system of interest. Proper selection of
the DC-SQUID biases is required to achieve a desired operating frequency and
sensitivity. The measured detector frequency of a prototype device (discussed
below) for various flux biases in the two DC-SQUID loops is shown in Fig. 2a.

Figure 2: a) The measured resonance frequency of a prototype device to applied
bias in the two DC-SQUID loops. The unbiased resonance frequency of this
device is $f_{0}=$ 6.91 GHz. The black line is a contour of constant frequency
at $f_{0}=6.83$ GHz and the linewidth is 21 MHz. b) Response $\left|\delta
f_{0}\right|$ in linewidths while moving along the contour. The signal in this
case was the flux change induced by the two states of the QFP magnetically
coupled to the SENSE SQUID. Note that for the low resonator quality factor
$Q_{r}$ and large detector bandwidth discussed here, the noise for this
readout is fixed and dominated by the readout electronics. In this case,
changes in responsivity directly impact the detector sensitivity.

A contour of constant frequency has been drawn on the surface at a frequency
shift of about $-$3.5 linewidths from the zero-flux resonance frequency. Only
the SENSE SQUID is coupled to the measurement system. For the schematic shown
in Fig. 1, the FASTR detector is coupled to a quantum flux parametron (QFP)
which is the last stage of a data shift register. Modulating the state of the
QFP while following this contour allows the device responsivity to be
measured. The result of this measurement is shown in Fig. 2b. Bias selection
then consists of first choosing a contour of constant frequency followed by a
traversal of the contour until the desired responsivity has been achieved.

Here we discuss the real-world example of utilizing a frequency-multiplexed
array of FASTR detectors for readout of the D-Wave Two quantum processor. We
first briefly discuss the readout constraints imposed by the D-Wave Two
quantum processor. We then present a parametric design and physical
realization of a single prototype FASTR detector suitable for processor
integration. While this specific application demonstrates a proof-of-principle
of a frequency- and sensitivity-tunable superconducting microresonator, it
should be kept in mind that this device can be broadly applied to a wide
variety of magnetic-flux sensing applications including readout of the family
of single flux quantum (SFQ) digital circuity.

##  II Readout constraints of the D-Wave Two quantum processor

The current D-Wave quantum processor utilizes the quantum annealing algorithm
to solve for the low lying energy states of the Ising spin Hamiltonian

$H=\sum\limits_{i}h_{i}\sigma_{i}^{z}+\sum\limits_{i<j}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$
(1)

where the $\sigma_{i}$ are Pauli matrices for the $i$th spin, $h_{i}$ are
local energy scaling factors, and the non-zero coupling terms $J_{ij}$ reflect
the underlying connectivity of the processor graph. For the D-Wave Two
processor, the graph consists of 64 repeated 8 qubit unit cells. Within a cell
the qubits are connected in a K4,4 planar topology. Each qubit is further
connected to two qubits in neighboring unit cells to yield 6 connections per
processor qubit. The processor topology, qubit parametric design, problem
specification (consisting of setting the $h_{i}$ and $J_{ij}$ terms in Eq. 1),
and the annealing algorithm have been previously described.harris:2010 ;
harris:2010b ; johnson:2011 ; bunyk:2014  It is important to note that during
the annealing procedure qubit entanglementlanting-ent:2014  is eventually
frozen out and all qubits are projected into a classical flux state before
readout.

Figure 3: Picture of the D-Wave Two processor chip, made up of 64 eight-qubit
unit cells. The shift register streets are highlighted yellow and lead to the
outside of the processor where it is straightforward to locate FASTR
detectors. Figure 4: Simplified schematic of a portion of the QFP-based shift
register. $Q_{1}$ and $Q_{2}$ are qubits and $\phi_{0...3}$ are flux biases
that control data flow.

Processor readout entails detecting the flux state of all qubits at the end of
a computation sequence. The readout architecture consists of a shift register
lattice interlaced between the 8-qubit unit cells. A picture of the processor
chip, with the shift register lattice highlighted in yellow, is shown in Fig.
3. The shift register eliminates the need for an individual detector per qubit
and is used to move state information from individual qubits out to detectors
on the chip perimeter. A simplified portion of the shift register is shown
schematically in Fig. 4.

A shift register stage is composed of a QFP that can be modulated between a
monostable unlatched and a bistable latched state using a local flux bias
$\phi_{n}$. In the latched state, information is stored in the QFP as either
circulating or counter-circulating current. This current is a flux source
which applies a flux bias to the neighboring shift register stages held in the
unlatched state. The next stage in the forward direction is subsequently
latched before the source stage is unlatched allowing information to be passed
down the shift register. Note that two unlatched stages are required between
every latched stage in order to protect against back-action from data further
ahead in the shift register. As indicated in Fig. 4, three QFP flux biases
$\phi_{1...3}$ are thus required to pass information along the shift register.
A separate QFP stage, biased by $\phi_{0}$, is directly connected to each
qubit and is used to selectively copy data from the qubit to the rest of the
shift register.

The topology and operating speed of the shift register set crucial constraints
for the FASTR readout array. For the processor shown in Fig. 3, there are 8
vertical and 8 horizontal shift registers. As data can be moved in either
direction along the shift register, there are 32 points where detectors can
easily be placed. Note, however, that the shift register can only move in one
direction at a time so that only half of these detectors can be operated
simultaneously. The detector on the other end of each line provides redundancy
for shift register breaks caused by inoperable QFPs. While in principle it is
possible to increase the shift register fan-out by utilizing more stages, this
option increases circuit complexity.

Figure 5: Simplified system electronics schematic for FASTR array readout. The
digital-to-analog and analog-to-digital converters are phase locked and
sampled at 2.5 GSPS.

The shift register operation speed is set by the bandwidth of the QFP flux
bias lines $\phi_{n}$. Currently these lines utilize a 30 MHz low-pass filter
to limit noise. As the shift register is composed of three biases
$\phi_{1...3}$, each shift register line can deliver data at $\sim$10 Mbits/s.
If only a single detector is attached to each shift register line, the
detector bandwidth should match or somewhat exceed 10 MHz.

Although the attainable readout speed for a FASTR detector is set by the
microresonator linewidth, more than a single linewidth of electronic bandwidth
is required per detector. Additional bandwidth is required to accommodate
shifting of the resonance frequency due to modulation of the signal coupled
into the SENSE SQUID. Even more electronic bandwidth is required to avoid
electronic crosstalk with neighboring microresonators, as the resonance
appreciably perturbs the microwave transmission beyond the nominal linewidth.

Given the preceding constraints, the FASTR readout for the D-Wave Two quantum
processor is designed to utilize resonances with 19.5 MHz linewidths,
modulated by 1 linewidth by the shift register data. Each detector is allotted
4 linewidths of electronic bandwidth to minimize crosstalk. A total of 2.5 GHz
of readout electronic bandwidth is therefore required to readout all 32
detectors on a single transmission line. Considering this, the current D-Wave
FASTR readout is designed to operate using a 2.5 GHz bandwidth with an array
center frequency of 6 GHz. A simplified schematic of the readout circuit is
shown in Fig. 5.

The readout design is similar to Kinetic Inductance Detector (KID) based
readouts with the notable feature of a large 2.5 GSPS sampling frequency of
the baseband digitizing electronics.bourrion:2011 ; yates:2009 ; mchugh:2012
In contrast to KID readouts used for imaging which continuously interrogate
the microresonator array, the FASTR readout only measures the array state
after completion of a shift register operation cycle. The $\sim$35$\%$ duty
cycle reduces the data rate to a level at which the Field Programmable Gate
Array (FPGA) used for signal processing can continuously sustain. As
previously mentioned, the center frequency of the 2.5 GHz of readout bandwidth
should be adjustable in order to compensate for any global shift in the
detector frequencies. For the current system, global shifts up to $\pm$750 MHz
can be tolerated by adjusting the local oscillator frequency. This limit is
imposed by the 4-8 GHz operating band of the commercial microwave components
used in the readout chain.

Other constraints are imposed when integrating FASTR devices (Fig. 1) into the
D-Wave Two quantum processor, which utilizes a multilayer fabrication stack
consisting of 6 niobium metal layers separated by planarized dielectric
layers. A particular challenge is to realize FASTR devices with this
multilayer stack that achieve sufficiently high intrinsic quality factors
$Q_{i}$. The linewidth requirement discussed previously dictates that the
microresonator quality factors $Q_{r}$ are of order
$f_{0}/\Delta{f}=6\mathrm{e}9/19.5\mathrm{e}6\sim 300$. In order achieve full
signal modulation and minimize processor heating due to dissipation by the
FASTR readout, it is desirable for $Q_{i}$ to be as high as possible.
Achieving the limit $Q_{i}\gg Q_{r}$ has the added advantage that $Q_{r}$ is
then simply determined by the coupling to the microwave circuit and can be
easily set in design with the coupling capacitance $C_{c}$.

The primary source of dissipation in superconducting microresonators has been
extensively studied,martinis:2005 ; oconnell:2008 ; gao:2008b  and low-loss
microresonators made of interdigitated capacitors on crystalline substrates
routinely achieve intrinsic quality factors exceeding a million.leduc:2010 ;
geerlings:2012 ; chang:2013  This low loss is achieved by engineering
resonators with the electric field energy constrained to vacuum or the
underlying low-loss crystalline substrate. The surface, where lossy dangling
bonds and oxides are present, are specifically avoided in these
devices.gao:2008c ; quintana:2014 ; dial:2015  For example, one technique that
is commonly used to improve the $Q_{i}$ of planar interdigitated capacitor is
to increase the finger spacing which decreases the participation ratio of the
surface defects.

Compared with single-layer microresonator designs on a crystalline substrate,
microresonators incorporating deposited dielectrics suitable for integration
in a multilayer fabrication stack typically exhibit much higher loss due to
the presence of two-level systems (TLS) in the amorphous
dielectric.martinis:2014  In this case, simple geometric changes such as
increasing the finger spacing of an interdigitated capacitor do not result in
an improved device $Q_{i}$ as the electric field is displaced but remains
within the lossy deposited dielectric. Fortunately, these dissipation pathways
become saturated above a critical electric field, such that $Q_{i}$ increases
with drive power. In this case, by minimizing the dielectric volume used for
fabricating the resonator dielectric, saturation can be achieved for a minimal
stored resonator energy. This allows high intrinsic quality factors and
minimal heating to be achieved for modest microwave drive powers.

Considering the need to minimize the volume of lossy dielectric, a parallel
plate geometry with a minimal dielectric thickness is particularly suitable
for the integrated FASTR design. Finally, a lumped-element geometry utilizing
a parallel-plate capacitor comes with the added advantage that, as opposed to
distributed resonator geometries, no harmonics are inherent in the design.
This allows the readout to occupy more than a single octave of bandwidth which
is an important consideration for scaling to large array sizes.

##  III Characterization of a prototype device

Figure 6: Prototype FASTR detector with three QFP shift register stages,
designed to have $C_{s}=1.7$ pF, $C_{c}=70$ fF, $Q_{c}=338$, $L_{g}=324$ pH,
$f_{0}=6.628$ GHz, and Josephson junction $I_{c}=11$ $\mu$A. The fabrication
stack is also shown, highlighting the 50 nm thick capacitor dielectric layer.

We have designed, fabricated, and measured a prototype FASTR detector suitable
for integration with the D-Wave Two quantum processor. The layout and
fabrication stack are shown in Fig. 6. Testing was performed at 10 mK. The
detector resonance frequency with zero flux in either SQUID loop was measured
to be 6.91 GHz. This is a large but still acceptable 300 MHz above the design
value and attributable to fabrication variation in the capacitor dielectric
layer thickness. While only a single prototype FASTR detector was
characterized in this initial measurement, the multiplexed readout system
shown in Fig. 5 was used for electrical readout while operating at the full
2.5 GSPS data rate.

Figure 7: Resonator intrinsic quality factor $Q_{i}$ and nonlinear shift $a$
as functions of the drive power, under typical bias conditions at $f_{0}=6.84$
GHz. These parameters have been extracted from the data shown in the inset.
The resonator should be operated at a sufficient power to provide a large
amplitude SNR required for high fidelity readout while avoiding intrinsic
device nonlinearity. For the prototype FASTR detector, the region with an SNR
$>$ 5 and weak device nonlinearity ($a$ $<$ 0.25) is indicated as the operable
region. The coupling quality factor $Q_{c}$ for this device was 329 and has an
unbiased resonance frequency of $f_{0}=6.910$ GHz.

The fabricated prototype device featured a parallel plate capacitor with a
nominally 50 nm thick sputtered silicon oxide dielectric. $S_{21}$
transmission data were taken and fit to the expression for a transmission line
shunted by a parallel RLC resonator

$S_{21}=1-\frac{Q_{r}}{Q_{c}}\frac{1}{1+2iQ_{r}x}$ (2)

where $x=\left(f-f_{0}\right)/f_{0}$ is the fractional detuning from
resonance. The results are shown in Fig. 7 for a variety of powers. The low-
power intrinsic quality factor of the resonator with this capacitor was
measured to be $Q_{i}\sim 1000$. The fit results show $Q_{i}$ changing with
power (red dots), exhibiting two-level system (TLS) saturation behavior with a
saturation field of $|\vec{E}|=50$ V/m. While the low-power $Q_{i}$ is not
much greater than the $Q_{r}\sim 300$ required for the 32-resonator design,
above the saturation threshold $Q_{i}$ increases dramatically.

Fig. 7 also shows how the resonance frequency changes with drive power (blue
triangles) in units of linewidths according to
$a=\left(f_{0}^{\prime}-f_{0}\right)Q_{r}/f_{0}^{\prime}$, where
$f_{0}^{\prime}$ is the low-power resonance frequency. This non-linear
response to stored energy is characterized by the Duffing oscillator
nonlinearity parameter
$a=\left(\alpha{Q_{r}}/4\right)\left(I/I_{c}\right)^{2}$.zmuidzinas:2012 ;
swenson:2013  Here, $\alpha=L_{J}/(L_{g}+L_{J})$ characterizes the
contribution of the Josephson Junction kinetic inductance $L_{J}=\Phi_{0}/2\pi
I_{c}$ to the total inductance, $I$ is the current through the resonator, and
bifurcation occurs at $a\sim 0.77$ linewidths. Keeping the microwave drive
power sufficiently low assures only minimal nonlinearity ($a<0.25$ linewidths)
and well-behaved detector operation. This requirement sets the upper bound to
the operable region labeled in Fig. 7. The lower bound is set by the a minimum
desired amplitude SNR of 5 which gives a bit error probability less than
$10^{-6}$. We have measured the system noise temperature to be 7.9 K, limited
primarily by the first stage amplifier and the insertion loss of the cryogenic
isolators used to prevent amplifier back-action from disturbing the processor.
The microwave drive power should exceed $P_{g}=$ -96 dBm to achieve the
desired SNR assuming the tone is fully modulated by the 1 linewidth resonance
shift. Full modulation is achieved by requiring $Q_{i}\gg Q_{c}$ and biasing
symmetrically around the resonance, which maps the two qubit states onto the
resonance curve as shown in Figure 8a.

Figure 8: a) The two QFP states, 0 and 1, shown at $\pm$1/2 linewidth from the
FASTR resonance $f_{0}$, which is achieved when both biases are properly
tuned. b) Complex plane plot of prepared qubit states read out by the FASTR
system, at a SNR$\sim$5\. The data have been translated and rotated for simple
state discrimination along $Q=0$.

##  IV FASTR readout of a qubit state and noise

As a proof-of-principle of full-system operation, the prototype FASTR detector
was used to readout a qubit via the three-stage shift register shown in Fig.
6. The detector was driven at $P_{g}=-98$ dBm ($a=0.05$ linewidths), slightly
below the SNR = 5 lower limit of the operable region indicated in Fig. 7. A
large positive flux bias applied to the qubit, followed by annealing, allows
for preparation in a known state. Using a large negative flux bias allows the
qubit to be prepared in the opposite state. After qubit annealing, the
resulting state can be passed along the shift register and read out with the
FASTR detector. Repeating this procedure many times with a known data pattern
allows the readout fidelity to be assessed. A plot of the measured fidelity
data is shown in Fig. 8b, where a linear transformation to the complex
transmission has been applied to simplify state discrimination. This process
of transformation and discrimination was applied to all subsequent qubit
measurements. From the fidelity data, a bit error probability of $\sim
10^{-5}$ was calculated, which matches well with expectation for the device
parameters and the microwave drive power utilized. Note that while an SNR $<$
5 is sufficient for verification of this prototype device, the measurement
power will be properly adjusted before full processor readout to ensure the
desired bit error probability of $<10^{-6}$ is achieved.

Figure 9: a) Qubit population measurement, read out with a FASTR, as qubit
flux bias was swept. b) Noise PSD measurement and fit at a sampling rate
$1/\tau_{s}=1/3.6\mu{s}$. The $1/f$ noise amplitude at 1 Hz was 11
$\mu\Phi_{0}/\sqrt{Hz}$, while the white noise was 0.8
$\mu\Phi_{0}/\sqrt{Hz}$.

The previous measurements validate the operation of a FASTR detector and the
readout system. It is now possible to repeat two of D-Wave’s standard qubit
characterization measurements: a qubit population measurement and a noise
power spectral measurement (PSD).lanting:2009 ; lanting:2014  The qubit
population measurement was made while sweeping the qubit flux bias through
degeneracy, with the results shown in Fig. 9a. The fit was done with the
equation

$P=\frac{1}{2}\left[1+\tanh\left(\frac{\Phi_{x}}{2W}\right)\right],$ (3)

where $\Phi_{x}$ is the externally applied flux and $2W$ is the width of the
transition. Inverting this fit gives a mapping between qubit population and
$\Phi_{x}$. When annealed at a fixed bias point near the center of the
transition, on average the qubit population should always return the same
value. However, noise will cause the population to wander. Repeated
measurement of the population and conversion into an equivalent external flux
allows a noise PSD to be calculated. The result of this measurement using the
prototype FASTR detector is shown in Fig. 9b. The noise spectra show a
combination of $1/f$ noise and white noise $w_{n}$, where
$w_{n}=\tau_{s}4W^{2}$. The $1/f$ noise component is due to intrinsic device
noise, and the white noise is due to the statistical nature of the qubit
measurement and drops with increased sample rate or a narrower transition
width. Note that traditional DC-SQUIDs are dissipative. Fast operation of
these devices results in an increased sample temperature and consequently a
wider transition width $W$. This heating can be reduced by increasing the time
per sample $\tau_{s}$. However in either case $w_{n}$ is
increased.lanting:2009 ; lanting:2014  Thus the low dissipation and fast duty
cycle allowed by a FASTR detector results in a low white noise level as
compared with a DC SQUID readout.

##  V Perspective on readout scaling

Scaling of the FASTR readout is constrained by the achievable device $Q_{i}$,
available readout bandwidth, and the number of analog wires required to tune
the individual resonance frequencies. Based on commercially available
electronics, the FASTR array was designed for a 2.5 GHz band centered around 6
GHz. Placing resonators at the ends of each shift register street, the number
of resonators for the D-Wave processor goes as $n_{res}=4\sqrt{n_{cells}}$.
Table 1 shows the scaling of the per-resonator bandwidth and required device
$Q_{i}$ for a fixed 2.5 GHz bandwidth readout design.

$N_{qubits}$ | $n_{cells}$ | $n_{res}$ |  $\Delta f$ (MHz) | $Q_{c}$ | $Q_{i}$ | $N$
---|---|---|---|---|---|---
512 | 64 | 32 | 19.5 |  $\sim$240-370 |  $>$3,700 | 4
1152 | 144 | 48 | 13.0 |  $\sim$360-560 |  $>$5,600 | 5
2048 | 256 | 64 | 9.8 |  $\sim$490-740 |  $>$7,400 | 6
3200 | 400 | 80 | 7.8 |  $\sim$610-930 |  $>$9,300 | 6
4608 | 576 | 96 | 6.5 |  $\sim$730-1110 |  $>$11,100 | 6
Table 1: Scaling of the FASTR readout scheme with the size of processor. $N$
is the number of additional analog wires required as discussed in the text.
The stated $Q_{i}$ assumes the need to meet the condition $Q_{i}\gg Q_{c}$ to
minimize chip heating and design effort.

Figure 7 shows that for the SiOx dielectric used for the prototype device, a
$Q_{i}$ of $\sim$6,000 was achieved in the operating regime. Increasing the
number of devices per line beyond $\sim$48 thus requires incorporation of
lower loss deposited dielectrics such as hydrogen-rich amorphous
silicon.mazin:2010b

A straightforward way to scale the readout is to increase the readout
bandwidth. Already, digitizing components with bandwidths exceeding 5 GHz and
sufficient dynamic range are commercially available. As the design here
utilizes a parallel-plate capacitor incorporated in a lumped-element resonator
geometry which does not exhibit harmonics, more than one octave of bandwidth
can be readily utilized. Perhaps the greatest difficulty for implementing
large bandwidth electronics is the very high data rates that must be handled
by the signal-processing hardware. Fortunately, very large FPGAs are becoming
available that should be sufficient for this task.

Each FASTR device requires two external flux biases. $2n_{res}$ lines will be
needed if each of these are supplied with individual analog lines. A
32-resonator system would require already a somewhat impractical 64 analog
lines. As a better alternative, the D-Wave Two processor uses superconducting
digital-to-analog flux converters to apply flux biases to thousands of devices
though an addressing scheme that scales as $O(\sqrt[3]{2n_{res}})$.bunyk:2014
For instance, a 96-resonator system would require the addition of only 6
analog lines using this addressing scheme, which is far more reasonable than
the 192 lines a per-device addressing scheme would require.

##  VI Conclusion

We have demonstrated a frequency- and sensitivity-tunable resonator detector
suitable for a wide variety of magnetic flux sensing applications. The
detector tuning parameters allow for compensation of detector-scale
fabrication variations, while shifting the center frequency of the readout
electronics allows for compensation of wafer-scale variations. Using a
combination of these two techniques, the resonator frequencies can be spaced
uniformly and a homogenized flux-sensitivity can be realized. In this way a
maximum packing efficiency can be achieved for a fixed electronic bandwidth.
The utility of this device includes readout of the D-Wave quantum processor,
SFQ digital circuitry, and indeed any magnetic flux sensing applications
requiring arrays of densely multiplexed detectors.

##  VII Acknowledgements

The authors would like to thank R. Neufeld and D. Walliman for providing
photographs of samples, as well as C. Enderud, C. Baron, and M. Babcock for
sample preparation and cryogenic support. A portion of this research was
carried out at the Jet Propulsion Laboratory, California Institute of
Technology, under a contract with the National Aeronautics and Space
Administration.

## References

  * (1) P. K. Day, H. G. LeDuc, B. A. Mazin, A. Vayonakis, and J. Zmuidzinas,  Nature 425, 817821 (2003).
  * (2) H. G. Leduc et al.,  Appl. Phys. Lett. 97, 102509 (2010).
  * (3) J. Gao et al.,  Appl. Phys. Lett. 98, 232508 (2011).
  * (4) T. Cecil et al.,  Appl. Phys. Lett. 101, 032601 (2012).
  * (5) P. J. de Visser et al.,  Phys. Rev. Lett. 106, 167004 (2011).
  * (6) R. Barends et al.,  Appl. Phys. Lett. 92, 223502 (2008).
  * (7) M. R. Vissers, M. P. Weides, J. S. Kline, M. Sandberg, and D. P. Pappas,  Applied Physics Letters 101, 022601 (2012).
  * (8) Catalano, A. et al.,  A&A 580, A15 (2015).
  * (9) A. Monfardini et al.,  A&A 521, A29 (2010).
  * (10) A. Monfardini et al.,  ApJ Supplement Series 194, 24 (2011).
  * (11) P. R. Maloney et al.,  Proc. SPIE Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy V 7741, 77410F (2010).
  * (12) L. J. Swenson et al.,  Proc. SPIE 8452, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VI , 84520P (2012).
  * (13) C. M. McKenney et al.,  Proc. SPIE 8452, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VI , 84520S (2012).
  * (14) B. A. Mazin et al.,  Proc. SPIE 7735, Ground-based and Airborne Instrumentation for Astronomy III , 773518 (2010).
  * (15) O. Quaranta, T. W. Cecil, L. Gades, B. Mazin, and A. Miceli,  Superconductor Sci. and Tech. 26, 105021 (2013).
  * (16) L. J. Swenson et al.,  Appl. Phys. Lett. 96, 263511 (2010).
  * (17) A. Cruciani et al.,  Journal of Low Temperature Physics 167, 311 (2012).
  * (18) D. C. Moore et al.,  Applied Physics Letters 100, 232601 (2012).
  * (19) M. Jerger et al.,  Europhys. Lett. 96, 40012 (2011).
  * (20) Y. Chen et al.,  Appl. Phys. Lett. 101, 182601 (2012).
  * (21) J. A. B. Mates, G. C. Hilton, K. D. Irwin, L. R. Vale, and K. W. Lehnert,  Appl. Phys. Lett. 92, 023514 (2008).
  * (22) J. D. Whittaker et al.,  Phys. Rev. B 90, 024513 (2014).
  * (23) M. R. Vissers et al.,  Applied Physics Letters 107, 062601 (2015).
  * (24) A. Palacios-Laloy et al.,  J. Low Temp. Phys. 151, 1034 (2008).
  * (25) P. Bunyk et al.,  IEEE Trans. Appl. Superconductivity 24, 1 (2014).
  * (26) R. Harris et al.,  Phys. Rev. B 81, 134510 (2010).
  * (27) R. Harris et al.,  Phys. Rev. B 82, 024511 (2010).
  * (28) M. W. Johnson et al.,  Nature 473, 194 (2011).
  * (29) T. Lanting et al.,  Phys. Rev. X 4, 021041 (2014).
  * (30) O. Bourrion et al.,  J. Inst. 6, P06012 (2011).
  * (31) S. J. C. Yates, A. M. Baryshev, J. J. A. Baselmans, B. Klein, and R. GÃŒsten,  Appl. Phys. Lett. 95, 042504 (2009).
  * (32) S. McHugh et al.,  Rev. Sci. Inst. 83, 044702 (2012).
  * (33) J. M. Martinis et al.,  Phys. Rev. Lett. 95, 210503 (2005).
  * (34) A. D. O’Connell et al.,  Appl. Phys. Lett. 92, 112903 (2008).
  * (35) J. Gao et al.,  Appl. Phys. Lett. 92, 212504 (2008).
  * (36) R. Barends et al.,  Appl. Phys. Lett. 92, 152505 (2008).
  * (37) C. Quintana et al.,  Appl. Phys. Lett. 105, 062601 (2014).
  * (38) O. Dial et al.,  arXiv:1509.03859 (2015).
  * (39) K. Geerlings et al.,  Appl. Phys. Lett. 100, 192601 (2012).
  * (40) J. B. Chang et al.,  Appl. Phys. Lett. 103, 012602 (2013).
  * (41) J. M. Martinis and A. Megrant,  arXiv:1410.5793 (2014).
  * (42) J. Zmuidzinas,  Annu. Rev. Cond. Mat. Phys. 3, 169 (2012).
  * (43) L. J. Swenson et al.,  J. Appl. Phys. 113, 104501 (2013).
  * (44) T. Lanting et al.,  Phys. Rev. B 79, 060509 (2009).
  * (45) T. Lanting et al.,  Phys. Rev. B 89, 014503 (2014).
  * (46) B. A. Mazin et al,  Appl. Phys. Lett. 96, 102504 (2010).


11institutetext: French National Cybersecurity Agency (ANSSI), Paris, France
11email<EMAIL_ADDRESS>
# A Source Separation Approach to Temporal Graph Modelling for Computer
Networks

Corentin Larroche

###### Abstract

Detecting malicious activity within an enterprise computer network can be
framed as a temporal link prediction task: given a sequence of graphs
representing communications between hosts over time, the goal is to predict
which edges should—or should not—occur in the future. However, standard
temporal link prediction algorithms are ill-suited for computer network
monitoring as they do not take account of the peculiar short-term dynamics of
computer network activity, which exhibits sharp seasonal variations. In order
to build a better model, we propose a source separation-inspired description
of computer network activity: at each time step, the observed graph is a
mixture of subgraphs representing various sources of activity, and short-term
dynamics result from changes in the mixing coefficients. Both qualitative and
quantitative experiments demonstrate the validity of our approach.

###### Keywords:

Cybersecurity Anomaly detection Temporal graphs.

##  1 Introduction

Besides preventive mechanisms such as firewalls and system hardening, securing
computer networks requires the ability to detect ongoing intrusions in order
to stop them in their tracks [25]. This need for effective intrusion detection
systems has motivated the emergence of several research topics, including
graph-based computer network monitoring, which represents communications
between the hosts of the network as a temporal graph [19, 23, 27]. Intrusion
detection then boils down to detecting anomalous edges in this graph. The
dominant method in the literature is to reframe this problem as a temporal
link prediction task, then solve this task using existing models originally
designed for social network analysis [26] or content recommendation [32].
However, the short-term temporal dynamics of computer network activity differ
from those of social networks and user–content interaction streams. In
particular, communication patterns within an enterprise network change
significantly over the course of a day and exhibit strong seasonality [31].
Therefore, models designed for social networks or recommender systems are not
well suited for computer network monitoring.

To take better account of these peculiarities, we formulate the following
hypothesis: at each time step, the observed activity is in fact the
superposition of several activity sources (e.g., office activity, network
administration, automated traffic), and the main cause for short-term dynamics
is a change in the respective importance of each source. This setting is
illustrated in Figure 1, and it provides a simple explanation for the steep,
global changes observed in communication graphs over time. Following this
hypothesis, modelling computer network activity can be seen as a source
separation problem: given a sequence of observed graphs, we aim to learn a
model for each activity source as well as the mixing coefficients for each
graph. Forecasting future activity then only requires predicting a small
number of mixing coefficients. In other words, such a model would have less
degrees of freedom than standard temporal link predictors. The expected upside
of this restrained expressivity is an increase in scalability,
interpretability and robustness to noise, in turn leading to more reliable
anomaly detection.

Figure 1: The intuition behind our approach: activity within a computer
network can be seen as a superposition of several activity sources, and the
time-varying influence of each source is the main cause for variation in the
global communication graph.

#### Contributions.

We propose a simple model for temporal link prediction in computer networks,
called Superposed Nonnegative Matrix Factorization (SNMF), of which we provide
an open-source implementation 111https://github.com/cl-
anssi/NetworkSourceSeparation. Qualitative and quantitative experiments
demonstrate the relevance of both SNMF itself and the underlying hypothesis on
the nature of the short-term dynamics of computer network activity. In
particular, SNMF achieves state-of-the-art performance on a public real-world
computer network monitoring dataset [17].

The rest of this paper is structured as follows. We formally define computer
network monitoring and its relation to temporal link prediction, then discuss
related work in Section 2. Section 3 introduces SNMF and the associated
inference procedure. Our experiments are described in Section 4, and we
finally discuss limitations and potential extensions of our work in Section 5.

##  2 Problem Statement and Related Work

Before describing our approach to computer network monitoring, we first
provide a formal definition of this task in Section 2.1 and review previous
contributions in Section 2.2.

###  2.1 Problem Statement — Computer Network Monitoring

#### Key notations.

Throughout this paper, we describe communications within a computer network as
a sequence of directed graphs $(\mathcal{G}_{t})_{t\geq{1}}$. These graphs all
share the same $N$ nodes, each of which represents one host. Each graph
$\mathcal{G}_{t}$ corresponds to a fixed-length time window, and there is an
edge $(i,j)$ in $\mathcal{G}_{t}$ if host $i$ has contacted host $j$ within
the $t$-th window. These edges are encoded as a binary adjacency matrix
$\mathbf{A}_{t}\in\\{0,1\\}^{N\times{N}}$. The set of all possible temporal
edges is denoted
$\mathcal{E}=\\{(i,j,t)\in[N]\times[N]\times\mathbb{N}^{*}:i\neq{j}\\}$, where
$[N]$ is short for $\\{1,\ldots,N\\}$. Note that self-loops are excluded.
Finally, the element-wise product (resp. quotient) between two matrices or
vectors is denoted $\odot$ (resp. $\div$), the Frobenius inner product between
two matrices is denoted $\langle\cdot\,,\cdot\rangle_{\mathrm{F}}$, and the
$N\times{N}$ all-one and identity matrices are called $\mathbf{1}_{N}$ and
$\mathbf{I}_{N}$, respectively.

#### Problem statement.

We define computer network monitoring as an anomaly ranking problem. Given a
training set consisting of $T\geq{1}$ graphs
$(\mathcal{G}_{1},\ldots,\mathcal{G}_{T})$ representing $T$ successive time
windows containing normal activity, our goal is to learn a link predictor
$h:\mathcal{E}\rightarrow\mathbb{R}$ such that $h(i,j,t)$ is higher for normal
edges than for anomalous edges. More specifically, the predictor $h$ should
minimize the false positive rate (i.e., the proportion of legitimate
connections erroneously flagged as anomalous) under some detection rate
constraint, meaning that edges resulting from malicious behavior should be
detected with high probability.

Since traces of truly malicious behavior are hard to obtain when training and
validating the model, the actual detection rate cannot be easily estimated. As
a consequence, this anomaly ranking problem is usually reframed as a temporal
link prediction task. In other words, the predictor $h$ should assign higher
scores to normal edges than to edges that do not occur as part of normal
activity. This connection with link prediction underpins most existing work on
computer network monitoring, which we further discuss in the next section.

###  2.2 Related Work — Temporal Graph Models

Early contributions on computer network monitoring [27, 28, 38] adopt a local
modelling approach: each host pair $(i,j)$ is treated independently, and its
usual activity patterns are inferred based solely upon its own past activity.
The main shortcoming of this approach is that it performs poorly on new edges,
which have no past activity by definition. As a consequence, interest quickly
shifted towards more sophisticated models designed for link prediction.

#### Latent space models.

Despite their differences, the most popular link prediction algorithms all
rely upon the following idea: each node $i$ is associated with two dense
embeddings $\mathbf{u}_{i},\mathbf{v}_{j}\in\mathbb{R}^{K}$, and the presence
of a link from $i$ to $j$ is predicted as
$h(i,j)=g(\mathbf{u}_{i},\mathbf{v}_{j})$ for some function
$g:\mathbb{R}^{K}\times\mathbb{R}^{K}\rightarrow\mathbb{R}$. This framework
encompasses latent space models introduced in the social network analysis
literature [12, 14, 15] as well as matrix factorization methods used in
recommender systems [8, 21, 33], up to more recent link prediction techniques
based on graph embedding algorithms [9, 29] or graph neural networks [7, 20].
We thus generically refer to all these models as latent space models. Note
that some models use only one embedding per node, which boils down to setting
$\mathbf{u}_{i}=\mathbf{v}_{i}$ for each node $i$. This is typically relevant
for undirected graphs. It is also possible to include additional information
such as node- or edge-related features, but we do not consider such extensions
here since they are rarely used in computer network monitoring applications
[1, 39].

A straightforward extension of latent space models to temporal graphs consists
in replacing fixed node embeddings with time series
$(\mathbf{u}_{i,t},\mathbf{v}_{i,t})_{t\geq{1}}$. This leads to a temporal
link predictor $h(i,j,t)=g(\mathbf{u}_{i,t},\mathbf{v}_{j,t})$, with node
embeddings typically modelled as a hidden Markov chain. Such models originated
in the social network analysis [36, 41] and recommender systems [10, 24]
communities, where they were used to model the long-term dynamics of temporal
graphs. They were also recently applied to link prediction and anomaly
detection in computer networks [19, 23, 35]. However, modelling network
dynamics by adjusting node embeddings is mostly relevant for social networks
and user–content interaction networks, where the evolution of the graph is
smooth and driven by node-specific factors (e.g., the social environment of an
individual in a social network). In contrast, activity in computer networks
varies sharply (see Section 4.1 for an illustration), and short-term dynamics
are mostly driven by global factors (e.g., office hours that apply to all
users). As a consequence, temporal link prediction methods originally designed
for social network analysis or recommender systems are not necessarily well
suited for computer network monitoring.

#### Tensor factorization and beyond.

Besides dynamic latent space models, another popular approach to temporal
graph modelling relies on Candecomp/Parafac (CP) tensor factorization [3, 4,
13], which Eren et al. [5] proposed to use for computer network monitoring.
The idea behind this approach is to represent a temporal graph as a three-mode
tensor, with modes representing the time step, the origin and the destination,
respectively. This leads to a link predictor
$h(i,j,t)=\mathbf{w}_{t}\cdot(\mathbf{u}_{i}\odot\mathbf{v}_{j})$, where
$\mathbf{w}_{t}\in\mathbb{R}^{K}$ is the embedding of time step $t$.
Interestingly, this model can be interpreted in source separation terms: each
adjacency matrix $\mathbf{A}_{t}$ is approximated as a linear combination of
rank-one matrices

$\hat{\mathbf{A}}_{t}=\sum_{k=1}^{K}w_{tk}\mathbf{U}_{k}\mathbf{V}_{k}^{\top},$

where $\mathbf{U}_{k}$ (resp. $\mathbf{V}_{k}$) is the $k$-th column of the
$N\times{K}$ embedding matrix
$\mathbf{U}=(\mathbf{u}_{1}\|\ldots\|\mathbf{u}_{N})^{\top}$ (resp.
$\mathbf{V}=(\mathbf{v}_{1}\|\ldots\|\mathbf{v}_{N})^{\top}$; $\|$ denotes
concatenation). Each of these matrices can then be interpreted as an activity
source, with the weights $w_{t1},\ldots,w_{tK}$ acting as time-varying mixing
coefficients. This corresponds to the setting illustrated in Figure 1, with
each activity source described using a very simple model (namely, rank-one
matrix factorization). Our approach, explained in the next section, builds
further upon this analogy by allowing the number of activity sources to be
different from (and typically smaller than) the embedding dimension $K$. In
other words, we make the number of activity sources smaller while using a more
expressive model to represent each source.

##  3 A Source Separation Approach to Temporal Graphs

We now move on to the description of our temporal link prediction and anomaly
detection algorithm. Section 3.1 describes our model and its inference
procedure, and we then address anomaly detection using the trained model in
Section 3.2.

###  3.1 Superposed NMF — Model Description and Inference

(a) Latent space model [23, 36]

(b) Tensor factorization [4, 5]

(c) SNMF

Figure 2: Three approaches to temporal graph modelling.

#### Model description.

Let $L\geq{1}$ be the number of activity sources, which is treated as a
hyperparameter. For each $\ell\in[L]$, we introduce origin and destination
embedding matrices
$\mathbf{U}_{\ell},\mathbf{V}_{\ell}\in\mathbb{R}_{+}^{N\times{K}}$. Then, our
link predictor is defined as

$h(i,j,t)=\sum_{\ell=1}^{L}w_{t\ell}\mathbf{u}_{i,\ell}\cdot\mathbf{v}_{j,\ell},$
(1)

where $\mathbf{u}_{i,\ell}$ and $\mathbf{v}_{j,\ell}$ denote the $i$-th and
$j$-th row of $\mathbf{U}_{\ell}$ and $\mathbf{V}_{\ell}$, respectively. The
initial mixing matrix $\mathbf{W}=(w_{t\ell})\in\mathbb{R}_{+}^{T\times{L}}$
is learned during training, and subsequent mixing coefficients are predicted
using a seasonal model (see Section 3.2).

Coming back to the source separation interpretation introduced in the previous
section, Equation 1 approximates the adjacency matrix $\mathbf{A}_{t}$ as a
weighted sum of nonnegative matrices, each of which is modelled through
nonnegative matrix factorization (NMF). We thus call this model Superposed
Nonnegative Matrix Factorization (SNMF). The only time-varying parameters are
the mixing coefficients $w_{t1},\ldots,w_{tL}$, where $L$ is typically less
than a dozen. This distinguishes SNMF from tensor factorization, as
illustrated in Figure 2: by setting the number of sources $L$ independently
from the embedding dimension $K$, we are able to allocate more parameters to
modelling each activity source without increasing the complexity of the
temporal model.

The expected benefits of a simpler model are threefold. First, it should be
more robust to noise in the training data, and thus more effective for link
prediction and anomaly detection. For instance, Figure 2a illustrates the
structure of a dynamic latent space model (i.e., a latent space model with
time-dependent node embeddings), which is much more complex and flexible than
both SNMF and tensor factorization. As we show empirically in Section 4.2,
this greater expressivity leads to lower link prediction and anomaly detection
performance when dealing with computer network activity. Secondly, simple
dynamics lead to smaller computational cost: it is indeed less expensive to
predict $L$ mixing coefficients than to adjust $2N$ node embeddings, each of
size $K$. Finally, provided our model infers meaningful activity sources, it
is more interpretable than dynamic latent space models or tensor
factorizations. We demonstrate this point empirically in Section 4.1.

#### Model inference.

Given a training sequence of adjacency matrices
$(\mathbf{A}_{1},\ldots,\mathbf{A}_{T})$, model inference boils down to
finding parameters $\mathbb{U}=(\mathbf{U}_{\ell})_{\ell=1}^{L}$,
$\mathbb{V}=(\mathbf{V}_{\ell})_{\ell=1}^{L}$, $\mathbf{W}$ minimizing

$\begin{split}J(\mathbb{U},\mathbb{V},\mathbf{W})=&\frac{1}{2}\sum_{t=1}^{T}\left\|(\mathbf{1}_{N}-\mathbf{I}_{N})\odot\left(\mathbf{A}_{t}-\sum_{\ell=1}^{L}w_{t\ell}\mathbf{U}_{\ell}\mathbf{V}_{\ell}^{\top}\right)\right\|_{\mathrm{F}}^{2}\\\
&+\lambda_{1}\left\|\mathbf{W}\right\|_{1}+\frac{\lambda_{2}}{2}\sum_{\ell=1}^{L}\left(\left\|\mathbf{U}_{\ell}\right\|_{\mathrm{F}}^{2}+\left\|\mathbf{V}_{\ell}\right\|_{\mathrm{F}}^{2}\right),\end{split}$
(2)

where $\lambda_{1},\lambda_{2}>0$ are regularization hyperparameters and
$\|\cdot\|_{1}$, $\|\cdot\|_{\mathrm{F}}$ denote the $L_{1}$ and Frobenius
norms, respectively. Note that we exclude the diagonal of the adjacency matrix
from the computation of the mean squared error as we assume there are no self-
loops. Besides, we use $L_{1}$ regularization on the mixing matrix in addition
to the usual $L_{2}$ penalty on the embedding matrices. Indeed, we want to
encourage sparsity in the mixing matrix as it reflects the domain-specific
assumption that most activity sources are intermittent: for instance, human
activity should not be observed outside office hours.

We minimize the training objective $J(\mathbb{U},\mathbb{V},\mathbf{W})$ using
a multiplicative update procedure similar to the one proposed by Lee and Seung
for standard NMF [22]. This procedure is described in Algorithm 1, and more
details on its derivation can be found in Appendix 0.A.

Data: Adjacency matrices $\mathbf{A}_{1},\ldots,\mathbf{A}_{T}$ (size
$N\times{N}$); embedding dimension $K$; number of sources $L$; regularization
hyperparameters $\lambda_{1},\lambda_{2}$.

Result: Mixing matrix $\mathbf{W}$; embedding matrix sequences $\mathbb{U}$
and $\mathbb{V}$.

Randomly initialize positive matrices
$\mathbf{W},(\mathbf{U}_{1},\ldots,\mathbf{U}_{L}),(\mathbf{V}_{1},\ldots,\mathbf{V}_{L})$;

repeat

/* Update the mixing coefficients */

foreach  _$\ell\in[L],t\in[T]$_ do

$w_{t\ell}\leftarrow
w_{t\ell}\frac{\left\langle\mathbf{A}_{t},\mathbf{U}_{\ell}\mathbf{V}_{\ell}^{\top}\right\rangle_{\mathrm{F}}}{\left\langle\mathbf{1}_{N}-\mathbf{I}_{N},\mathbf{U}_{\ell}\mathbf{V}_{\ell}^{\top}\odot\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{U}_{\ell^{\prime}}\mathbf{V}_{\ell^{\prime}}^{\top}\right\rangle_{\mathrm{F}}+\lambda_{1}}$;

end foreach

/* Update the origin embedding matrices */

foreach  _$\ell\in[L]$_ do

$\mathbf{M}\leftarrow\sum_{t=1}^{T}w_{t\ell}\left(\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\left(\mathbf{1}_{N}-\mathbf{I}_{N}\right)\odot\left(\mathbf{U}_{\ell^{\prime}}\mathbf{V}_{\ell^{\prime}}^{\top}\right)\right)\mathbf{V}_{\ell}+\lambda_{2}\mathbf{U}_{\ell}$;

$\mathbf{U}_{\ell}\leftarrow\mathbf{U}_{\ell}\odot\left(\sum_{t=1}^{T}w_{t\ell}\mathbf{A}_{t}\mathbf{V}_{\ell}\right)\div\mathbf{M}$;

end foreach

/* Update the destination embedding matrices */

foreach  _$\ell\in[L]$_ do

$\mathbf{M}\leftarrow\sum_{t=1}^{T}w_{t\ell}\left(\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\left(\mathbf{1}_{N}-\mathbf{I}_{N}\right)\odot\left(\mathbf{U}_{\ell^{\prime}}\mathbf{V}_{\ell^{\prime}}^{\top}\right)^{\top}\right)\mathbf{U}_{\ell}+\lambda_{2}\mathbf{V}_{\ell}$;

$\mathbf{V}_{\ell}\leftarrow\mathbf{V}_{\ell}\odot\left(\sum_{t=1}^{T}w_{t\ell}\mathbf{A}_{t}^{\top}\mathbf{U}_{\ell}\right)\div\mathbf{M}$;

end foreach

until  _convergence_ ;

return _
$\mathbf{W},\,\mathbb{U}=(\mathbf{U}_{\ell})_{\ell=1}^{L},\,\mathbb{V}=(\mathbf{V}_{\ell})_{\ell=1}^{L}$
_

Algorithm 1 Multiplicative update procedure for SNMF.

###  3.2 Anomaly Detection Using the Trained Model

Given a new temporal edge $(i,j,t)\in\mathcal{E}$ (with $t>T$), computing its
anomaly score $h(i,j,t)$ as defined in Equation 1 requires knowing the weights
$w_{t1},\ldots,w_{tL}$, which are not readily available at test time. We thus
adopt a seasonal modelling approach: letting $\tau\geq{1}$ denote the period
of the model, we predict $\mathbf{w}_{t}$ as

$\hat{\mathbf{w}}_{t}=\frac{1}{|\mathcal{T}(t,\tau)|}\sum_{t^{\prime}\in\mathcal{T}(t,\tau)}\mathbf{w}_{t^{\prime}},\quad\text{where
}\mathcal{T}(t,\tau)=\\{t^{\prime}\in[t-1]:t\equiv t^{\prime}\mod\tau\\}.$

Typically, a reasonable value for the period $\tau$ is one week. This is the
value we use in our experiments (see Section 4.2).

In addition, after fully observing and scoring a new graph $\mathcal{G}_{t}$,
we learn its vector of mixing coefficients $\mathbf{w}_{t}$ by minimizing the
objective function of Equation 2 for the adjacency matrix $\mathbf{A}_{t}$,
with $\mathbb{U}$ and $\mathbb{V}$ fixed. This allows the model to
continuously expand its set of historical mixing coefficients, which should in
turn lead to more reliable predictions using the seasonal model. Note that we
have not yet implemented any forgetting mechanism, which could be useful when
deploying the model for a long time (see Section 5 for further discussion of
long-term dynamics and associated challenges).

##  4 Experiments

In order to evaluate our model and test the underlying hypothesis on the
nature of the short-term dynamics of computer network activity, we perform two
kinds of experiments. In Section 4.1, we train SNMF on a small, well-
documented dataset, which allows us to interpret and qualitatively evaluate
the inferred activity sources. We then evaluate the performance of our model
on link prediction and anomaly detection tasks for a larger, more complex and
realistic dataset, and compare it with several state-of-the-art baselines in
Section 4.2.

Table 1: Datasets used in our experiments. Dataset | Nodes | Time steps | Edges
---|---|---|---
Min. | Med. | Max.
VAST | 1427 | 339 | 0 | 463 | 15,087
LANL | 12,702 | 720 | 15,147 | 33,980.5 | 59,944

###  4.1 Qualitative Evaluation — Network Traffic Analysis

First of all, we need to assess whether our hypothesis on the dynamics of
computer network activity is correct: are there such things as activity
sources that can be recovered by our model? One way to answer this question is
to inspect the inferred sources and mixing coefficients and confront them with
our assumptions.

#### Dataset description.

To perform this qualitative evaluation, we use a dataset that was originally
generated for the Mini-Challenge 3 of the VAST 2013 competition [40]. This
dataset contains two weeks of simulated network traffic involving an
enterprise network and various external hosts. Note that while it contains
some attacks, our main goal here is not to detect them: we only want to
determine whether our model yields a meaningful decomposition of the temporal
graph. Moreover, even though the synthetic nature of the data raises the
question of their realism, it also provides us with a detailed knowledge of
the structure of the network, the functional role of each host and the types
of activity going on at any time. This allows us to confront the patterns
extracted by SNMF with the ground truth.

We divide the dataset into one hour windows. For each window, we build a graph
whose nodes are the hosts, identified by their IP addresses. The temporal edge
$(i,j,t)$ exists if $i$ has initiated at least one network communication with
$j$ in time window $t$. See Table 1 for a description of the obtained temporal
graph.

#### Experiment results.

We train SNMF with the following hyperparameters: embedding dimension $K=5$,
number of activity sources $L=4$, regularization coefficients
$\lambda_{1}=10^{-3}$ and $\lambda_{2}=10^{-5}$. For each source $\ell$, we
then extract the embedding matrices $\mathbf{U}_{\ell},\mathbf{V}_{\ell}$ and
perform $k$-means clustering on their concatenation
$(\mathbf{U}_{\ell}\|\mathbf{V}_{\ell})\in\mathbb{R}_{+}^{N\times(2K)}$. The
number of clusters is chosen by maximizing the silhouette score [16]. Finally,
we compute the predicted adjacency matrix
$\mathbf{U}_{\ell}\mathbf{V}_{\ell}^{\top}$ and turn it into a graph by
creating an edge $(i,j)$ for each coefficient
$\mathbf{u}_{i,\ell}\cdot\mathbf{v}_{j,\ell}$ greater than a threshold
$\theta_{\ell}$. This threshold is set through manual analysis of the
distribution of the coefficients.

Figure 3:  Activity sources found in the VAST dataset (top) and evolution of
their mixing coefficients over time (bottom). The intensity of each link
reflects the total number of edges between the two clusters.

The obtained clustered graphs are displayed along with the inferred mixing
coefficients in Figure 3. The most striking observation is that the timeline
can be divided into two clearly different parts, separated by a period of
inactivity. This is consistent with the underlying scenario: after facing
several attacks in the first week, the system administrator takes the network
offline for three days in order to strengthen its defenses. In particular, an
intrusion prevention system (IPS) is deployed and subsequently blocks most of
the malicious traffic. These changes are clearly reflected in the mixing
coefficients, which highlights the ability of SNMF to identify relevant
patterns in the data. For instance, source 1, which contains inbound traffic
from external hosts into the enterprise network, disappears in the second
week. In contrast, sources 3 and 4, which exhibit a simple client/server
structure (traffic from internal workstations to internal servers and external
hosts), are dominant after the activity gap. Moreover, sources 3 and 4 exhibit
a clear seasonal pattern, suggesting office-related activity, while source 1
appears to be more of a background noise partly caused by traffic coming from
outside the enterprise network. Overall, Figure 3 thus demonstrates the
existence of distinct activity sources that can be recovered to some extent by
SNMF.

###  4.2 Quantitative Evaluation — Temporal Link Prediction and Anomaly
Detection

We now aim to evaluate the performance of our model for temporal link
prediction and malicious activity detection in computer networks. Indeed,
while our primary goal is to detect edges corresponding to malicious behavior,
the ability of our model to predict future edges also reflects its conformity
with the true nature of the short-terms dynamics of computer network activity.

#### Dataset description.

This quantitative evaluation is carried out using the "Comprehensive, Multi-
Source Cyber-Security Events" dataset released by the Los Alamos National
Laboratory [17, 18]. This dataset contains event logs collected from a real-
world enterprise network over 58 days. Moreover, a red team exercise took
place during this time frame, meaning that security experts tried to breach
the network in order to assess its security. The remote logons performed by
the red team are labelled, and we thus aim to detect them. We consider the
first 30 days of data and divide them into one hour windows. For each window,
we build a graph whose nodes are the hosts of the enterprise network. The edge
$(i,j,t)$ exists if there is at least one remote logon from $i$ to $j$ in
window $t$. See Table 1 for a description of the obtained temporal graph. We
use the first seven days for training, the eighth day as a validation set and
the next 22 days for evaluation.

#### Evaluation protocol and baselines.

Besides detecting malicious logons, we also evaluate the temporal link
prediction performance of SNMF. More specifically, we adopt the evaluation
protocol of Poursafaei et al. [30], who propose to detect positive temporal
edges among three types of negative edges:

  * •

random negative edges, obtained by randomly rewiring positive edges;

  * •

historical negative edges, sampled from the set of edges observed during
training but not in the current time window;

  * •

inductive negative edges, sampled from the set of edges observed in the test
set but not in the training set, and not in the current time window.

For each time window $t$ in the test set, we generate as many negative edges
of each type as there are positive edges in the graph $\mathcal{G}_{t}$.

We evaluate SNMF against several baselines. First, two dynamic models are
included: Poisson tensor factorization (PTF [5]) and the bilinear mixed-
effects model proposed by Lee et al. (BME [23]). These two methods represent
the two main approaches to dynamic latent space modelling identified in
Section 2.2: adjusting all node embeddings using global time-dependent
parameters, or allowing each node embedding to evolve on its own. We also
include two static models: hierarchical Poisson factorization (HPF [35]) and
the GL-GV model introduced by Bowman et al. [1], which relies on the node2vec
graph embedding algorithm [9]. Next up, SedanSpot [6] is a generic anomaly
detection algorithm for temporal graphs: unlike other baselines, it was not
specifically proposed as a computer network monitoring method. Finally, the
two naive temporal link prediction methods introduced by Poursafaei et al.
[30] are also evaluated: EdgeBank∞ memorizes all edges seen so far and
considers them normal, while its variant EdgeBankw only memorizes edges seen
in the recent past. The length of the memorized window for EdgeBankw is set to
one week in our experiments. As for latent space models, their hyperparameters
are set by maximizing the average area under the ROC curve (AUC) for the three
link prediction tasks on the validation set. In particular, the embedding
dimension $K$ is selected from the set $\\{10,20,30,40,50\\}$. For SNMF, the
number of sources $L$ can take the values 2, 3, 4 or 5, and the embedding
dimension for each source is set to $\lfloor{K}/{L}\rfloor$ for fair
comparison with the baselines.

Table 2: Mean and standard deviation of the area under the ROC curve for each task.  Method | Anomaly | Random | Historical | Inductive
---|---|---|---|---
SNMF | 99.1$\pm$0.1 | 98.4$\pm$0.0 | 76.9$\pm$0.2 | 98.2$\pm$0.1
PTF [5] | 97.6$\pm$0.9 | 98.6$\pm$0.0 | 68.5$\pm$0.2 | 96.7$\pm$0.3
BME [23] | 90.3$\pm$0.2 | 98.5$\pm$0.0 | 73.3$\pm$0.0 | 96.2$\pm$0.0
HPF [35] | 97.7$\pm$0.3 | 99.1$\pm$0.0 | 69.6$\pm$0.1 | 97.5$\pm$0.0
GL-GV [1] | 87.0$\pm$2.4 | 95.8$\pm$1.0 | 61.2$\pm$0.5 | 74.6$\pm$1.9
SedanSpot [6] | 63.6$\pm$7.3 | 51.2$\pm$2.2 | 54.7$\pm$1.4 | 53.2$\pm$2.7
EdgeBank∞ [30] | 96.2$\pm$0.0 | 97.2$\pm$0.0 | 56.0$\pm$0.0 | 98.3$\pm$0.0
EdgeBankw [30] | 96.0$\pm$0.0 | 97.0$\pm$0.0 | 58.0$\pm$0.0 | 98.1$\pm$0.0
Figure 4: Mean and standard deviation of the NDCG@1% score (anomaly detection
task).

#### Evaluation results.

Table 2 reports the mean and standard deviation over 10 runs of the AUC for
each task and each evaluated method. First of all, SNMF outperforms all
baselines on the anomaly detection task, which is especially important as the
main goal of our work is to detect malicious activity. Note, however, that the
highly imbalanced nature of the anomaly detection task makes the AUC
insufficient as an evaluation metric. Indeed, in a real-world computer network
monitoring setting, only an extremely low false positive rate would be
acceptable, whereas the AUC is computed on the whole ROC curve. We thus also
compute the normalized discounted cumulative gain for the top 1% of the
anomaly ranking (NDCG@1%), which accounts for the greater importance of the
highest-ranking temporal edges. The results, displayed in Figure 4, confirm
that SNMF outperforms all baselines.

Regarding link prediction tasks, SNMF performs best on historical negative
edges, which also sustains our initial hypothesis on short-term dynamics of
computer network activity. Indeed, distinguishing positive edges from
historical negatives requires not only memorizing edges seen during training,
but also accurately modelling the corresponding temporal patterns. In
particular, the fact that SNMF outperforms PTF confirms that a small number
$L<K$ of activity sources is sufficient to model these patterns. As for
inductive negative edges, while SNMF is the best-performing latent space
model, none of these models actually outperforms the naive baseline EdgeBank∞.
The most likely cause for this is that new edges observed at test time in the
LANL dataset are noisy events that do not reoccur frequently. As a
consequence, distinguishing them from positive edges is easy, and modelling
temporal patterns does not bring any improvement. A similar phenomenon occurs
for random negatives, which explains why static matrix factorization performs
best on this task.

Figure 5: Mean and standard deviation of the area under the ROC curve for each
task, as a function of the number of sources $L$ and the total dimension
$L\cdot{K}$.

#### Sensitivity analysis.

Finally, in order to study the respective influence of the embedding dimension
$K$ and the number of sources $L$ on the accuracy of SNMF for each task, we
display the AUC as a function of $L$ and $K$ in Figure 5. The obtained plot
for the anomaly detection task is noteworthy: increasing $K$ degrades
performance, and this decrease gets steeper as $L$ increases. In other words,
simpler models perform better, which results from the somewhat overt nature of
malicious activity in the LANL dataset: the anomalous logons correspond to
edges that are never observed as part of normal activity, which makes them
easy to detect even without a deep understanding of the temporal dynamics.
However, such knowledge would be useful for detecting stealthier attacks.

As for link prediction, increasing $L$ up to 4 improves performance in the
historical setting, but this improvement then plateaus for $L=5$. This
confirms that a few activity sources are enough, thereby demonstrating the
relevance of our approach. Finally, in the random and inductive settings,
simply memorizing frequently reoccurring edges leads to satisfactory
performance, which explains why increasing $K$ improves the AUC while
increasing $L$ slightly degrades it.

##  5 Discussion and Future Work

We now discuss limitations and potential extensions of our work, identifying
several leads for future research.

#### Extension to more complex data and models.

One important characteristic of computer network monitoring is that it
produces rich and complex data. While we focus on communications between hosts
and their temporal component, one important avenue of research deals with
jointly modelling the many other facets of these communications, such as the
associated user accounts or network protocols [5]. Note, however, that the
core of our approach could easily be extended in this direction: typically,
instead of predicting temporal edges $(i,j,t)$, we could predict hyperedges
$(v_{1},\ldots,v_{d},t)$ using the same source separation approach, provided
that the model for each activity source is adapted to such data.

In addition, we represent network activity as a sequence of graphs, each of
which stands for a given time window. However, this temporal aggregation step
leads to some information loss, and it could thus be interesting to directly
model individual interactions between hosts instead. The stream of
interactions could for instance be modelled as a marked point process, in a
way similar to previous contributions on recommender systems [37] and computer
network monitoring [34]. The source separation approach would then consist in
modelling the stochastic intensity of this point process as a time-dependent
linear combination of simpler functions representing the activity sources.

#### Better prediction of the mixing coefficients.

As described in Section 3.2, we use a simple seasonal model to predict the
vector of mixing coefficients $\mathbf{w}_{t}$ at test time. However, more
flexible and sophisticated methods could also be designed. For instance, an
autoregressive model could be used to learn the temporal dynamics of
$\mathbf{w}_{t}$. More generally, the sequence $(\mathbf{w}_{t})_{t\geq{1}}$
can be seen as a low-dimensional multivariate time series, and any time series
forecasting method can thus be used to predict the mixing coefficients.
Another possible approach, similar to the work of Gutflaish et al. [11], is to
learn to predict $\mathbf{w}_{t}$ based on some explicit features of the
current time step (e.g. hour of day, day of week, number of observed events).

#### Combining short- and long-term dynamics.

Finally, we emphasize that our model only addresses short-term dynamics, i.e.,
variations in the observed activity within a given day or week. However,
computer networks are also prone to long-term dynamics, which are driven by
two main causes. First, existing nodes can change their behavior: for
instance, a server can start hosting a new application, thus receiving more
traffic from other hosts. Secondly, the node set can evolve over time as new
hosts are added to the network. These two problems are well-known in the
recommender systems literature, where they are typically referred to as
concept drift and cold start, respectively. While concept drift can be
mitigated by partially retraining the model on a regular basis [10],
addressing the cold start problem typically involves initializing a new node’s
embedding based on those of old nodes whose behavior is similar [11]. Both of
these procedures could easily be applied to our model.

##  6 Conclusion

We propose a simple model for the short-term dynamics of communications within
a computer network. This model builds upon the idea that the observed activity
results from the superposition of a small number of activity sources, and that
short-term dynamics are mainly caused by the temporal variation of the
respective weights of these sources. Both qualitative and quantitative
evaluation support this hypothesis. In particular, our model outperforms
state-of-the-art baselines at detecting anomalous edges resulting from
malicious behavior.

These results lead to the following insights. First, computer network activity
exhibits specific temporal dynamics, which differ significantly from those
observed in social networks or user–content interaction streams. As a
consequence, temporal link prediction models designed for such data should not
be applied directly to computer networks. Secondly, these dynamics, despite
their peculiarity, appear to be reasonably simple and predictable. Building
upon simple models might thus be the best way forward.

##  Appendix 0.A Derivation of the Inference Procedure

Training the SNMF model requires minimizing the function $J$ defined in
Equation 2, which amounts to a minimization problem subject to nonnegativity
constraints. The Karush-Kuhn-Tucker (KKT) conditions [2] translate to

$\left\\{\begin{array}[]{lll}\frac{\partial J}{\partial w_{t\ell}}\geq
0&\quad\text{and}\quad w_{t\ell}\frac{\partial J}{\partial w_{t\ell}}=0&\text{
for all }t\in[T],\,\ell\in[L]\\\ \frac{\partial J}{\partial u_{ik,\ell}}\geq
0&\quad\text{and}\quad u_{ik,\ell}\frac{\partial J}{\partial
u_{ik,\ell}}=0&\text{ for all }i\in[N],\,k\in[K],\,\ell\in[L]\\\
\frac{\partial J}{\partial v_{jk,\ell}}\geq 0&\quad\text{and}\quad
v_{jk,\ell}\frac{\partial J}{\partial v_{jk,\ell}}=0&\text{ for all
}j\in[N],\,k\in[K],\,\ell\in[L]\\\ \end{array}\right.$ (3)

where $u_{ik,\ell}$ (resp. $v_{jk,\ell}$) is the coefficient of
$\mathbf{U}_{\ell}$ (resp. $\mathbf{V}_{\ell}$) at position $(i,k)$ (resp.
$(j,k)$). The partial derivatives are

$\left\\{\begin{array}[]{ll}\frac{\partial J}{\partial
w_{t\ell}}&=-\sum_{1\leq i\neq j\leq
N}\mathbf{u}_{i,\ell}\cdot\mathbf{v}_{j,\ell}\left(A_{ij,t}-\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{u}_{i,\ell^{\prime}}\cdot\mathbf{v}_{j,\ell^{\prime}}\right)+\lambda_{1}\\\
\frac{\partial J}{\partial u_{ik,\ell}}&=-\sum_{t=1}^{T}\sum_{j\neq
i}w_{t\ell}v_{jk,\ell}\left(A_{ij,t}-\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{u}_{i,\ell^{\prime}}\cdot\mathbf{v}_{j,\ell^{\prime}}\right)+\lambda_{2}u_{ik,\ell}\\\
\frac{\partial J}{\partial v_{jk,\ell}}&=-\sum_{t=1}^{T}\sum_{i\neq
j}w_{t\ell}u_{ik,\ell}\left(A_{ij,t}-\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{u}_{i,\ell^{\prime}}\cdot\mathbf{v}_{j,\ell^{\prime}}\right)+\lambda_{2}v_{jk,\ell}\\\
\end{array}\right.$ (4)

Plugging Equation 4 into the equality conditions of Equation 3 then yields

$\left\\{\begin{array}[]{lll}w_{t\ell}\sum_{1\leq i\neq j\leq
N}A_{ij,t}\mathbf{u}_{i,\ell}\cdot\mathbf{v}_{j,\ell}&=&w_{t\ell}\sum_{1\leq
i\neq j\leq
N}\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{u}_{i,\ell^{\prime}}\cdot\mathbf{v}_{j,\ell^{\prime}}+\lambda_{1}\\\
u_{ik,\ell}\sum_{t=1}^{T}w_{t\ell}\sum_{j\neq
i}A_{ij,t}v_{jk,\ell}&=&u_{ik,\ell}(\sum_{t=1}^{T}w_{t\ell}\sum_{j\neq
i}v_{jk,\ell}\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{u}_{i,\ell^{\prime}}\cdot\mathbf{v}_{j,\ell^{\prime}}\\\
&&+\lambda_{2}u_{ik,\ell})\\\ v_{jk,\ell}\sum_{t=1}^{T}w_{t\ell}\sum_{i\neq
j}A_{ij,t}u_{ik,\ell}&=&v_{jk,\ell}(\sum_{t=1}^{T}w_{t\ell}\sum_{i\neq
j}u_{ik,\ell}\sum_{\ell^{\prime}=1}^{L}w_{t\ell^{\prime}}\mathbf{u}_{i,\ell^{\prime}}\cdot\mathbf{v}_{j,\ell^{\prime}}\\\
&&+\lambda_{2}v_{jk,\ell})\\\ \end{array}\right.$

which leads to the multiplicative updates of Algorithm 1.

## References

  * [1] Bowman, B., Laprade, C., Ji, Y., Huang, H.H.: Detecting lateral movement in enterprise computer networks with unsupervised graph AI. In: RAID (2020)
  * [2] Boyd, S., Vandenberghe, L.: Convex optimization. Cambridge University Press (2004)
  * [3] Carroll, J.D., Chang, J.J.: Analysis of individual differences in multidimensional scaling via an n-way generalization of "Eckart-Young" decomposition. Psychometrika 35(3), 283–319 (1970)
  * [4] Dunlavy, D.M., Kolda, T.G., Acar, E.: Temporal link prediction using matrix and tensor factorizations. ACM Trans. Knowl. Discov. Data 5(2), 1–27 (2011)
  * [5] Eren, M.E., Moore, J.S., Alexandro, B.S.: Multi-dimensional anomalous entity detection via Poisson tensor factorization. In: ISI (2020)
  * [6] Eswaran, D., Faloutsos, C.: Sedanspot: Detecting anomalies in edge streams. In: ICDM (2018)
  * [7] Gilmer, J., Schoenholz, S.S., Riley, P.F., Vinyals, O., Dahl, G.E.: Neural message passing for quantum chemistry. In: ICML (2017)
  * [8] Gopalan, P., Hofman, J.M., Blei, D.M.: Scalable recommendation with hierarchical poisson factorization. In: UAI (2015)
  * [9] Grover, A., Leskovec, J.: node2vec: Scalable feature learning for networks. In: KDD (2016)
  * [10] Gultekin, S., Paisley, J.: A collaborative Kalman filter for time-evolving dyadic processes. In: ICDM (2014)
  * [11] Gutflaish, E., Kontorovich, A., Sabato, S., Biller, O., Sofer, O.: Temporal anomaly detection: calibrating the surprise. In: AAAI (2019)
  * [12] Handcock, M.S., Raftery, A.E., Tantrum, J.M.: Model-based clustering for social networks. J. R. Stat. Soc. A Stat. 170(2), 301–354 (2007)
  * [13] Harshman, R.A.: Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multimodal factor analysis. UCLA Work. Pap. Phon. 16, 1–84 (1970)
  * [14] Hoff, P.D.: Bilinear mixed-effects models for dyadic data. J. Am. Stat. Assoc. 100(469), 286–295 (2005)
  * [15] Hoff, P.D., Raftery, A.E., Handcock, M.S.: Latent space approaches to social network analysis. J. Am. Stat. Assoc. 97(460), 1090–1098 (2002)
  * [16] Kaufman, L., Rousseeuw, P.J.: Finding groups in data: an introduction to cluster analysis. Wiley (2009)
  * [17] Kent, A.D.: Comprehensive, multi-source cyber-security events. Los Alamos National Laboratory (2015). https://doi.org/10.17021/1179829
  * [18] Kent, A.D.: Cybersecurity data sources for dynamic network research. In: Dynamic Networks in Cybersecurity. Imperial College Press (2015)
  * [19] King, I.J., Huang, H.H.: EULER: Detecting network lateral movement via scalable temporal link prediction. In: NDSS (2022)
  * [20] Kipf, T.N., Welling, M.: Variational graph auto-encoders. arXiv preprint arXiv:1611.07308 (2016)
  * [21] Koren, Y., Bell, R.M., Volinsky, C.: Matrix factorization techniques for recommender systems. Computer 42(8), 30–37 (2009)
  * [22] Lee, D.D., Seung, H.S.: Algorithms for non-negative matrix factorization. In: NeurIPS (2000)
  * [23] Lee, W., McCormick, T.H., Neil, J., Sodja, C., Cui, Y.: Anomaly detection in large-scale networks with latent space models. Technometrics 64(2), 241–252 (2022)
  * [24] Lu, Z., Agarwal, D., Dhillon, I.S.: A spatio-temporal approach to collaborative filtering. In: RecSys (2009)
  * [25] Marchette, D.J.: Computer intrusion detection and network monitoring: a statistical viewpoint. Springer (2001)
  * [26] Marjan, M., Zaki, N., Mohamed, E.A.: Link prediction in dynamic social networks: A literature review. In: CIST (2018)
  * [27] Neil, J., Hash, C., Brugh, A., Fisk, M., Storlie, C.B.: Scan statistics for the online detection of locally anomalous subgraphs. Technometrics 55(4), 403–414 (2013)
  * [28] Neil, J., Uphoff, B., Hash, C., Storlie, C.: Towards improved detection of attackers in computer networks: New edges, fast updating, and host agents. In: ISRCS (2013)
  * [29] Perozzi, B., Al-Rfou, R., Skiena, S.: Deepwalk: Online learning of social representations. In: KDD (2014)
  * [30] Poursafaei, F., Huang, S., Pelrine, K., Rabbany, R.: Towards better evaluation for dynamic link prediction. In: NeurIPS (2022)
  * [31] Price-Williams, M., Turcotte, M.J., Heard, N.: Time of day anomaly detection. In: EISIC (2018)
  * [32] Rabiu, I., Salim, N., Da’u, A., Osman, A.: Recommender system based on temporal models: a systematic review. Appl. Sci. 10(7), 2204 (2020)
  * [33] Salakhutdinov, R., Mnih, A.: Probabilistic matrix factorization. In: NeurIPS (2007)
  * [34] Sanna Passino, F., Heard, N.A.: Mutually exciting point process graphs for modeling dynamic networks. J. Comput. Graph. Stat. pp. 1–15 (2022)
  * [35] Sanna Passino, F., Turcotte, M.J., Heard, N.A.: Graph link prediction in computer networks using Poisson matrix factorisation. Ann. Appl. Stat. 16(3), 1313–1332 (2022)
  * [36] Sarkar, P., Moore, A.W.: Dynamic social network analysis using latent space models. In: NeurIPS (2005)
  * [37] Shang, J., Sun, M.: Local low-rank Hawkes processes for temporal user-item interactions. In: ICDM (2018)
  * [38] Turcotte, M., Heard, N., Neil, J.: Detecting localised anomalous behaviour in a computer network. In: IDA (2014)
  * [39] Turcotte, M., Moore, J., Heard, N., McPhall, A.: Poisson factorization for peer-based anomaly detection. In: ISI (2016)
  * [40] Whiting, M., Cook, K., Paul, C., Whitley, K., Grinstein, G., Nebesh, B., Liggett, K., Cooper, M., Fallon, J.: VAST challenge 2013: Situation awareness and prospective analysis. In: VAST (2013)
  * [41] Zhu, L., Guo, D., Yin, J., Ver Steeg, G., Galstyan, A.: Scalable temporal latent space inference for link prediction in dynamic social networks. IEEE Trans. Knowl. Data Eng. 28(10), 2765–2777 (2016)