# Supplemental Material: Observation of wave-packet branching through an
engineered conical intersection

Christopher S. Wang<EMAIL_ADDRESS>Nicholas E. Frattini  Benjamin
J. Chapman  Shruti Puri  S. M. Girvin  Michel H. Devoret  Robert J. Schoelkopf
Departments of Physics and Applied Physics, Yale University, New Haven, CT
06511, USA.  Yale Quantum Institute, Yale University, New Haven, CT 06520,
USA.

##  I Constructing a model reaction Hamiltonian containing a conical
intersection

In this section, we describe the general formalisms behind obtaining various
model and/or ab-initio molecular Hamiltonians that involve strong vibronic
coupling. We begin with a brief review of adiabatic potential energy surfaces
and highlight the difficulties that arise in the vicinity of conical
intersections. This motivates the use of diabatic electronic states, which
form the basis of the Hamiltonians that we consider Domcke2004 .

The standard molecular Hamiltonian is:

$\hat{H}_{\textrm{mol}}=\hat{T}_{n}(\partial_{\vec{R}})+\underbrace{\hat{T}_{e}(\partial_{\vec{r}})+\hat{V}(\vec{r},\vec{R})}_{\hat{H}_{e}}$
(1)

where $\hat{T}$ and $\hat{V}$ correspond to kinetic and Coulomb potential
energies, respectively, and the subscripts $n$ and $e$ denote nuclei and
electrons. $\vec{r}$ and $\vec{R}$ represents the positions of the electrons
and nuclei, respectively.

The conventional approach based on the Born-Oppenheimer approximation begins
with noting that the electron masses $m_{e}$ are much smaller than that of the
nuclear masses $m_{N}$. This motivates momentarily dropping
$\hat{T}_{n}(\partial_{\vec{R}})\propto 1/m_{N}$ altogether, leaving behind a
reduced Hamiltonian of the electrons only with parametric dependence on the
nuclear coordinates $\vec{R}$. By choosing an ansatz for the total molecular
wave-function
$\Psi(\vec{r},\vec{R})=\Sigma_{i}\varphi_{i}(\vec{r},\vec{R})\chi_{i}(\vec{R})$,
we get a reduced electronic Schrodinger equation:

$\hat{H}_{e}\varphi_{i}(\vec{r},\vec{R})=E_{i}(\vec{R})\varphi_{i}(\vec{r},\vec{R})$
(2)

where $E_{i}(\vec{R})$ is the potential energy surface for the $i^{th}$
electronic eigenstate. Here, we have identified a complete set of adiabatic
electronic eigenfunctions $\varphi_{i}(\vec{r},\vec{R})$. Returning to the
full Schrodinger equation

$\hat{H}_{\textrm{mol}}\Psi(\vec{r},\vec{R})=(\hat{T}_{n}(\partial_{\vec{R}})+\hat{H}_{e})\sum_{i}\varphi_{i}(\vec{r},\vec{R})\chi_{i}(\vec{R}),$
(3)

we can obtain a reduced equation for the nuclear motion by integrating over a
complete set of adiabatic electronic eigenfunctions
$\varphi^{*}_{j}(\vec{r},\vec{R})$. This brings out terms such as:

$\frac{\bra{\varphi_{j}(\vec{r},\vec{R})}\hat{\nabla}_{R}\hat{H}_{e}\ket{\varphi_{i}(\vec{r},\vec{R})}}{E_{j}(\vec{R})-E_{i}(\vec{R})}$
(4)

which arise from the application of
$\hat{T}_{n}(\partial_{\vec{R}})\propto-\frac{1}{2m_{N}}\hat{\nabla}^{2}_{R}$
on $\varphi_{i}(\vec{r},\vec{R})$. These are commonly referred to as non-
adiabatic coupling terms (NACT) in the literature. It is clear that in the
vicinity of a conical intersection, these terms diverge as the denominator
becomes very small and the adiabatic electronic basis fails to be an
appropriate basis for calculations and analyses.

Given the aforementioned issue, one can consider instead a diabatic electronic
basis $\phi_{k}(\vec{r})$ such that the molecular wave-function can be
expressed as:

$\Psi(\vec{r},\vec{R})=\sum_{k}\phi_{k}(\vec{r})\chi^{\prime}_{k}(\vec{R})$
(5)

where the diabatic states are, by definition, diagonal in the nuclear kinetic
energy operator. Off-diagonal couplings between electronic states must of
course exist, but now they arise via the potential
$\bra{\phi_{j}}\hat{V}(\vec{r},\vec{R})\ket{\phi_{i}}$ and do not involve
wave-function derivatives.

This forms the basis for a general vibronic coupling Hamiltonian:

$\hat{H}_{\textrm{vc}}=\sum_{n}\ket{\phi_{n}}(\hat{T}_{n}+W_{nn}(\vec{R}))\bra{\phi_{n}}+\sum_{n\neq
m}\ket{\phi_{n}}W_{nm}(\vec{R})\bra{\phi_{m}}$ (6)

In our experiment, we consider a minimal model for a two-dimensional linear
vibronic coupling (LVC) Hamiltonian, where we have:

$\displaystyle\hat{T}_{n}$
$\displaystyle=\frac{\hat{p}^{2}_{x}}{2m_{x}}+\frac{\hat{p}^{2}_{y}}{2m_{y}}$
(7) $\displaystyle\hat{W}$
$\displaystyle=\left(\begin{array}[]{cc}\frac{1}{2}\omega^{2}_{x}\hat{x}^{2}+\frac{1}{2}\omega^{2}_{y}\hat{y}^{2}+g_{x}\hat{x}&g_{y}\hat{y}\\\
g_{y}\hat{y}&\frac{1}{2}\omega^{2}_{x}\hat{x}^{2}+\frac{1}{2}\omega^{2}_{y}\hat{y}^{2}-g_{x}\hat{x}\end{array}\right)$
(10)

for two normal coordinates $\hat{x}$ and $\hat{y}$. By re-casting these
coordinates into creation and annihilation operators
$\hat{x}\propto(\hat{a}+\hat{a}^{\dagger}),\hat{y}\propto(\hat{b}+\hat{b}^{\dagger})$,
we arrive at Eq. (1) of the main text.

##  II Engineering the conical intersection interaction

Here, we show how we engineer the reaction Hamiltonian (Eq. 1 of the main
text). As discussed in the main text, this Hamiltonian consists of two
simultaneous conditional displacement interactions of a single qubit to two
different cavity modes, where the qubit coupling axes are orthogonal. To
simplify the derivation, we will begin by focusing on how we enact just one of
these interactions. As we will see, adding additional interactions to other
cavity modes is relatively straightforward, and the qubit coupling axis is
freely adjustable in the effective $x-y$ plane.

We expand upon the derivation provided by SHG2016  by incorporating the finite
anharmonicity of the transmon mode. We begin with the static Hamiltonian of a
transmon mode $\hat{q}$ dispersively coupled to a cavity mode $\hat{c}$:

$\hat{H}_{\textrm{static}}/\hbar=\omega_{c}\hat{c}^{\dagger}\hat{c}+\omega_{q}\hat{q}^{\dagger}\hat{q}-\frac{\alpha_{q}}{2}\hat{q}^{\dagger}\hat{q}^{\dagger}\hat{q}\hat{q}-\chi\hat{c}^{\dagger}\hat{c}\hat{q}^{\dagger}\hat{q}$
(11)

where $\alpha_{q}$ is the transmon anharmonicity and $\chi$ is the dispersive
shift. At a high level, we will see that the conditional displacement
interaction arises by transforming the cross-Kerr interaction between the
transmon and the cavity. Thus, our approach will be to consider how driving
each mode transforms the static interaction. Specifically, we drive the system
with one tone coupled to the transmon and two coupled to the cavity:

$\hat{H}_{d}/\hbar=2\varepsilon_{R}\textrm{cos}[(\omega_{q}+\Delta_{R})t](\hat{q}+\hat{q}^{\dagger})-2i\varepsilon_{1}\textrm{sin}(\omega_{1}t+\varphi_{1})(\hat{c}-\hat{c}^{\dagger})-2i\varepsilon_{2}\textrm{sin}(\omega_{2}t+\varphi_{2})(\hat{c}-\hat{c}^{\dagger})$
(12)

such that the full system Hamiltonian is described by
$\hat{H}=\hat{H}_{\textrm{static}}+\hat{H}_{d}$. For convenience, we re-group
the terms such that we can write
$\hat{H}=\hat{H}_{q}(\hat{q},\hat{q}^{\dagger})+\hat{H}_{c}(\hat{c},\hat{c}^{\dagger})+\hat{H}_{\textrm{int}}$,
where
$\hat{H}_{\textrm{int}}/\hbar=-\chi\hat{c}^{\dagger}\hat{c}\hat{q}^{\dagger}\hat{q}$.
We first go into the rotating frame of the transmon drive via
$\hat{H}\rightarrow\hat{U}\hat{H}\hat{U}^{\dagger}+i\dot{\hat{U}}\hat{U}^{\dagger}$,
where $\hat{U}=e^{i(\omega_{q}+\Delta_{R})t\hat{q}^{\dagger}\hat{q}}$:

$\hat{H}/\hbar=\underbrace{-\Delta_{R}\hat{q}^{\dagger}\hat{q}-\frac{\alpha_{q}}{2}\hat{q}^{\dagger}\hat{q}^{\dagger}\hat{q}\hat{q}+\varepsilon_{R}(\hat{q}+\hat{q}^{\dagger})}_{\hat{H}_{\textrm{q}}/\hbar}+\hat{H}_{c}/\hbar-\chi\hat{c}^{\dagger}\hat{c}\hat{q}^{\dagger}\hat{q}$
(13)

noting that we have performed the rotating-wave approximation (RWA) and
discarded terms rotating at $\mathcal{O}(\omega_{q})$. Furthermore, the cross-
Kerr term remains unaffected since it is proportional to
$\hat{q}^{\dagger}\hat{q}$. Now, we diagonalize $\hat{H}_{q}$ and re-express
it in the resulting eigenbasis:

$\hat{H}_{\textrm{q}}/\hbar=\sum_{i}\epsilon_{i}\ket{i}\bra{i}$ (14)

where we label $i\in\\{+,-,\tilde{f},...\\}$ in correspondence with the fact
that we will be working in a regime where the two lowest driven eigenstates
strongly resemble those of a standard qubit that is driven on resonance, but
now incorporate a weak dressing with higher levels of the transmon. We
identify the Rabi frequency to be the energy difference between the lowest two
eigenstates $\epsilon_{+}-\epsilon_{-}=\Omega_{R}$ and define an effective
anharmonicity as
$\epsilon_{-}-\epsilon_{\tilde{f}}=\Omega_{R}+\tilde{\alpha}$. At this stage,
we turn to numerics and construct a unitary basis transformation between the
undriven and driven transmon eigenstates for a finite truncation of the
transmon Hilbert space. We then re-express the cross-Kerr interaction in the
driven basis, giving us:

$\hat{H}/\hbar=\sum_{i}\epsilon_{i}\ket{i}\bra{i}+\hat{H}_{c}/\hbar-\chi\hat{c}^{\dagger}\hat{c}\sum_{jk}u_{jk}\ket{j}\bra{k}$
(15)

We can further simplify this by going into the frame of the driven transmon
$\hat{U}=e^{i\hat{H}_{\textrm{Rabi}}t/\hbar}$ which performs the
transformations $\ket{j}\rightarrow e^{i\epsilon_{j}t/\hbar}\ket{j}$,
resulting in:

$\hat{H}/\hbar=\hat{H}_{c}/\hbar-\chi\hat{c}^{\dagger}\hat{c}\underbrace{\sum_{jk}u_{jk}e^{i(\epsilon_{j}-\epsilon_{k})t/\hbar}\ket{j}\bra{k}}_{\hat{q}^{\dagger}\hat{q}}$
(16)

where $u_{jk}=u^{*}_{kj}$. We consider the terms associated with the lowest
three levels explicitly:

$\displaystyle\begin{split}\hat{q}^{\dagger}\hat{q}=\,\,&u_{++}\ket{+}\bra{+}+u_{--}\ket{-}\bra{-}+u_{\tilde{f}\tilde{f}}\tilde{\ket{f}}\tilde{\bra{f}}+u_{+-}e^{i\Omega_{R}t}\ket{+}\bra{-}+u_{-+}e^{-i\Omega_{R}t}\ket{-}\bra{+}\\\
&+u_{-\tilde{f}}e^{i(\Omega_{R}+\tilde{\alpha})t}\ket{-}\tilde{\bra{f}}+u_{\tilde{f}-}e^{-i(\Omega_{R}+\tilde{\alpha})t}\tilde{\ket{f}}\bra{-}+u_{+\tilde{f}}e^{i(2\Omega_{R}+\tilde{\alpha})t}\ket{+}\tilde{\bra{f}}+u_{\tilde{f}+}e^{-i(2\Omega_{R}+\tilde{\alpha})t}\tilde{\ket{f}}\bra{+}\end{split}$
(17)

At this stage, we pause and turn to simplify $\hat{H}_{c}$. First, we choose
to parameterize the two drive frequencies
$\omega_{1/2}=\omega_{c}-\Delta_{c}\mp\Omega_{R}$. By going into the rotating
frame at the average of the drive frequencies
$\hat{U}=e^{i(\omega_{c}-\Delta_{c})t\hat{c}^{\dagger}\hat{c}}$, we arrive at:

$\hat{H}/\hbar=\hat{H}_{\textrm{int}}/\hbar+\Delta_{c}\hat{c}^{\dagger}\hat{c}-\varepsilon_{1}(\hat{c}e^{-i\Omega_{R}t+i\varphi_{1}}+\hat{c}^{\dagger}e^{i\Omega_{R}t-i\varphi_{1}})-\varepsilon_{2}(\hat{c}e^{i\Omega_{R}t+i\varphi_{2}}+\hat{c}^{\dagger}e^{-i\Omega_{R}t-i\varphi_{2}})$
(18)

Finally, we assume that the drive strengths are equal
$\varepsilon_{2}=-\varepsilon_{1}=\varepsilon$ and parameterize the drive
phases as their sum and differential components
$\varphi_{\Sigma}=(\varphi_{1}+\varphi_{2})/2$ and
$\varphi_{\delta}=(\varphi_{1}-\varphi_{2})/2$. This allows us to further
simplify our Hamiltonian to:

$\hat{H}/\hbar=\hat{H}_{\textrm{int}}/\hbar+\Delta_{c}\hat{c}^{\dagger}\hat{c}-2i\varepsilon\textrm{sin}(\Omega_{R}t-\varphi_{\delta})(\hat{c}e^{i\varphi_{\Sigma}}-\hat{c}^{\dagger}e^{-i\varphi_{\Sigma}})$
(19)

We can observe here that the sum phase of the two sidebands contributes simply
as a static rotation of $\hat{c}$, therefore we can always align to this frame
by experimentally adjusting this phase. Hence, we will set
$\varphi_{\Sigma}=0$ here on out to simplify our expressions. At this stage,
we aim to eliminate this time-dependent drive term by performing a
displacement transformation
$\hat{U}=e^{\xi(t)\hat{c}^{\dagger}-\xi^{*}(t)\hat{c}}$. This is achieved by
choosing
$\xi(t)=\frac{2\varepsilon}{\Omega_{R}}\textrm{cos}(\Omega_{R}t+\varphi_{\delta})=\xi_{0}(e^{i(\Omega_{R}t+\varphi_{\delta})}+e^{-i(\Omega_{R}t+\varphi_{\delta})})$
where $\xi_{0}=\frac{\varepsilon}{\Omega_{R}}$, which also transforms
$\hat{c}\rightarrow\hat{c}$ \+ $\xi(t)$. This gives:

$\displaystyle\begin{split}\hat{H}/\hbar&=\Delta_{c}(\hat{c}^{\dagger}+\xi^{*}(t))(\hat{c}+\xi(t))-\chi(\hat{c}^{\dagger}+\xi^{*}(t))(\hat{c}+\xi(t))\sum_{jk}u_{jk}e^{i(\epsilon_{j}-\epsilon_{k})t/\hbar}\ket{j}\bra{k}\\\
&=\Delta_{c}(\hat{c}^{\dagger}\hat{c}+\xi(t)(\hat{c}+\hat{c}^{\dagger})+\xi^{2}_{0})-\chi(\hat{c}^{\dagger}\hat{c}+\xi_{0}(e^{i(\Omega_{R}t+\varphi_{\delta})}+e^{-i(\Omega_{R}t+\varphi_{\delta})})(\hat{c}+\hat{c}^{\dagger})+\xi^{2}_{0})\sum_{jk}u_{jk}e^{i(\epsilon_{j}-\epsilon_{k})t/\hbar}\ket{j}\bra{k}\end{split}$
(20)

By substituting the expansion for $\hat{q}^{\dagger}\hat{q}$, discarding terms
that rotate at $\Omega_{R}$ and higher, and neglecting constant offsets, we
are left with an effective static interaction Hamiltonian:

$\hat{H}/\hbar=\Delta_{c}\hat{c}^{\dagger}\hat{c}-\chi\xi_{0}u_{+-}(e^{-i\varphi_{\delta}}\ket{+}\bra{-}+e^{i\varphi_{\delta}}\ket{-}\bra{+})(\hat{c}+\hat{c}^{\dagger})-\chi\hat{c}^{\dagger}\hat{c}(u_{++}\ket{+}\bra{+}+u_{--}\ket{-}\bra{-}+u_{\tilde{f}\tilde{f}}\tilde{\ket{f}}\tilde{\bra{f}})$
(21)

Importantly, this approximation requires larger Rabi frequencies as we drive
harder to induce larger desired interaction strengths (see subsection A. for a
more detailed analysis). Finally, by neglecting the final term (see subsection
B.), we arrive at the conditional displacement Hamiltonian between the
transmon and one cavity mode:

$\hat{H}/\hbar=\Delta_{c}\hat{c}^{\dagger}\hat{c}-g(\textrm{cos}(\varphi_{\delta})\hat{\sigma}_{x}+\textrm{sin}(\varphi_{\delta})\hat{\sigma}_{y})(\hat{c}+\hat{c}^{\dagger})$
(22)

where $g=\chi\xi_{0}u_{+-}\approx\frac{\chi\xi_{0}}{2}$ and we have defined
our Pauli operators such that
$\hat{\sigma}_{z}=\ket{+}\bra{+}-\ket{-}\bra{-}$. Here, we have formally
identified the qubit that we will use within the larger driven transmon
Hilbert space. As we can see, the coupling axis of the qubit fully depends on
the differential phase $\varphi_{\delta}$ of the cavity sidebands relative to
the qubit phase (which we have defined as zero here), which can be easily
adjusted experimentally without invoking additional Hamiltonian terms. Note
that throughout our derivations, we have assumed no accidental frequency
collisions that bring unintended Hamiltonian terms into resonance.

This scheme extends relatively straightforwardly to multiple cavity modes
dispersively coupled to the same transmon. Incorporating a pair of sidebands
on each additional cavity is sufficient to activate a conditional displacement
involving the driven transmon, as long as the resonance condition is
satisfied. It is also worth noting that we have assumed that each sideband
couples only to a single cavity mode - in practice, finite crosstalk
complicates the calibration procedure of activating all sidebands together. In
our experiment, we enact the aforementioned Hamiltonian Eq. (22) for two
cavity modes $\hat{c}\in\\{\hat{a},\hat{b}\\}$ coupled to orthogonal axes of
the qubit.

###  II.1 Optimizing the static cross-Kerr

Here, we consider (Eq. 20) and ask the questions: How large do we need the
Rabi frequency to be in order to safely discard all the rotating terms? Does
the answer to this question inform any design choices with regard to our
static Hamiltonian? To answer this question, we consider all of the different
terms that rotate at $\Omega_{R}$, neglecting any phases:

$\displaystyle\big{(}\Delta_{c}\xi_{0}(\hat{c}+\hat{c}^{\dagger})-\chi
u_{+-}\hat{c}^{\dagger}\hat{c}\hat{\sigma}_{\mp}-\chi\xi^{2}_{0}u_{+-}\hat{\sigma}_{\mp}\big{)}e^{\pm
i\Omega_{R}t}$

Note that there are also terms that rotate at
$\tilde{\alpha},2\Omega_{R},\Omega_{R}+\tilde{\alpha},2\Omega_{R}+\tilde{\alpha}$,
and $3\Omega_{R}+\tilde{\alpha}$, but since the prefactors will all be of the
same order, we consider the smallest rotating frequency for the most stringent
condition. Importantly, we also require that $\tilde{\alpha}>\Omega_{R}$,
otherwise other terms involving $\tilde{\ket{f}}$ will be activated and we can
no longer restrict ourselves to a qubit subspace. This sets a limit on how
large of a Rabi frequency can be used for a fixed anharmonicity $\alpha_{q}$.
The above terms reveal that our conditions for the RWA are:

$\Omega_{R}\gg\\{\Delta_{a}|\xi_{0}\langle(\hat{c}+\hat{c}^{\dagger})\rangle|,\frac{\chi}{2}\langle\hat{c}^{\dagger}\hat{c}\rangle,\frac{\chi}{2}\xi^{2}_{0}\\}$
(23)

which notably depends on the state of the cavity. It is clear from this that
as the conditional displacement interaction strength
$g\approx\frac{\chi\xi_{0}}{2}$ increases, the approximation becomes less
valid. However, we can instead rewrite the condition for a fixed $g$:

$\Omega_{R}\gg\\{\Delta_{a}|\xi_{0}\langle(\hat{c}+\hat{c}^{\dagger})\rangle|,\frac{g}{\xi_{0}}\langle\hat{c}^{\dagger}\hat{c}\rangle,g\xi_{0}\\}$
(24)

which reveals that there is indeed an optimal value for $\xi_{0}$ given a
fixed $g$. For considering photon numbers
$\langle\hat{c}^{\dagger}\hat{c}\rangle\sim\mathcal{O}(1)$ and
$g\approx\Delta_{c}$, we best satisfy all these conditions by choosing
$\xi_{0}\approx 1$. This, in turn for a fixed $g$, suggests that we should
roughly target a static cross-Kerr strength of $\chi\approx 2g$.

###  II.2 Choice of static detuning of the Rabi drive

The final term in Eq. (21) represents an effective cross-Kerr interaction
between cavity photons and the driven transmon eigenstates. For a true two-
level system driven on resonance, which is a good approximation for transmons
in the regime that the Rabi frequency is much weaker than the anharmonicity
$\Omega_{R}\ll\alpha_{q}$, one finds that $u_{++}=u_{--}=1/2$ which results in
a static frequency shift of the cavity and hence a nulled cross-Kerr. As the
Rabi frequency $\Omega_{R}$ approaches the anharmonicity $\alpha_{q}$,
however, $u_{++}\neq u_{--}$ for a drive that is on resonance owing to the
hybridization of the driven eigenstates with higher energy levels of the
transmon. This results in a residual cross-Kerr which can be interpreted as a
slight shift in the Rabi frequency due to the presence of photons in the
cavity. This is problematic as it both changes the resonance condition of the
interaction and biases our measurement scheme as a function of the cavity
photon distribution.

Figure 1: Optimizing the static detuning. By numerically diagonalizing
$\hat{H}_{q}/\hbar$ for various values of $\Delta_{R}$, we can plot the
dimensionless factor that contributes to the residual cross-Kerr for
$\Omega_{R}/2\pi=80$ MHz and $\alpha_{q}/2\pi=244$ MHz. We find an optimal
value of $\Delta_{R}/2\pi\approx$ 7 MHz.

By adding an additional static detuning knob $\Delta_{R}$ on the Rabi drive,
we can determine an optimal working configuration that nulls this effective
cross-Kerr. We show this optimization in FIG. 1. The presence of this static
detuning thus dictates that we perform an adiabatic preparation of our driven
qubit eigenstates. This has the further benefit of eliminating leakage events
associated with large Rabi frequencies and finite transmon anharmonicity, but
has an additional challenge which we address in Section V. D.

##  III Measuring the driven qubit

In order to measure the qubit along its original driven basis, i.e.
$\hat{\sigma}_{x}$, we need to know the total phase accumulated on the excited
state in the lab frame in order to perform the final $\pi/2$ rotation (in
order to decode onto our measurement basis) around the correct axis. As
depicted in Fig. 2b of the main text, this phase should be equal to the Rabi
frequency multiplied by the total free evolution duration plus an additional
offset phase that is fixed for a given set of ramp times and drive amplitudes.
Since our decoding rotation is performed with a local oscillator that is
locked to the $\ket{g}\leftrightarrow\ket{e}$ transition frequency of the
undriven transmon, the proper phase to perform this rotation is
$(\Omega_{R}-\Delta_{R})\tau+\varphi_{0}$, which notably depends on time.

Figure 2: State tomography and decoherence of transverse Bloch vectors. a. A
Ramsey-style pulse sequence is used to calibrate the decode rate. The final
$\pi/2$ rotation has a phase that depends on both the delay time $\tau$ and
the programmed decode rate $\Omega_{d}$. b. The 2D plot reveals an optimal
decode rate which is inferred to be equal to the true Rabi frequency offset by
the static detuning $\Omega_{R}-\Delta_{R}$. c. The decode rate can then be
fixed to the optimal value to extract a driven transverse relaxation time
$T_{2\rho}\approx$ 27 $\mu$s via fitting to an exponential decay function
(solid line).

FIG. 2 shows a calibration experiment where we scan the decode rate and the
delay time in order to reveal the Rabi frequency (offset by the static
detuning), which can then be fit to extract a driven transverse decoherence
time $T_{2\rho}$ Gustavsson2012 .

##  IV Experimental Setup and Device Parameters

A table of the static system parameters is given in TABLE I. A schematic of
the wiring diagram for this experiment is depicted in FIG. 3. We comment below
(subsection A.) on notable features not explicitly shown in the wiring
diagram.

###  IV.1 Wiring Diagram

An FPGA-based quantum controller synchronizes multiple modules that contain
DACs and ADCs for generating the pulses (I & Q control) and digitizing the
readout signals, respectively. RF switches are only open while pulses are
played on a corresponding mode. The control line for Alice is split between
the resonant drive (left) and sideband pumps (right), which are never played
simultaneously (see Fig. 2b of the main text), thus a closed loop where both
RF switches are open together is never formed. The bandstop filter on the pump
line is centered on Alice’s resonance frequency, suppressing pump-induced
noise that may lead to heating of the cavity mode. DC blocks are placed around
each active component, as well as on each line at the boundary of the cryostat
(i.e. the line separating 300 K and 4 K in the schematic). All components in
the 20 mK region are thermally anchored to the mixing chamber plate via OFHC
copper braids. We use Josephson Parametric Converters (JPCs) as quantum
limited amplifiers - only coupling to the signal port is shown.

Figure 3: Experimental wiring diagram. System quantity | Parameter | Value
---|---|---
Transmon frequency | $\omega_{q}/2\pi$ | 4850 MHz
Transmon anharmonicity | $\alpha_{q}/2\pi$ | 244 MHz
Transmon relaxation | $T^{q}_{1}$ | 80 $\mu$s
Transmon decoherence | $T^{q}_{2}$ | 7 $\mu$s
Alice frequency | $\omega_{a}/2\pi$ | 5436 MHz
Alice linewidth | $\kappa_{a}/2\pi$ | 0.23 kHz
Bob frequency | $\omega_{b}/2\pi$ | 6506 MHz
Bob linewidth | $\kappa_{b}/2\pi$ | 320 kHz
Alice - Transmon coupling | $\chi_{aq}/2\pi$ | 295 kHz
Bob - Transmon coupling | $\chi_{bq}/2\pi$ | 210 kHz
Ancilla frequency | $\omega_{qa}/2\pi$ | 4509 MHz
Ancilla relaxation | $T^{qa}_{1}$ | 60 $\mu$s
Ancilla decoherence | $T^{qa}_{2}$ | 10 $\mu$s
Alice - Ancilla coupling | $\chi_{a,qa}/2\pi$ | 845 kHz
Table 1: List of system parameters.

###  IV.2 Temperature stabilization

As the derivations of Section II. suggest, the resonance condition for
enacting our reaction Hamiltonian (Eq. (3) of the main text) relies on
matching the sideband detunings to the Rabi frequency. The Rabi frequency
depends linearly on the amplitude of the Rabi drive (roughly speaking, using
the two-level approximation that
$\Omega_{R}=\sqrt{\epsilon^{2}_{R}+\Delta^{2}_{R}}$ and we operate in a regime
where $\epsilon_{R}\gg\Delta_{R}$, resulting in
$\Omega_{R}\approx\epsilon_{R}(1+\frac{\Delta^{2}_{R}}{2\epsilon^{2}_{R}})$),
and thus is susceptible to amplitude fluctuations such as those caused by
variations in the gain at any stage of our microwave control chain. A dominant
source of these variations is due to ambient temperature fluctuations in the
lab. To this end, we suppress these fluctuations by anchoring the components
(particularly a Marki LXP IQ-mixer and MITEQ low noise amplifier, see
turquoise box in FIG. 3) to a Thorlabs optical breadboard and placing the
breadboard in a cardboard box. We then actively stabilize the temperature of
the air in the box via an op-amp based PID feedback controller that heats an
Ohmite ceramic resistor ($R$ = 2.5 $\Omega$) based on a differential
measurement of the temperature using a 100 k$\Omega$ thermistor referenced to
a set point. FIG. 4 shows the typical performance of our stabilization and
correlates the temperature variations with the amplitude variations as
measured via the Rabi frequency. Over the course of 24 hours, we achieve an
absolute temperature stability within 50 mK and a relative amplitude stability
of 100 kHz / 72.9 MHz $\approx 10^{-3}$, suggesting that we have a relative
amplitude sensitivity of 1% per 500 mK. We note that the timscale for a
typical calibration and measurement of a dataset presented in this paper is
roughly a few hours, meaning we can operate in a window where the relative
amplitude stability can be much better than $10^{-3}$.

Figure 4: Tracking system stability over time. (Top panel) Simultaneous
measurements of the Rabi frequency and the ambient temperature where the
active microwave components are held reveals correlations between temperature
drifts and amplitude drifts. (Middle panel) Extracting driven coherences
$T_{2\rho}$ suggests a stable amplitude noise spectrum within an acquisition
time $\tau_{\textrm{acq}}$ = 3 minutes (bottom left panel: a typical time-
domain Ramsey trace), with a single instance where the amplitude drift was
large (bottom right panel), as confirmed via looking at the raw data binned 10
shots at a time over $\tau_{\textrm{acq}}$.

##  V Calibration Procedures

In this section, we provide an overview on the calibration procedures in our
experiment. We calibrate the interaction between each individual cavity
sideband and our driven qubit before proceeding to combine sidebands to enact
conditional displacements. As more drives are added, care is taken to re-
optimize the resonance conditions given that the presence of each additional
drive will result in a Stark shift of each mode, particularly the transmon
whose Rabi frequency we are very sensitive to.

###  V.1 Individual sideband interactions

The sideband interactions have the form:

$\displaystyle\hat{H}_{\textrm{red}}/\hbar$
$\displaystyle=g\hat{c}\hat{\sigma}^{-}_{x}+g^{*}\hat{c}^{\dagger}\hat{\sigma}^{+}_{x}$
(25) $\displaystyle\hat{H}_{\textrm{blue}}/\hbar$
$\displaystyle=g\hat{c}\hat{\sigma}^{+}_{x}+g^{*}\hat{c}^{\dagger}\hat{\sigma}^{-}_{x}$
(26)

for a general cavity annihilation operator $\hat{c}$ and
$\hat{\sigma}^{\pm}_{x}=\ket{\mp}\bra{\pm}$ are the raising and lowering
operators of the two driven qubit eigenstates that adiabatically connect to
the ground and first excited state of the transmon. For our system, we have
$\hat{c}\in\\{\hat{a},\hat{b}\\}$. The qualitative behavior of each individual
sideband interacting with the driven qubit will be different given that we are
operating in the regime where $g>\kappa_{a}$ and $g\lesssim\kappa_{b}$. The
former will result in either creating and annihilating two excitations
simultaneously $\ket{+,0}\leftrightarrow\ket{-,1}$ (red sideband) or a
coherent exchange between an excitation in the driven qubit and a photon in
the cavity $\ket{-,0}\leftrightarrow\ket{+,1}$ (blue sideband) Lu2017 . The
latter will stabilize the qubit in either the excited state $\ket{-}$ (red
sideband) or the ground state $\ket{+}$ (blue sideband) Murch2012 .

In order to calibrate the interaction strengths, we will operate in the
restricted two-dimensional subspace of the joint Hilbert space of the cavity
and qubit as described above. This allows us to simplify our analysis and
replace the qubit raising and lowering operators $\hat{\sigma}^{\pm}_{x}$ with
general bosonic creation and annihilation operators $\hat{d}^{(\dagger)}$. We
then follow Refs. Pfaff2017 ; Grimm2020  and capture the full range of
dynamics by solving the equations of motion for $\hat{d}$ under
$\hat{H}_{\textrm{blue}}$ and incorporating a cavity damping rate $\kappa$. We
also include a static detuning term $\delta\hat{c}^{\dagger}\hat{c}$ to
capture the effect of sweeping the pump frequency that enables the
interaction. The resulting field has the form:

$\hat{d}(t)=\frac{\hat{d}(0)}{\Omega}e^{-\frac{\kappa_{\textrm{eff}}t}{4}}\bigg{(}\Omega\textrm{cosh}\big{(}\frac{\Omega
t}{4}\big{)}+\kappa_{\textrm{eff}}\textrm{sinh}\big{(}\frac{\Omega
t}{4}\big{)}\bigg{)}$ (27)

where $\Omega=\sqrt{\kappa^{2}_{\textrm{eff}}-(4g)^{2}}$ and
$\kappa_{\textrm{eff}}=\kappa+2i\delta$.

For each interaction, we prepare our system in either $\ket{-,0}$ (blue
sidebands) or $\ket{+,0}$ (red sidebands) and scan the frequency of the cavity
sideband and the delay time for a given pump amplitude. By measuring the qubit
population, we extract $\langle\hat{d}^{\dagger}(t)\hat{d}(t)\rangle$ and can
fit the resulting data using Eq. (27). For the exchange interaction under the
blue sideband, we have
$\langle\hat{d}^{\dagger}(0)\hat{d}(0)\rangle_{\textrm{blue}}=1$ for the qubit
initially in its excited state, giving us an expression for
$\langle\hat{d}^{\dagger}(t)\hat{d}(t)\rangle_{\textrm{blue}}$. For the red
sideband, the features are qualitatively identical, with the exception that
the qubit starts out in the ground state, giving us
$\langle\hat{d}^{\dagger}(t)\hat{d}(t)\rangle_{\textrm{red}}=1-\langle\hat{d}^{\dagger}(t)\hat{d}(t)\rangle_{\textrm{blue}}.$
Notably, this assumes that the effective interaction strength $g$ is
independent of the pump detuning $\delta$, which is strictly not true but
should be a very good approximation in our regime given that the scale of the
chevron features, set by $g$, is much smaller than the absolute detuning from
the cavity resonance $\sim\Omega_{R}$. We allow for an overall amplitude,
global offset, and time offset in our fit, leaving $g$ and $\kappa$ to be the
only free parameters. In the case of Bob, we first perform this fit for a
range of interaction strengths $g\approx\kappa_{b}$ and extract a decay rate
$\kappa_{b}/2\pi\approx 320$ kHz. For the remainder of the calibrations where
$g<\kappa_{b}$, we fix this quantity and let the interaction strength $g$ be
the only free parameter to be fitted. The full calibration for different pump
amplitudes is shown in FIG. 5.

Figure 5: Calibrating individual exchange interactions a. Interaction
strengths $g$ for individual sideband interactions between the driven qubit
and Alice (top panel) or Bob (bottom panel). $g$ is extracted via a fit using
Eq. (27) for the red sidebands (red squares) and blue sidebands (blue
squares). Error bars are smaller than the markers. $\xi$ is calibrated via
independent measurements of $\chi_{a/b}$ and the qubit’s Stark shift as a
function of pump amplitude. Solid lines represent the expected interaction
strength in the ideal case $g=\frac{\chi_{a/b}\xi}{2}$. Deviations of the data
from a linear relation for Alice coincides with regions where the Stark shift
becomes nonlinear in pump power. Slight disagreement between the data and
theory expectation for Bob may stem from an inaccurate estimate of $\chi_{b}$
and from a sensitivity of the extracted fit values to a true value of
$\kappa_{b}$. b. Raw calibration data of the qubit population for a fixed pump
amplitude (dashed circles in a.) for Alice (top) and Bob (bottom). The red
sideband (left) either coherently creates and annihilates correlated photons
in the driven qubit and the resonator ($g>\kappa_{a}$, Alice) or stabilizes
the driven qubit in the first excited state ($g\lesssim\kappa_{b}$, Bob). The
opposite is true for the blue sideband (right), which either coherently
exchanges a photon between the driven qubit and resonator (Alice) or
stabilizes the driven qubit in the ground state (Bob). Simultaneous
measurement of the photon population in Alice via its ancilla qubit reveals
the aforementioned correlations. The sideband detunings are referenced to the
cavity frequencies, showing that we are operating at Rabi frequencies
$\Omega_{R}/2\pi\approx 80$ MHz.

###  V.2 Calibrating conditional displacements

The combination of simultaneous red and blue cavity sidebands enacts a
conditional displacement interaction. This requires that two conditions are
fulfilled: 1) the interaction strengths of each individual sideband
interaction be equal and 2) the frequency difference between the two sidebands
equals twice the Rabi frequency. If these two conditions are met, then we can
model the interaction using Eq. (20).

In practice, the presence of each sideband will Stark shift both the transmon
and cavity modes. Thus, to capture the dominant effect of all of these Stark
shifts (which influences the resonance condition), we perform individual
sideband calibrations with the opposite sideband on but detuned by an amount
larger than the interaction strength we are using (i.e. by an additional 2 MHz
in our experiments). We scan the pump amplitudes and match the individual
sideband strengths before bringing both sidebands into resonance. This relies
on the assumption that over a variation of $\sim$2 MHz, the relative change in
the cavity Stark shift, which influences the dimensionless pump strength that
determines the interaction strength, is negligible. Finally, we fine tune the
difference frequency of the two sidebands while keeping the average value
fixed (which fixes $\Delta_{c}$ in Eq. (20)) in order to account for any
change in the Rabi frequency which we are very sensitive to. This sensitivity
is revealed by measuring the transverse relaxation time $T_{2\rho}$, and
choosing a calibration point where this value is maximized, suggesting that
the resonance conditions are fulfilled as best as possible. This calibration
procedure gives us the data shown in Fig. 2 of the main text.

For Alice, the dynamics of an initial vacuum state $\ket{\alpha_{0}=0}_{a}$
evolving under a conditional displacement interaction will be a circular
trajectory in phase space around the location of the ground state
$\alpha_{g}=g_{x}/\Delta_{a}$. By measuring the population in $\ket{n=0}_{a}$,
we are effectively measuring the overlap of a coherent state with itself as it
oscillates in time. This justifies the use of a simple model, where the state
autocorrelation function is $\langle\beta e^{-i\Delta
t}|\beta\rangle=e^{|\beta|^{2}(e^{-i\Delta_{a}t}-1)}$. The corresponding
probability is:

$P_{0}=|\langle\beta e^{-i\Delta
t}|\beta\rangle|^{2}=e^{2|\beta|^{2}(\textrm{cos}[\Delta_{a}t]-1)}$ (28)

To make the connection with our model, we choose
$\beta=\alpha_{g}=g_{x}/\Delta_{a}$.

###  V.3 Calibrating cavity displacements along the interaction axis

As described in Section II., the cavity phase of the conditional displacement
(i.e. the phase which defines the position operator
$\hat{x}\propto(\hat{c}e^{i\varphi_{\Sigma}}+\hat{c}^{\dagger}e^{-i\varphi_{\Sigma}})$)
is determined by the sum phase of the red and blue sidebands. Given that we
are turning on the conditional displacement interaction suddenly
($\tau^{\textrm{sb}}_{\textrm{ramp}}\ll 1/g$), the phase of our initial
displacement operation $\hat{\mathcal{D}}(\alpha_{0})$ on Alice will determine
the location in the driven phase space where the wave-packet begins.
Displacements whose phase is aligned to the conditional displacement cavity
phase will prepare wave-packets along the position axis, whereas care needs to
be taken to prepare wave-packets with various momentums that are located at
one of the two ground state positions. Since we are interested in modeling
scenarios where a wave-packet arrives on a potential energy surface via
optical excitation in a Franck-Condon region, we prepare coherent states with
no initial momentum along the reaction coordinate.

In order to calibrate the displacement phase, we begin with a calibrated
conditional displacement where we have extracted $g_{x}$ and $\Delta_{x}$,
which gives us a value for the ground state amplitude
$\alpha_{g}=g_{x}/\Delta_{x}$ (Fig. 2c of the main text). Note that this does
not rely on any displacement phase since we are starting off in a vacuum
state. Next, we scan the phase of an initial displacement of $2\alpha_{g}$ and
turn on the conditional displacement interaction for various delay times. The
optimal phase will be the one where we recover revivals that are half a period
out of phase from those in Fig. 2c of the main text. This can be interpreted
as follows. A vacuum state in the lab frame $\ket{0}_{\textrm{lab}}$ looks
like a displaced state $\ket{-\alpha_{g}}_{\textrm{disp}}$ with respect to the
displaced ground state
$\ket{\alpha_{g}}_{\textrm{lab}}=\ket{0}_{\textrm{disp}}$, and thus will
oscillate around the ground state, reaching $\ket{\alpha_{g}}_{\textrm{disp}}$
after half a period. By determining the phase that enables us to prepare
$\ket{\alpha_{g}}_{\textrm{disp}}$ (which is
$\ket{2\alpha_{g}}_{\textrm{lab}}$ in the lab frame and will return to the
vacuum state $\ket{-\alpha_{g}}_{\textrm{disp}}$ after half a period) we can
prepare any state along the position axis, including the displaced ground
state.

Figure 6: Echoing away residual entanglement. Simultaneous measurement of
Alice’s vacuum projector via a transmon ancilla (left) and the driven qubit
(right) while a conditional displacement interaction is active. The pulse
sequence used for this calibration is shown in Fig. 2b of the main text
without (top) and with (bottom) the $\pi$ rotation and echo delay
$\tau_{\textrm{ED}}$ and only sidebands on Alice. There is no initial
displacement amplitude ($\alpha_{0}=0$). Time dynamics of Alice’s vacuum
projector reveals coherent revivals as expected (as in Fig. 2c of the main
text), independent of the decode rate on the driven qubit. Scanning the decode
rate without the echo sequence results in distortions to the driven Ramsey
data that are correlated with the cavity photon distribution. Implementing the
echo sequence eliminates this effect, suggesting that the systems remain
unentangled at all delay times. Here, we use an optimized value of the delay
time $\tau_{\textrm{ED}}=144$ ns.

###  V.4 Echoing the residual entanglement during ramping of the Rabi drive

The adiabatic preparation of our driven qubit eigenstates as motivated in
Section II. B. has two benefits: for a fixed transmon anharmonicity, we can 1)
use larger Rabi frequencies while cancelling the residual cross-Kerr and 2)
avoid leakage events to higher transmon levels (up to natural heating rates of
the dressed eigenstates) associated with a resonant drive. The primary
consequence of this approach is the undesired interaction between the qubit
and cavity photons during the ramp time of the Rabi drive.

Qualitatively, this undesired interaction stems from ramping between the
static interaction Hamiltonian
$\hat{H}_{\textrm{int}}/\hbar=-\chi\hat{c}^{\dagger}\hat{c}\hat{q}^{\dagger}\hat{q}$
and the driven Hamiltonian where the cross-Kerr interaction is nulled. Thus,
if there are photons in the cavity during either the ramp on or off of the
Rabi drive, they will entangle with the superposition states of the qubit that
we are manipulating. We can avoid this effect during the ramp on of the Rabi
drive by performing our displacement operation after the Rabi drive is fully
ramped on, i.e. during the time when the cross-Kerr interaction is nulled
(hence the relative position of the displacement operation in Fig. 2b of the
main text). For addressing the entanglement during the ramp off of the Rabi
drive, we implement a simple and short echo sequence of the qubit to un-do the
interaction. This works because the entanglement is fully determined by $\chi$
and the ramp time $\tau^{R}_{\textrm{ramp}}$ and not the cavity photon
distribution. This is important as we do not want a scheme which depends on
the cavity state that we are manipulating. FIG. 6 shows how implementing this
protocol eliminates spurious features that arise from this entanglement when
performing a decode calibration experiment (FIG. 2b) when a conditional
displacement interaction with Alice is active. In practice, this calibration
is only done with respect to photons in Alice. Given that we operate Bob in a
regime where $g_{y}<\kappa_{b}$, the photon distribution in Bob remains
relatively small and thus any residual entanglement effects are negligible.

##  VI Extended data & post-selecting on leakage events

In this section we present additional data (FIG. 7) that supports what is
shown in the main text, specifically focusing on leakage statistics.

Figure 7: Extended data and leakage statistics. a,b,c. Probability of leaking out of the qubit subspace to $\ket{f}$ for the data presented in Figs. {2b, 2c, 3b} of main text, respectively. d. The same dataset as appears in Fig. 4b of the main text (left), but including two additional initial states with $\ket{\alpha_{0}}_{a}=\ket{\alpha_{g}/2}_{a}$ (blue) and $\ket{3\alpha_{g}/2}_{a}$ (red), with the associated leakage probabilities (right). All of the leakage probabilities are equivalent to the percentage of data that is post-selected away for the rest of the data in the paper, which is always less than 5%. Dataset |  Drives (+ Rabi) | Figure |  $\Omega_{R}[\Delta_{c}]/2\pi$ (MHz) |  $T_{2\rho}$ [$\alpha_{0}$] ($\mu$s) | Leakage
---|---|---|---|---|---
Coherent state revivals | Alice sidebands | 2c | 82.532 [0.457] | not measured | SM Fig. 7a
82.532 [0.355]
82.492 [0.246]
Measurement-induced dephasing | Bob sidebands | 2d | 79.915 [0] | 36.6 | SM Fig. 7b
79.9175 [0.4] | 40.5
79.9225 [0.8] | 36.3
Coherent revivals | Alice & Bob sidebands | 3b | 80.045 | see Fig. 3b | SM Fig. 7c
\+ state preparation
Conical intersection dynamics | Alice & Bob sidebands | 4b | 81.013 | 51.7 [0] | SM Fig. 7d
51.1 [$\alpha_{g}/2$]
57.8 [$\alpha_{g}$]
53.9 [$3\alpha_{g}/2$]
48.8 [$2\alpha_{g}$]
Table 2: Extended data information. Quantities in brackets correspond to
identifiers within a dataset. Different values of $\Delta_{c}$,
$c\in\\{a,b\\}$ require a fine-tuning calibration on the drive configuration
to match Rabi frequency shifts, but preparing different initial states does
not.

Our measurement of the transmon is able to distinguish between
$\\{\ket{g},\ket{e},\ket{f}\\}$ on a single-shot basis, which gives us
information on leakage events outside the qubit manifold that we use for our
experiments. In the absence of decoherence, our adiabatic preparation scheme
should ideally eliminate any leakage to the second excited state $\ket{f}$ and
higher (assuming there are no accidental resonances induced by the drives). In
practice, any relaxation or heating between undriven transmon levels will lead
to transitions between the driven eigenstates that have finite support across
multiple undriven basis states. In all of the data presented, we post-select
away outcomes where the transmon is measured to be in $\ket{f}$. We compile
the post-selection statistics in FIG. 7, and note that overall the leakage
probabilities are small.

##  VII Dissipation analysis and Zeno dynamics

In this section we describe the oscillator-induced dissipation of the qubit
resulting from the evolution under the master equation,

$\displaystyle\dot{\hat{\rho}}=-i[\hat{H}_{\textrm{Zeno}},\hat{\rho}]+\kappa_{b}\left(\hat{b}\hat{\rho}\hat{b}^{\dagger}-\frac{1}{2}\hat{b}^{\dagger}\hat{b}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{b}^{\dagger}\hat{b}\right)$
(29)

with

$\displaystyle\hat{H}_{\textrm{Zeno}}=E(x)\hat{\sigma}_{x}+\Delta_{y}\hat{b}^{\dagger}\hat{b}+g\hat{\sigma}_{z}(\hat{b}^{\dagger}+\hat{b})$
(30)

from Eq. (5) of the main text. The dynamics of the system are effectively the
same as Zeno dynamics of a driven qubit generated by the competition between
$E(x)\hat{\sigma}_{x}$, which tries to lock the state of the qubit to an
eigenstate of $\hat{\sigma}_{x}$, and
$g\hat{\sigma}_{z}(\hat{b}^{\dagger}+\hat{b})$, which in combination with
photon loss in the oscillator tries to project the state of the qubit onto an
eigenstate of $\hat{\sigma}_{z}$. Note that here we choose the dissipation to
be along $\hat{\sigma}_{z}$ as opposed to $\hat{\sigma}_{y}$. The resulting
dissipative dynamics is crucially dependent on the parameters
$E(x),\Delta_{y},g,\kappa_{b}$. While an exact analytic expression for the
time-dependence of the qubit density matrix cannot be obtained in general,
here we consider a few extreme parameter regimes which can be easily analyzed
and provide a window to the vast range of qubit dynamics possible.

This simplest case is when $E(x)=0$. In this case, it is possible to use
Fokker-Planck equations to calculate an effective dephasing rate. However,
when $E(x)\neq 0$, this approach will fail to give an analytical expression
for the dephasing rate. In this case, we will invoke additional constraints
which will allow us to apply the Born-Markov approximation and derive an
effective master equation for the qubit.

###  VII.1 $E(x)\neq 0$

To get a simple intuitive understanding of the oscillator-induced dephasing
when $E(x)\neq 0$, we consider two cases: (a) $\Delta_{y}=0$, $\kappa_{b}\gg
g$ and (b) $|2\Delta_{y}-E(x)|\ll|2\Delta_{y}+E(x)|$,
$|2\Delta_{y}+E(x)|\gg\kappa_{b}\gg|2\Delta_{y}-E(x)|,g$.

Case (a): In this case, the oscillator mode can be eliminated using the
standard Born-Markov approximation and a simple effective master equation for
the qubit is obtained,

$\displaystyle\hat{\rho}_{q}=-i[\hat{H}_{q},\hat{\rho}_{q}]+\kappa_{q}\left(\hat{\sigma}_{z}\hat{\rho}_{q}\hat{\sigma}_{z}-\hat{\rho}_{q}\right)$
(31)

with $\hat{H}_{q}=E(x)\hat{\sigma}_{x}$,
$\kappa_{q}=g^{2}\kappa_{b}/(4E(x)^{2}+\kappa_{b}^{2}/4)$. In this case, the
state of the qubit will be an equal mixture of $\ket{+}$ and $\ket{-}$ or
$\ket{g}$ and $\ket{e}$.

Case (b): In this case, the standard Born-Markov approximation yields the
following master equation,

$\displaystyle\hat{\rho}_{q}$ $\displaystyle=-i[\hat{H}_{q},\hat{\rho}_{q}]$
$\displaystyle+\kappa_{q}\left(\ket{-}\bra{+}\hat{\rho}_{q}\ket{+}\bra{-}-\frac{1}{2}\ket{+}\bra{+}\hat{\rho}_{q}-\frac{1}{2}\hat{\rho}_{q}\ket{+}\bra{+}\right)$
(32)

with $\hat{H}_{q}=E(x)\hat{\sigma}_{x}$,
$\kappa_{q}=g^{2}\kappa_{b}/((2E(x)-\Delta_{y})^{2}+\kappa_{b}^{2}/4)$. Unlike
in case (a), here, the steady state of the qubit will be the pure state
$\ket{+}$. The difference between the master equations Eq. (31) and Eq. (32)
can be understood by observing that when
$|2\Delta_{y}+E(x)|\gg\kappa_{b}\gg|2\Delta_{y}-E(x)|$, then the transition
between $\ket{+}\otimes|0\rangle$ and $\ket{-}\otimes|1\rangle$ becomes more
likely than that between $\ket{+}\otimes|1\rangle$ to
$\ket{-}\otimes|0\rangle$ because the former, happening at frequency
$2\Delta_{y}-E(x)$ lies within the bandwidth $\kappa_{b}$. The first term in
the tensor product refers to the state of the qubit and the second term refers
to the vacuum $\ket{0}$ and single-photon Fock state $\ket{1}$ of the
oscillator mode. Consequently, if the photon is subsequently lost from the
oscillator, the qubit will be preferentially projected to the state $\ket{-}$.

Thus, the steady-state value of $\langle\hat{\sigma}_{x}\rangle$ is 0 for case
(a) and $-1$ for case (b). Clearly, we see that the qubit dynamics depends on
the parameters $\\{E(x),\Delta_{y},g,\kappa_{b}\\}$, and in certain cases can
have a steady state value of $\langle\hat{\sigma}_{x}\rangle$ between 0 and
$-1$.

We note here that by considering our full model (Eq. (1) of the main text)
with not only dissipation on the coupling mode $\kappa_{b}$ but also
dissipation on the tuning mode $\kappa_{a}$, we can expect a more complex
landscape of steady state dynamics which may qualitatively fall into the above
regime.

###  VII.2 $E(x)=0$

In this case, we follow the approach in Ref. Gambetta2008  and begin by
writing the density matrix of the qubit and oscillator as

$\displaystyle\hat{\rho}=\ket{g}\bra{g}\otimes\hat{\rho}_{gg}+\ket{e}\bra{e}\otimes\hat{\rho}_{ee}+\ket{g}\bra{e}\otimes\hat{\rho}_{ge}+\ket{e}\bra{g}\otimes\hat{\rho}_{eg}$
(33)

where $\ket{g,e}$ represent the qubit states and $\hat{\rho}_{ij}$ acts on the
oscillator Hilbert space conditioned on the qubit state. From Eq. (29), (30)
with $E(x)=0$ and Eq. (33) we get,

$\displaystyle\dot{\hat{\rho}}_{ge}$
$\displaystyle=-ig(\hat{b}^{\dagger}+\hat{b}){\hat{\rho}}_{ge}-ig{\hat{\rho}}_{ge}(\hat{b}^{\dagger}+\hat{b})-i\Delta_{y}(\hat{b}^{\dagger}\hat{b}{\hat{\rho}}_{ge}-{\hat{\rho}}_{ge}\hat{b}^{\dagger}\hat{b})$
$\displaystyle+\kappa_{b}\hat{b}{\hat{\rho}}_{ge}\hat{b}^{\dagger}-\frac{\kappa_{b}}{2}\hat{b}^{\dagger}\hat{b}{\hat{\rho}}_{ge}-\frac{\kappa_{b}}{2}{\hat{\rho}}_{ge}\hat{b}^{\dagger}\hat{b}$
(34)

Next, we use the positive-P representation ${\hat{\rho}}_{ge}=\int
P(\alpha,\alpha^{*},t)d\alpha d\alpha^{*}$ and write the effective Fokker-
Planck equation,

$\displaystyle\frac{\partial P}{\partial t}$
$\displaystyle=-2igP(\alpha+\alpha^{*})+i(g+\Delta_{y}\alpha)\frac{\partial
P}{\partial\alpha}$ $\displaystyle+i(g-\Delta_{y}\alpha^{*})\frac{\partial
P}{\partial\alpha^{*}}+\kappa_{b}P+\frac{\kappa_{b}\alpha}{2}\frac{\partial
P}{\partial\alpha}$
$\displaystyle+\frac{\kappa_{b}\alpha^{*}}{2}\frac{\partial
P}{\partial\alpha^{*}}$ (35)

Now we must solve the above equation with some given initial condition. In the
setup of interest we start with the oscillator mode in vacuum so that
$P(\alpha,\alpha^{*},0)=\delta^{2}(\alpha)=\lim_{\varepsilon\rightarrow
0}(1/\pi\varepsilon)\exp(-|\alpha|^{2}/\varepsilon)$ and the qubit in the
state $\ket{+}$. We can make the Gaussian ansatz,
$P(\alpha,\alpha^{*},t)=\exp(-a(t)+b(t)\alpha+c(t)\alpha^{*}-d(t)|\alpha|^{2})$
and write the equivalent differential equations for $a,b,c,d$ to get,

$\displaystyle-\dot{a}$ $\displaystyle=igb+igc+\kappa_{b}$
$\displaystyle-\dot{d}$ $\displaystyle=-\kappa_{b}d$ $\displaystyle\dot{b}$
$\displaystyle=-2ig+\left(i\Delta_{y}+\frac{\kappa_{b}}{2}\right)b-igd$
$\displaystyle\dot{c}$
$\displaystyle=-2ig+\left(-i\Delta_{y}+\frac{\kappa_{b}}{2}\right)b-igd$ (36)

These time-dependent equations can be easily solved with initial conditions
now written as $a(0)=\ln\pi\varepsilon,b(0)=0,c(0)=0,d(0)=1/\varepsilon$ but
the analytic expressions are considerably simplified when $\Delta_{y}=0$. Once
we get $a,b,c,d$ we are able to reconstruct $P$ and hence $\hat{\rho}_{ge}$.
The relevant quantity of interest is the time-dependence of the expectation
value of the operator $\ket{g}\bra{e}$,

$\displaystyle\langle\ket{g}\bra{e}(t)\rangle=Tr[\hat{\rho}_{ge}]=\int
P(\alpha,\alpha^{*},t)d^{2}\alpha$ (37)

For $\Delta_{y}=0$ and an initial condition where the qubit is in $\ket{+}$
and the oscillator in vacuum, we have

$\displaystyle\langle\ket{g}\bra{e}(t)\rangle=\frac{1}{2}\exp\left[-\frac{8g^{2}t}{\kappa_{b}}+\frac{16g^{2}}{\kappa_{b}^{2}}(1-e^{-\kappa_{b}t/2})\right]$
(38)

Next, we consider the rate equation for ${\hat{\rho}}_{gg}$,

$\displaystyle\dot{\hat{\rho}}_{gg}$
$\displaystyle=-ig(\hat{b}^{\dagger}+\hat{b}){\hat{\rho}}_{gg}+ig{\hat{\rho}}_{gg}(\hat{b}^{\dagger}+\hat{b})-i\Delta_{y}(\hat{b}^{\dagger}\hat{b}{\hat{\rho}}_{gg}-{\hat{\rho}}_{gg}\hat{b}^{\dagger}\hat{b})$
$\displaystyle+\kappa_{b}\hat{b}{\hat{\rho}}_{gg}\hat{b}^{\dagger}-\frac{\kappa_{b}}{2}\hat{b}^{\dagger}\hat{b}{\hat{\rho}}_{gg}-\frac{\kappa_{b}}{2}{\hat{\rho}}_{gg}\hat{b}^{\dagger}\hat{b}$
(39)

Like before, we use the positive-P representation ${\hat{\rho}}_{gg}=\int
P(\alpha,\alpha^{*},t)d\alpha d\alpha^{*}$, so that

$\displaystyle\frac{\partial P}{\partial t}$ $\displaystyle=-ig\frac{\partial
P}{\partial\alpha^{*}}+ig\frac{\partial
P}{\partial\alpha}-i\Delta_{y}\frac{\partial\alpha^{*}P}{\partial\alpha^{*}}+i\Delta_{y}\frac{\partial\alpha
P}{\partial\alpha}$ $\displaystyle+\frac{\kappa_{b}\alpha}{2}\frac{\partial
P}{\partial\alpha}+\frac{\kappa_{b}\alpha^{*}}{2}\frac{\partial
P}{\partial\alpha^{*}}$ (40)

We can again make the Gaussian ansatz
$P(\alpha,\alpha^{*},t)=\exp(-a(t)+b(t)\alpha+c(t)\alpha^{*}-d(t)|\alpha|^{2})$,
and solve the corresponding differential equations for $a,b,c,d$ to get
$P(\alpha,\alpha^{*},t)=\delta^{2}(\alpha-\alpha(t))$ where
$\alpha(t)=2ig(1-\exp(-\kappa_{b}t/2))/\kappa_{b}$ (for $\Delta_{y}=0$). Thus,
$\langle|g\rangle\langle g|(t)\rangle=1/2$. Similarly we can show that
$\langle|e\rangle\langle e|(t)\rangle=1/2$. Hence the diagonal terms of the
qubit density matrix do not decay with time — only the off-diagonal terms do.
Thus, in this case the qubit undergoes pure dephasing due to its coupling with
the oscillator mode and at steady state becomes an equal mixture of $\ket{g}$
and $\ket{e}$.

## References

  * (1) Conical Intersections: Electronic Structure, Dynamics & Spectroscopy, edited by W. Domcke, D. R. Yarkony, and H. Köppel (World Scientific Pub Co. Inc., Singapore, 2004).
  * (2) S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh, K.B. Whaley, and I. Siddiqi, Quantum Dynamics of Simultaneously Measured Non-Commuting Observables, Nature 538, 491 (2016).
  * (3) S. Gustavsson et al., Driven Dynamics and Rotary Echo of a Qubit Tunably Coupled to a Harmonic Oscillator, Phys. Rev. Lett. 108, 170503 (2012).
  * (4) Y. Lu, S. Chakram, N. Leung, N. Earnest, R. K. Naik, Z. Huang, P. Groszkowski, E. Kapit, J. Koch, and D. I. Schuster, Universal Stabilization of a Parametrically Coupled Qubit, Phys. Rev. Lett. 119, 150502 (2017).
  * (5) K. W. Murch, U. Vool, D. Zhou, S. J. Weber, S. M. Girvin, and I. Siddiqi, Cavity-Assisted Quantum Bath Engineering, Phys. Rev. Lett. 109, 183602 (2012).
  * (6) W. Pfaff, C. J. Axline, L. D. Burkhart, U. Vool, P. Reinhold, L. Frunzio, L. Jiang, M. H. Devoret, and R. J. Schoelkopf, Controlled Release of Multiphoton Quantum States from a Microwave Cavity Memory, Nat. Phys. 13, 882 (2017).
  * (7) A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S.M. Girvin, S. Shankar, and M. H. Devoret, Stabilization and operation of a Kerr-cat qubit, Nature 584, 205 (2020).
  * (8) J. Gambetta, A. Blais, M. Boissonneault, A. A. Houck, D. I. Schuster, and S. M. Girvin, Quantum Trajectory Approach to Circuit QED: Quantum Jumps and the Zeno Effect, Phys. Rev. A 77, 012112 (2008).


# Random Projections of Sparse Adjacency Matrices

Frank Qiu Statistics Department; University of California, Berkeley

###### Abstract

We analyze a random projection method for adjacency matrices, studying its
utility in representing sparse graphs. We show that these random projections
retain the functionality of their underlying adjacency matrices while having
extra properties that make them attractive as dynamic graph representations.
In particular, they can represent graphs of different sizes and vertex sets in
the same space, allowing for the aggregation and manipulation of graphs in a
unified manner. We also provide results on how the size of the projections
need to scale in order to preserve accurate graph operations, showing that the
size of the projections can scale linearly with the number of vertices while
accurately retaining first-order graph information. We conclude by
characterizing our random projection as a distance-preserving map of adjacency
matrices analogous to the usual Johnson-Lindenstrauss map.

##  1 Introduction

The adjacency matrix is a popular and flexible graph representation, encoding
a graph’s edgeset in an explicit and easy to access manner. Furthermore, many
natural graph operations correspond to linear algebraic operations on the
adjacency matrix, such as edge addition-deletion, edge or vertex queries, and
edge composition. The properties of the adjacency matrix and related
quantities like the graph Laplacian also yield important insights about the
underlying graph [19] [1].

However, a major defect of the adjacency matrix is that its size scales
quadratically with the number of vertices. In many real world applications
where the underlying graph is extremely large, this presents a fundamental
problem. For example, graphs in the Stanford Large Network Collection include
graphs derived from social, communication, and road networks that each have
millions of vertices. The adjacency matrices of such large graphs would
require trillions of parameters to store and work with, which is infeasible in
practice. Hence, there has been much interest in graph compression techniques.
There is a particular interest in sparse graph techniques, since a majority of
real world graphs tend to have sparsely connected vertices.

In a previous work [18], we introduced a graph embedding method that
represents a graph’s edgeset as a random sum. This method can be viewed as a
random projection technique for adjacency matrices, and in this paper we
discuss its application to large sparse graphs. In particular, we show that
these random projections have the same space and operational time complexity
as popular sparse graph representations like the Compressed Sparse Row format
(CSR). Moreover, they retain all of the linear-algebraic graph operations of
adjacency matrices, giving a flexible and computationally-efficient
representation. Interestingly, the random projection described in this paper
is an inner-product preserving projection on the space of adjacency matrices,
analogous to the Johnson-Lindenstrauss map for vectors[7].

##  2 Related Work

We proposed our random projection method in previous work [18], framing it as
a graph embedding method in the spirit of hyperdimensional computing(HDC)
[8][10][6]. HDC graph representations encode vertices as high-dimensional
vectors, bind these vectors to generate edge embeddings, and sum these edge
embeddings to represent the entire edgeset: a bind-and-sum approach
[15][12][17][16][9]. Our method can be classified as a bind-and-sum method,
since we assign a random vector to each vertex and represent an edge by
binding the source and target vertices via the tensor(outer) product. However,
our method is purely a random projection and does not seek to learn a good
node or graph embedding, instead leveraging pseudo-orthogonality to maintain
properties of the adjacency matrix.

Our randomized projections can also be viewed as a ”reverse” adjacency matrix
factorization. These techniques seek to find an informative and compressed
factorization of the adjacency matrix (or adjacency tensor for typed graphs).
For example, the RESCAL algorithm [14] applied to a $n\times n$ adjacency
matrix $A$ learns a $d\times n$ entity matrix $E$ and $d\times d$ relational
matrix $R$ such that $A\approx ERE^{T}$. Here, the columns of $E$ are the
entity (vertex) embeddings that encode information about each entity, and the
matrix $R$ encodes information about entity-entity interactions. We work with
the source-target factorization of an adjacency matrix $A=CD^{T}$, where $C$
and $D$ are $n\times k$ matrices whose $i^{th}$ columns are the respective
coordinate vectors of the source and target vertex of the $i^{th}$ edge. Our
method encodes each vertex as a random vector and then replaces the columns of
$C$ and $D$ with their corresponding random vertex codes.

As a random projection of matrices, our method falls within the field of
random numerical linear algebra [13]. These are random algorithms and
techniques for solving problems in linear algebra, such as solving linear
systems and finding maximum/minimum eigenvalues. These techniques aim to offer
advantages in speed and memory over their deterministic counterparts, and our
method aims to provide both in representing and working with a sparse
adjacency matrix. Among these algorithms, the FastRP algorithm[4] is another
example of a random projection method. It computes node embeddings by taking a
weighted sum of the power of the adjacency matrix and then randomly projects
them into a $d\times n$ matrix whose columns form the node embeddings. While
the goal of FastRP is to generate node embeddings that capture multi-hop
information, our method seeks to informatively embed the adjacency matrix
itself.

##  3 Intuition Behind the Projection

Given directed graph $G=(V,E)$ and an arbitrary ordering of its vertex set
$V$, its adjacency matrix can be expressed as the sum of outer products of
coordinate vectors:

$A=\sum_{ij:A_{ij}=1}e_{i}e_{j}^{T}$

In this form, it is easy to see that many graph operations correspond to
linear operations on the adjacency matrix. For example, adding/deleting the
edge $(v_{s},v_{t})$ corresponds to adding/subtracting the matrix
$e_{s}e_{t}^{T}$ from $A$; querying if the graph contains the edge
$(v_{s},v_{t})$ corresponds to the product $e_{s}^{T}Ae_{t}$; the $k^{th}$
matrix powers correspond to $k$-length paths. One property these fundamental
graph operations share is that they only require orthonormality of the vertex
vectors. If we swap the coordinate vectors with any set of orthonormal
vectors, the transformed adjacency matrix retains all of its functionality.

For example, let $B=[b_{1}\cdots b_{n}]$ be any orthogonal matrix, whose
columns form an orthonormal basis. Changing bases to $B$, our transformed
adjacency matrix $A_{B}$ takes the form:

$A_{B}=\sum_{ij:A_{ij}=1}b_{i}b_{j}^{T}=BAB^{T}$

In this form, we can perform all the usual adjacency matrix operations by
swapping coordinate vectors with their counterparts in $B$. For example,
adding/deleting the edge $(v_{s},v_{t})$ now corresponds to adding/subtracting
the matrix $b_{s}b_{t}^{T}$ from $A_{B}$; the edge query of $(v_{s},v_{t})$
becomes $b_{s}A_{B}^{T}b_{t}$; matrix powers still correspond to finite length
paths in the sense that $b_{i}^{T}A_{B}^{k}b_{j}=1$ if and only if there is a
$k$-length path from $v_{i}$ to $v_{j}$.

From the above discussion, we see that orthonormality is the key property
required for exact graph functionality. Hence, relaxing to approximate
orthonormality allows us to compress adjacency matrices while retaining graph
functionality in a minimally noisy manner. We make use of random pseudo-
orthonormal vectors, which are vectors whose dot products are negligible with
high probability. While we can pack at most $d$ orthonormal vectors in a
$d$-dimensional space, we can pack many more pseudo-orthonormal vectors in the
same space. This allows us to compress adjacency matrices, and in the case of
sparse matrices this compression effect is especially pronounced.

##  4 Random Projection of Adjacency Matrices

Consider a directed graph $G=(V,E)$ with $|V|=n$. An ordering of the vertex
set $V$ induces the $n\times n$ adjacency matrix $A$:

$A_{ij}=\begin{cases}1\quad\text{if $(v_{i},v_{j})\in E$}\\\
0\quad\text{otherwise}\end{cases}$

We then construct a random $d\times n$ projection matrix $P_{V}$, where the
columns of $P_{V}$ are sampled i.i.d. from the uniform distribution on the
$d$-dimension unit sphere $\mathbb{S}^{d-1}$. The random projection
$\pi_{V}(A)$ is then given by the following equation:

$\pi_{V}(A)=P_{V}AP_{V}^{T}$

Expressing the adjacency matrices as the sum of outer products between
coordinate vectors $A=\sum e_{i}e_{j}^{T}$, the projection swaps the $i^{th}$
coordinate vector with the $i^{th}$ column $p_{i}$ of $P_{V}$: $\sum
e_{i}e_{j}^{T}\mapsto\sum p_{i}p_{j}^{T}$. In light of the discussion of
Section 3, this can be viewed as a pseudo-orthogonal basis change.

###  4.1 Graph Operations on Random Projections

We saw in Section 3 that operations in a new orthonormal basis require
swapping coordinate vectors with the new basis vectors. Similarly, we can
perform all the usual graph operations with the projection $\pi_{V}(A)$ by
substituting the appropriate random code for each vertex. For example, to
query if the edge $(v_{i},v_{j})$ is in the edgeset, we would usually compute
the product $e_{i}^{T}Ae_{j}$. This returns a 1 if $(v_{i},v_{j})$ is an edge
and 0 otherwise. The randomized analogue then takes the form
$p_{i}^{T}\pi_{V}(A)p_{j}$. This quantity is close to 1 with high probability
if $(v_{i},v_{j})$ is an edge and close to 0 with high probability otherwise.

In general, every graph operation has an analogue on the projected matrices by
substituting in the random vertex codes. Rather than using the coordinate
vector $e_{i}$ associated with vertex $v_{i}$, when working with the
randomized projection $\pi_{V}(A)$ we use the $i^{th}$ column vector $p_{i}$
of the projection matrix $P_{V}$ instead.

###  4.2 Changing the Vertex Set and Graph Aggregation

One attractive property about our random projections is that they transform
nicely under changes to the underlying vertex set. Suppose we wish to expand
our graph by adding a new vertex $v_{n+1}$ and new edges
$\\{(v_{i},v_{n+1})\\}_{i}$. With the usual $n\times n$ adjacency matrix, this
would require expanding to a $n+1\times n+1$ matrix and filling in the
appropriate entries of the added column and row. However, for the projected
matrix we need to only generate a new random vector $p_{n+1}$ and add the
appropriate edges to the projection: $\pi_{V}(A)+\sum_{i}p_{i}p_{n+1}^{T}$.

Our random matrix projection can be viewed as a projection method for the set
of graphs whose vertices are in a fixed vertex set $V$. In light of the above
discussion, this projection can be naturally extended to sets containing $V$
or restricted to subsets of $V$ by expanding or restricting the set of random
vertex codes respectively. This property is particularly attractive for
dynamic graph representations, where the edge and vertex set of the graph
change over time. The dimension of the projected matrices stays the same under
changes to the vertex set, and the addition/deletion of edges involving new
vertices is simple. We only need to keep in mind the total capacity of the
projection space: how many edges we can add to a projected matrix before
accuracy in graph operations begins to break down. We analyze this behavior in
Section 5.

###  4.3 Translation between Different Projections

Once we assign a random code to each vertex and construct our projection
matrix $P_{V}$, that vertex-vector codebook is fixed. However, we might desire
to reassign a new random vector to each vertex. This translation procedure has
a simple analogue for projected matrices. Suppose we have a vertex set $V$
with two different random projection matrices $P_{V}$ and $Q_{V}$. The
$i^{th}$ columns of $P_{V}$ and $Q_{V}$ are the random vectors assigned to
vertex $v_{i}$ under each projection respectively. In order to swap the vector
$p_{i}$ with $q_{i}$, we construct the translation matrix
$T_{(P,Q)}=Q_{V}P_{V}^{T}$. By pseudo-orthonormality, we see that
$T_{(P,Q)}p_{i}\approx q_{i}$ for every $i$. Hence, if $A_{P}$ is the
projection of $A$ under $P_{V}$ and $A_{Q}$ is the projection under $Q_{V}$,
we have the following relation:

$A_{Q}\approx T_{(P,Q)}AT_{(P,Q)}^{T}$

###  4.4 Graph Subsets and Aggregation

Given a subset $S\subseteq V$, the subgraph generated generated by $S$,
denoted $G_{S}\subset G$, is the subgraph generated by all edges involving
vertices in $S$. Given the projected adjacency matrix $\pi(A_{G})$, we can
extract the projection of the subgraph $\pi(A_{G_{S}})$ by the following
procedure. Let $P_{S}$ be the $d\times|S|$ matrix whose columns are the random
codes assigned to each vertex of $S$, and define $T_{S}=P_{S}P_{S}^{T}$. By
pseudo-orthonormality, we have the following approximate relation:

$\pi(A_{G_{S}})\approx T_{S}\pi(A_{G})T_{S}^{T}$

This subset procedure, along with graph translation, can be applied to
aggregate multiple graphs into a single large graph. For example, suppose we
have two disjoint graphs $G$ and $H$. Their aggregate graph $G\cup
H=(V_{G}\cup V_{H},E_{G}\cup E_{H})$ is the result of combining their vertex
and edge sets together. We saw earlier that edge addition corresponds to
adding the projected edges to our projected matrices. Therefore, given
projections $\pi(A_{G})$ and $\pi(A_{H})$, the projection of their aggregate
is the sum of their projections:

$\pi(A_{G}\cup A_{H})=\pi(A_{G})+\pi(A_{H})$

This can be combined in tandem with the subsetting procedure to combine
selected subgraphs from a collection of graphs. Even when their vertex codes
are different, the translation procedure mentioned previously allows use to
translate them all into a single code before aggregation.

One interesting application of this suite of procedures is allowing a divide-
and-conquer approach to storing a large sparse graph. Graph operations on the
projected matrices depend on the size of the projection space (see Section 5),
so breaking a large graph into subgraphs allows us to store their projections
in a smaller projection space. Operations involving an individual subgraph
also happen in a lower-dimensional space and have decreased time complexity.
However, when the need arises to operate on information from a pool of these
subgraphs, we can use the above subset and pooling procedures to easily
generate a wide suite of projections that correspond to those generated by
their graph aggregates.

The above examples demonstrate that we can represent graphs of varying size
and different vertex sets in the same projection space. Translation and
pooling operations are of fixed dimension, allowing for a unified way to
manipulate and combine a broad range of graphs. This suggests that random
projections are also ideal for applications where graph aggregation is an
important operation such as the divid-and-conquer situation described above.

##  5 Scaling Properties of Random Projections

We now study how the size of the random projection space needs to scale with
the underlying graph. Importantly, we show that the size scales with size of
the edgeset rather than the vertex set. This makes our method particularly
amenable to sparse graphs, where the size of the edgeset is proportional to
the size of the vertex set.

###  5.1 $m$-Order Graph Operations

First, we need to define $m$-order graph operations. This is important because
each order requires a different scaling of the projection space. We say a
graph operation has order $\boldsymbol{m}$ if it can be expressed as an
operation involving the $m^{th}$ power of the adjacency matrix. For example,
the edge query is a first order graph operation because it is a function of
the adjacency matrix: $q((v_{i},v_{j}),G)=e_{i}^{T}Ae_{j}$.

###  5.2 First and Second Order Scaling

We first characterize how the size of the projection needs to scale in order
to retain accurate first and second order graph operations. Intuitively, these
correspond to edge information of the 1-hop and 2-hop neighborhoods of the
graph respectively. We give two informal statements on how the dimension $d$
of the projection space $\mathbb{R}^{d\times d}$ must scale with the number of
edges, subject to a constraint on the vertex connectivity (number of edges
each vertex participates in). An account of the technical results justifying
these statements is given in Appendix A.

######  Property 1 (First Order Scaling).

Let $G$ be a graph with $k$ edges with maximum node degree $O(\sqrt{k})$. A
random projection into $\mathbb{R}^{d\times d}$ needs to have
$d=\Omega(\sqrt{k})$ in order to retain accurate first order operations.

######  Property 2 (Second Order Scaling).

Let $G$ be a graph with $k$ edges with maximum node degree $O(k^{1/3})$. A
random projection into $\mathbb{R}^{d\times d}$ needs to have
$d=\Omega(k^{2/3})$ in order to retain accurate second order operations.

For sparse graphs, the node degree condition is satisfied for all or a
majority of vertices. Indeed, a closer inspection of the proofs in Appendix A
show that as long as the vertices of the relevant edges satisfy this or a much
looser bound on the node degree, the results still hold. This allows us to
include sparse graphs that might have central vertices, which act like hubs
and connect to many other vertices. If we are only interested in first order
properties of the graph, the first result implies that we can compress a large
sparse adjacency matrix using $n^{2}$ parameters into a smaller matrix using
$\Omega(n)$ parameters, resulting a drastic compression for large sparse
graphs.

###  5.3 $m$-Order Scaling

While first and possibly second order operations comprise a bulk of the
interesting graph operations, for completeness we have the following informal
scaling statement for $m$-order graph operations.

######  Property 3 ($m$-Order Scaling).

Let $G$ be a graph with $k$ edges with maximum node degree $O(k^{1/(m+1)})$. A
random projection into $\mathbb{R}^{d\times d}$ needs to have
$d=\Omega(k^{\frac{m}{m+1}})$ in order to retain accurate $m$-order
operations.

##  6 Edge Representations of Sparse Matrices

In this section, for various sparse matrix representations we analyze both
their size complexity as well as the time complexity of graph operations. We
also consider the time complexity of numerical sparse matrix methods when
appropriate. In this manner, we hope to contextualize both the strengths and
weaknesses of our random projections relative to other sparse matrix
representations. Table 1 summarizes the results of this section.

###  6.1 Alternate Sparse Matrix Representations

We first consider the coordinate list representation. This represents a sparse
matrix as a list of 3-tuples $(r,c,v)$ corresponding to each non-zero entry,
where $r$ and $c$ are the row and column indices while $v$ is the value. A
related representation is the Dictionary of Keys (DoK) representation, which
is similar to the coordinate list representation except now each row-column
index $(r,c)$ serves as a key with value $v$. Finally, the CSR format
represents a matrix using three arrays containing the non-zero values, their
column indices, and the number of non-zero entries above each row
respectively. The CSR format can be easily computed from the previous two
representations and vice versa.

###  6.2 Space Complexity

A sparse graph with $n$ vertices has $O(n)$ edges. Hence, from their
descriptions all three of the alternate sparse matrix representations require
$O(n)$ parameters. For random projections, if we are only concerned with
representing first-order properties, Theorem A.1 establishes that the $d\times
d$ projections must have $d=\Omega(n)$, meaning the projections are of size
$O(n)$.

###  6.3 Graph Operation Complexity

We analyze some graph operations where the complexity differs based on the
representation. For data structure operations, we use their complexity in
Python as stated in the Python Wiki. Throughout this section, our random
projection are $d\times d$ matrices.

####  6.3.1 Edge Query

To look up if an edge $(v_{i},v_{j})$ exists in a coordinate list, we merely
check if the list contains a tuple whose first two entries are $(i,j)$. Since
the coordinate list has length $O(n)$, list lookup has complexity $O(n)$. As a
three list format, the CSR representation has the same complexity. However,
the ”in” operator for dictionaries is $O(1)$, meaning edge lookup in DoK is
$O(1)$. The edge query for random projections is the product
$p_{i}^{T}\pi_{V}(A)p_{j}$, which has $O(d^{2})$. If our project space is
calibrated to preserve only first-order operations, then this edge query has
complexity $O(n)$.

####  6.3.2 Edge Composition

Edge composition is a complicated case. For each of the sparse matrix
representatiosn, edge compositions requires a nested procedure where, for each
edge $(i,j)$, one needs to identify all edges $(j,k)$ and return the edge
$(i,k)$, which is naively $O(n^{2})$. Importantly, this naive algorithm is not
parallelizable, making it especially painful to compute. Alternatively, edge
composition naturally corresponds to the second power of the adjacency matrix.
For a sparse adjacency matrix with $O(n)$ edges, numerical methods for sparse
matrices have complexity at most $O(n^{2})$ [2].

As for random projections, the naive algorithm for the multiplication of two
$d\times d$ matrices has complexity $O(d^{3})$[5] [11]. In order to retain
accuracy for a second-order operations, property 2 states that
$d=\Omega(n^{2/3})$. Thus, computing the second matrix power of our
projections has naive complexity $O(n^{2})$ while using methods like
Strassen’s algorithm[11] can speed it up even further.

###  6.4 Fast Numerical Linear Algebra

One important aspect of our random projections is they benefit from numerical
methods for linear algebra, such as parallel computation. This means that the
time complexity of graph operations can be substantially reduced through these
methods. The other sparse matrix representations do not enjoy these benefits,
as we saw with the edge composition example. Hence, in practice the time
complexity of graph operations with random projections can be significantly
reduced, and they benefit from advances in numerical linear algebra.

| Random Projections | DoK | CL | CSR | Sparse NLA
---|---|---|---|---|---
Space Complexity | $O(n)$ | $O(n)$ | $O(n)$ | $O(n)$ | N/A
Edge Query | $O(n)$ | $O(1)$ | $O(n)$ | $O(n)$ | N/A
Edge Composition | $O(n^{2})^{*}$ | $O(n^{2})$ | $O(n^{2})$ | $O(n^{2})$ | N/A
Matrix Addition | $O(n)$ | N/A | N/A | N/A | $O(n)$
Matrix Multiplication | $O(n^{2})^{*}$ | N/A | N/A | N/A | $O(n^{2})$
Table 1: Table of space/time complexity for various graph and matrix
operations. The underlying sparse graph has $n$ vertices and $O(n)$ edges.
From left to right, the methods considered are random projections, dictionary
of keys (DoK), coordinate list (CL), compressed sparse row (CSR), and sparse
numerical linear algebra (Sparse NLA). For matrix multiplication, we use the
time complexity of the naive algorithm rather than other advanced methods
[11]. We assume our random projections are calibrated to first order
operations except for edge composition/matrix multiplication($*$) where we
assume calibration to second order operations.

##  7 Johnson-Lindenstrauss Analogue for Graphs

In the introduction, we stated that our random projection method is a norm-
preserving map analogous to the one given by Johnson-Lindenstrauss (JL) lemma
for vectors. Here, we first discuss a natural inner product on the space of
adjacency matrices. We then state results showing how our random projections
preserve this inner product and its associated norm for a finite set of
adjacency matrices, characterizing our method as a JL map for the space of
graphs.

###  7.1 A Natural Inner Product for Adjacency Matrices

In our previous work [18], we noted that there is a natural inner product on
the space of adjacency matrices that counts the number of shared edges.
Consider two graphs $G$ and $H$ whose vertex sets are contained in some larger
set $V$. Given an ordering of $V$, we have the induced adjacency matrices
$A_{G}$ and $A_{H}$. Define the inner product between their adjacency matrices
as:

$\langle A_{G},A_{H}\rangle\coloneqq tr(A_{G}^{T}A_{H})$

Expressing each adjacency matrix as the sum of coordinate outer products, we
see that this function counts the number of edges common to both graphs:

$\displaystyle tr(A_{G}^{T}A_{H})$ $\displaystyle=tr([\sum_{(i,j)\in
E_{G}}e_{i}e_{j}^{T}]^{T}[\sum_{(k,l)\in E_{H}}e_{k}e_{l}^{T}])$
$\displaystyle=tr(\sum_{(i,j)\in E_{G},(j,l)\in E_{H}}e_{i}e_{l}^{T})$
$\displaystyle=\sum_{(i,j)\in E_{G}\cap E_{H}}1$ $\displaystyle=|E_{G}\cap
E_{H}|$

One can check this is indeed an inner product, and in fact it generates the
Frobenius norm: $tr(A^{T}A)=||A||_{F}^{2}$.

Note that we assumed the vertex sets of both graphs were contained in some
larger set $V$. This is not an issue for a finite set of graphs $G_{i}$ with
finite vertex sets $V_{i}$, as we can define $V=\cup V_{i}$. The adjacency
matrix of a graph with respect to this large set $V$ is a block matrix whose
single block is equal to the original adjacency matrix.

###  7.2 Random Projections are a JL Map

In light of graph inner product defined in Section 7.1, we almost have a JL-
type result since the Frobenius norm of matrices coincides with the Euclidean
norm if we regard matrices in $\mathbb{R}^{d\times d}$ as vectors in
$\mathbb{R}^{d^{2}}$. We confirm this by proving that, with high probability,
our random projection method is a map that satisfies the conditions of the JL
Lemma.

######  Theorem 7.1 (Random Projections are a JL map).

Consider a set of $N$ graphs with adjacency matrices $A_{1},\cdots,A_{n}$. For
small $\epsilon>0$ and $d=\Omega(\frac{\log N}{\epsilon^{2}})$, with high
probability our random projections satisfy:

$(1-\epsilon)||A_{i}-A_{j}||^{2}\leq||\pi(A_{i})-\pi(A_{j})||^{2}\leq(1+\epsilon)||A_{i}-A_{j}||^{2}\quad\forall
i,j$ (1)

##  8 Discussion

We presented a random projection method for sparse adjacency matrices. These
random projections exploit pseudo-orthonormality of random vectors to
drastically compress the adjacency matrix while still retaining all of its
functionality. While the exact scaling depends on the graph properties we wish
to preserve, Theorem A.1 shows that we can compress $n\times n$ sparse
matrices into $\Omega(\sqrt{n})\times\Omega(\sqrt{n})$ matrices while
preserving first order graph operations. These random projections also enjoy
properties that their underlying adjacency matrices do not. Sections 4.2 and
4.3 show that random projection allow us to represent graphs of varying size
and different vertex sets in the same projection space. This common space is
equipped with modification and aggregation operations that apply to all graphs
in the space, and random projections provide a unified way for representing
and working with graphs. The complexity analysis of 6 shows that these random
projections are competitive with existing sparse matrix representations and
numerical methods, and they can take advantage of numerical techniques for
speeding up linear algebra operations. All these properties suggest that
random adjacency matrix projections are a dynamic, flexible, and expressive
graph compression technique well suited to a variety of applications where
large sparse matrices occur. One interesting application of our random
compression technique would be to the field of graph neural networks [20].
Such networks use the adjacency matrix during the linear portion of its
computations, and random projections could help extend such networks to both
large sparse graphs as well as time-varying, dynamic graphs.

## References

  * [1] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural Computation, 15:1373–1396, 2003.
  * [2] Keivan Borna and Sohrab Fard.  A note on the multiplication of sparse matrices.  Open Computer Science, 4(1):1–11, 2014.
  * [3] Stéphane Boucheron, Gábor Lugosi, and Olivier Bousquet.  Concentration Inequalities, pages 208–240.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.
  * [4] Haochen Chen, Syed Fahad Sultan, Yingtao Tian, Muhao Chen, and Steven Skiena.  Fast and accurate network embeddings via very sparse random projection, 2019.
  * [5] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein.  Introduction to Algorithms, Third Edition.  The MIT Press, 3rd edition, 2009.
  * [6] Ross W. Gayler and Simon D. Levy.  A distributed basis for analogical mapping.  2009\.
  * [7] William B. Johnson and Joram Lindenstrauss.  Extensions of lipschitz mappings into a hilbert space.  Contemporary Mathematics, 26, 1984.
  * [8] Pentti Kanerva.  Sparse Distributed Memory.  MIT Press, Cambridge, MA, USA, 1988.
  * [9] Jaeyoung Kang, Minxuan Zhou, Abhinav Bhansali, Weihong Xu, Anthony Thomas, and Tajana Rosing.  Relhd: A graph-based learning on fefet with hyperdimensional computing.  2022 IEEE 40th International Conference on Computer Design (ICCD), pages 553–560, 2022.
  * [10] Denis Kleyko, Dmitri A. Rachkovskij, Evgeny Osipov, and Abbas Jawdat Rahim.  A survey on hyperdimensional computing aka vector symbolic architectures, part ii: Applications, cognitive models, and challenges.  ArXiv, abs/2112.15424, 2021.
  * [11] Yan Li, Sheng-Long Hu, Jie Wang, and Zheng-Hai Huang.  An introduction to the computational complexity of matrix multiplication.  Journal of the Operations Research Society of China, 8(1):29–43, 2020.
  * [12] Yunpu Ma, Marcel Hildebrandt, Volker Tresp, and Stephan Baier.  Holistic representations for memorization and inference.  In UAI, 2018.
  * [13] Per-Gunnar Martinsson and Joel Tropp.  Randomized numerical linear algebra: Foundations and algorithms, 2021\.
  * [14] Maximilian Nickel, Xueyan Jiang, and Volker Tresp.  Reducing the rank in relational factorization models by including observable patterns.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
  * [15] Maximilian Nickel, Lorenzo Rosasco, and Tomaso A. Poggio.  Holographic embeddings of knowledge graphs.  In AAAI, 2016.
  * [16] Igor O. Nunes, Mike Heddes, Tony Givargis, Alexandru Nicolau, and Alexander V. Veidenbaum.  Graphhd: Efficient graph classification using hyperdimensional computing.  2022 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 1485–1490, 2022.
  * [17] Prathyush Poduval, Haleh Alimohamadi, Ali Zakeri, Farhad Imani, M. Hassan Najafi, Tony Givargis, and Mohsen Imani.  Graphd: Graph-based hyperdimensional memorization for brain-like cognitive learning.  Frontiers in Neuroscience, 16, 2022.
  * [18] Frank Qiu.  Graph embeddings via tensor products and approximately orthonormal codes, 2022.
  * [19] Ulrike von Luxburg.  A tutorial on spectral clustering, 2007.
  * [20] Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph neural networks: A review of methods and applications.  AI Open, 1:57–81, 2020.

##  Appendix A Scaling Theorems

We provide theorems justifying the properties listed Section 5. These are
mainly an adaptation of results found in our previous work [18]. Throughout
this section, we abuse notation and use the same symbol $u$ to denote both the
vertex and its random code. Similarly, we use the same symbol $G$ to denote
the graph and its random projection.

###  A.1 Auxiliary Lemma for Spherical Random Vectors

We will need the following lemma [18].

######  Lemma A.1.

Let $X$ be the dot product between two vectors sampled uniformly and
independently from the $d$-dimensional unit sphere $\mathbb{S}^{d-1}$. Then
$E(X)=0$ and $E(X^{2})=\frac{1}{d}$.

###  A.2 Property 1: First-Order Scaling

The following theorem and its proof is an abbreviated adaptation of Theorems
10.7 in our previous work [18]. Given an edge $(v_{i},v_{j})$ and adjacency
matrix $A$, the edge query $e_{i}^{T}Ae_{j}$ returns 1 if $(v_{i},v_{j})$ is
an edge of the graph and 0 otherwise. Our goal is prove the edge query
analogue for our random projection returns the correct value with high
probability, showing that our random projection accurately retains first-order
edge information. An edge query is a true query if the queried edge exists in
the graph and a false query if it doesn’t.

######  Theorem A.1.

Let $G$ be a graph with $k+1$ edges and maximum node degree $\frac{l}{2}$. For
a random projection into $\mathbb{R}^{d\times d}$, we have the following:

  1. 1.

If $T$ denotes the result of a true query, then:

$P(|T-1|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\frac{l}{d}+\frac{k-l}{d^{2}}+\frac{\epsilon}{3})})$
(2)

  2. 2.

If $F$ denotes the result of a false query, then:

$P(|F|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\frac{l}{d}+\frac{k+1-l}{d^{2}}+\frac{\epsilon}{3})})$
(3)

###### Proof.

Consider the random projection of a graph with $k+1$ edges:

$G=uv^{T}+\sum_{i=1}^{k}p_{i}q_{i}^{T}$ (4)

where $u$, $v$, $p_{i}$, $q_{i}$ are random vectors drawn i.i.d. from the
uniform distribution on the $d$-dimensional unit sphere.

We first prove equation 2 by considering the result of querying our random
projection for the exist of the edge $(u,v)$, which is a valid edge. We assume
$m$ of the other edges share one vertex with $(u,v)$, and WLOG we assume they
all share vertex $u$. The result of true query $T$ can be expressed as:

$T=u^{T}Gv=1+\sum_{i=1}^{m}\langle p_{i},v\rangle+\sum_{j=1}^{k-m}\langle
u,q_{j}\rangle\langle v,r_{j}\rangle$ (5)

From Lemma A.1, both sums have mean 0 with the first sum having variance
$\frac{m}{d}$ and the second having variance $\frac{k-m}{d^{2}}$. Hence, the
total variance is $\frac{m}{d}+\frac{k-m}{d^{2}}$, and applying Bernstein’s
inequality gives:

$P(|T-1|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\frac{m}{d}+\frac{k-m}{d^{2}}+\frac{\epsilon}{3})})$

By assumption, we can bound $m\leq\frac{l}{2}$, and plugging in the worst case
of $m=\frac{l}{2}$ gives equation 2.

Similarly, now suppose we query by a spurious edge $(s,t)$. Assuming $q$ of
the other edges share one vertex with $(s,t)$, we use the same argument as
above to derive the false query bound given by equation 3. ∎

Letting $Q$ denote a general edge query and $\sigma^{2}=Var(Q)$, Theorem A.1
can be summarized as:

$P(|Q-EQ|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\sigma^{2}+\frac{\epsilon}{3})})$

Rewriting this bound in terms of $\sigma$:

$P(|Q-EQ|>C\sigma)\leq
2\exp(\frac{C^{2}\sigma^{2}}{2\sigma^{2}+\frac{2\sigma}{3}})\approx
2\exp(\frac{-3C^{2}}{2\sigma})$

If $l=O(\sqrt{k})$ then $\sigma^{2}=\frac{O(\sqrt{k})}{d}+\frac{O(k)}{d^{2}}$.
Hence, if $d=\Omega(\sqrt{k})$ then $Q$ is close to $EQ$ with high
probability. Since $EQ=1$ for a true query and $EQ=0$ for a false query, $Q$
is close to the correct value with high probability. This justifies the
statement of property 1.

###  A.3 Property 2: Second-Order Scaling

This following theorem and its proof is an abbreviated adaptation of Theorems
11.7 from our previous work [18]. Two edges are composable if the target
vertex of one is the source vertex of another, and an edge $(u,v)$ is in the
second power of the adjacency matrix if and only if it is the composition of
two composable edges. We need to show that the second matrix power accurately
represents edge information. To this end, we prove the edge query returns the
correct result with high probability when applied to the second power of the
random projection.

######  Theorem A.2.

Let $G$ be a graph with composable edges $(u,v)$ and $(v,w)$ along with $k-2$
nuisance edges, and assume $G$ has maximum node degree $\frac{l}{2}$. For a
random projection into $\mathbb{R}^{d\times d}$, we have the following results
for edge queries involving projection’s second matrix power:

  1. 1.

If $T$ denotes the result of a true query of the second matrix power, then:

$P(|T-1|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\frac{2l+1}{d}+\frac{kl-3l-3}{d^{2}}+\frac{k^{2}-k(l+2)+l+1}{d^{3}}+\frac{\epsilon}{3})}))$
(6)

  2. 2.

If $F$ denotes the result of a false query of the second matrix power, then:

$P(|F|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\frac{kl}{d^{2}}+\frac{k^{2}-kl}{d^{3}}+\frac{\epsilon}{3})}))$
(7)

###### Proof.

For the nuisance edge, assume $m_{1}$ have common vertex $v$ and $m_{2}$ have
common source $u$ or target $w$, and let $m=m_{1}+m_{2}$. WLOG we may assume
all $m_{2}$ edges have common source vertex $u$, and our projection can be
expressed as:

$G=uv^{T}+vw^{T}+\sum_{i=1}^{m_{1}}vp_{i}^{T}+\sum_{j=1}^{m_{2}}up_{j}^{T}+\sum_{l=1}^{k-m-2}q_{l}r_{l}^{T}$
(8)

From equation 8, we see the second matrix power should only contain the edge
$(u,w)$.

We first prove equation 6 and query the second matrix power $G^{2}$ for the
presence of the edge $(u,w)$. Letting $\epsilon$ denote the dot product of two
independent spherical vectors, the result of this edge query is a sum of terms
that are products of i.i.d. $\epsilon$’s. Let $n_{1}=k(m_{2})-m_{2}-m-3$ and
$n_{2}=k^{2}-k(m_{2}+2)+m_{2}+1$. Grouping terms by how many $\epsilon$’s they
contain, we write our true query $T$ as:

$T=1+\sum_{i=1}^{m+1}\epsilon_{i}+\sum_{j=1}^{n_{1}}\epsilon_{j_{1}}\epsilon_{j_{2}}+\sum_{h=1}^{n_{2}}\epsilon_{h_{1}}\epsilon_{h_{2}}\epsilon_{h_{3}}$

Using Lemma A.1, we see the variance $\sigma^{2}$ of the error terms is
$\sigma^{2}=\frac{m+1}{d}+\frac{n_{1}}{d^{2}}+\frac{n_{3}}{d^{3}}$. An
application of Bernstein’s inequality gives:

$P(|T-1|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{2(\frac{m+1}{d}+\frac{n_{1}}{d^{2}}+\frac{n_{2}}{d^{3}}+\frac{\epsilon}{3})})$

Since $m_{1},m2\leq l$ and $m\leq 2l$, plugging in the worst case of
$m_{1}=m_{2}=l$ and $m=2l$ gives equation 6.

Similarly, now we consider querying $G^{2}$ by a spurious edge $(s,t)$. We
assume that $m$ of the edges have either common source vertex $s$ or target
vertex $t$. Then, a similar computation as above shows that we can express the
false query $F$ as:

$F=\sum_{i=1}^{km}\epsilon_{i_{1}}\epsilon_{i_{2}}+\sum_{j=1}^{k^{2}-km}\epsilon_{j_{1}}\epsilon_{j_{2}}\epsilon_{j_{3}}$

An application of Bernstein’s inequality and a worst-case bound gives equation
7. ∎

As with Theorem A.1, the bounds of Theorem A.2 can be expressed in terms of
the query variance $\sigma^{2}$:

$P(|Q-EQ|>C\sigma)\leq
2\exp(\frac{C^{2}\sigma^{2}}{2\sigma^{2}+\frac{2\sigma}{3}})\approx
2\exp(\frac{-3C^{2}}{2\sigma})$

If $l=O(k^{1/3})$, then in both cases all terms in the variance can be
expressed as powers of $O(\frac{k^{2/3}}{d})$. Hence, as long as
$d=\Omega(k^{2/3})$ then we see that the edge query $Q$ is close to its
correct value $EQ$ with high probability. This justifies the statement of
Property 2.

###  A.4 Property 3: $m$-Order Scaling

We give a proof sketch is adapted from Section 12.1 of our previous work[18].
We aim to analyze the accuracy recovering edge information from the $m^{th}$
matrix power of our random projections and how the dimension of the projection
space $\mathbb{R}^{d\times d}$ needs to scale with $m$.

As in the proof of Theorem A.2, let $\epsilon$ denote the dot product of
distinct spherical vectors. From the proofs of Theorems A.1 and A.2, we need
to control the edge query variance to ensure that true and false queries are
close to their expected values with high probability. Intuitively, if the node
degree is small then performing an edge query on the $m^{th}$ matrix power
will result in a majority of the $k^{m}$ error terms being products of $m+1$
independent $\epsilon$’s. Such terms will be mean 0 and variance
$\frac{1}{d^{m+1}}$. Since the true and false query scores are 1 and 0
respectively, for accurate recovery we need the edge query variance to be less
than 1. As a sum of independent terms, the variance of the noise term will be
approximately $\frac{k^{m}}{d^{m+1}}$, implying that
$d=\Omega(k^{\frac{m}{m+1}})$ in order to retain accuracy. The node degree
bound of $l=O(k^{1/(m+1)})$ ensures that the variance can be expressed as sums
of powers of $O(\frac{k^{m/(m+1)}}{d})$.

##  Appendix B JL Lemma for Graphs

Here, we aim to prove Theorem 7.1 by establishing intermediate results.
Importantly, along the way we show how our random projection method preserves
the graph inner product of Section 7.1.

###  B.1 Random Projections Preserve Inner Products

######  Theorem B.1.

Consider two graphs $G$ and $H$ with $n_{1}$ and $n_{2}$ edges respectively.
Suppose they have $k$ edges in common. Among all $n_{1}n_{2}$ pairs
$(e,e^{\prime})\in E_{G}\times E_{H}$, suppose $q$ of these edge pairs share
exactly one vertex. Let $\pi(G)$ and $\pi(H)$ denote the random projections of
their adjacency matrices $A_{G}$ and $A_{H}$ into $\mathbb{R}^{d\times d}$.
For any $\epsilon>0$, we have:

$P(|\langle\pi(G),\pi(H)\rangle-\langle A_{G},A_{H}\rangle|>\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{\frac{q}{d}+\frac{n_{1}n_{2}-k-q}{d^{2}}+\frac{\epsilon}{3}})$
(9)

###### Proof.

We can express their random projections as:

$\displaystyle\pi(G)=\sum^{k}_{i=1}a_{i}b^{T}+\sum_{j=1}^{n_{1}-k}c_{j}d_{j}^{T}$
$\displaystyle\pi(H)=\sum^{k}_{i=1}a_{i}b^{T}+\sum_{l=1}^{n_{2}-k}e_{l}f_{l}^{T}$

where $a,b,c,d,e,f$ are random spherical vectors. Computing their inner
product:

$tr(\pi(G)^{T}\pi(H))=k+tr(\pi(G)^{T}\sum_{l=1}^{n_{2}-k}e_{l}f_{l}^{T})+tr(\pi(H)^{T}\sum_{j=1}^{n_{1}-k}c_{j}d_{j}^{T})=k+E$
(10)

We proceed to bound the error term $E$ in equation 10. By assumption, $q$ of
the terms in $E$ share one vertex with the remaining terms sharing no
vertices:

$\displaystyle
E=\sum_{i=1}^{q}\epsilon_{i}+\sum_{j=1}^{n_{1}n_{2}-k-q}\epsilon_{j}\epsilon^{\prime}_{j}$

where each $\epsilon$ is the dot product of independent spherical vectors.
Using Lemma A.1 and Bernstein’s inequality[3], we get the following bound:

$P(|E|\geq\epsilon)\leq
2\exp(\frac{-\epsilon^{2}}{\frac{q}{d}+\frac{n_{1}n_{2}-k-q}{d^{2}}+\frac{\epsilon}{3}})$

∎

###  B.2 Random Projections are a JL Map for Graphs

######  Theorem B.2.

Let $G$ be a graph with $k$ edges. For $d<k$ and small $\epsilon>0$, its
random projection into $\mathbb{R}^{d\times d}$ satisfies:

$P(|\lVert\pi(G)\lVert^{2}-k|>k\epsilon)\leq 2\exp{-d\epsilon^{2}}$ (11)

###### Proof.

Of the $k^{2}$ pairs $(e,e^{\prime})\in E_{G}\times E_{G}$, suppose that $q$
of them share exactly one vertex. Theorem B.1 states that:

$P(|\lVert\pi(G)\lVert^{2}-k|>k\epsilon)\leq
2\exp(\frac{-k^{2}\epsilon^{2}}{\frac{q}{d}+\frac{k^{2}-k-q}{d^{2}}+\frac{k\epsilon}{3}})$

As $q\leq k^{2}-k$, the worst case bound occurs when $q=k^{2}-k$ which gives
equation 11. Note that if $\epsilon$ is large, the third term
$\frac{k\epsilon}{3}$ dominates and gives the trivial bound
$2\exp(-k\epsilon)$. ∎

######  Theorem B.3.

Let $G$ and $H$ be two graphs with $\lVert A_{G}-A_{H}\lVert^{2}=m$. For $d<k$
and small $\epsilon>0$, the random projections of their adjacency matrices
into $\mathbb{R}^{d\times d}$ satisfies the following:

$P(|\lVert\pi(G)-\pi(H)\lVert-m|>m\epsilon)\leq 2\exp(-d\epsilon^{2})$ (12)

###### Proof.

If we include signed edges, the graph $G-H$ has $m$ total edges. An
application of Theorem B.2 gives the result. ∎

####  B.2.1 Proof of Theorem 7.1

###### Proof.

From Theorem B.3, our random projection satisfies the following inequalities
with high probability:

$(1-\epsilon)||A_{G}-A_{H}||^{2}\leq||\pi(G)-\pi(H)||^{2}\leq(1+\epsilon)||A_{G}-A_{H}||^{2}$
(13)

If we want this to hold for a set of $N$ sparse adjacency matrices, a union
bound over all $\binom{N}{2}$ pairs shows equation 13 holds with probability
at least $1-2\binom{N}{2}\exp{-d\epsilon^{2}}=1-N(N-1)\exp{-d\epsilon^{2}}$.
For a fixed probability threshold $T$ of violating the inequalities 13, let us
choose the optimal $d$ given $N$, denoted $d_{opt}$. That is, $d_{opt}$ is the
smallest integer $d$ such that $N(N-1)\exp{-d\epsilon^{2}}\leq T$. The optimal
$d_{opt}$ satisfies the following inequalities:

$\displaystyle N(N-1)\exp{-d_{opt}\epsilon^{2}}\leq T$ $\displaystyle
N(N-1)\exp{-(d_{opt}-1)\epsilon^{2}}>T$

Combining the two inequalities gives:

$\frac{\log T-2\log N}{\epsilon^{2}}\leq d_{opt}<\frac{\log T-2\log
N}{\epsilon^{2}}+1$

Thus, we have $d_{opt}=\Omega(\frac{\log N}{\epsilon^{2}})$, which matches the
scaling of the usual Johnson-Lindenstrauss lemma and establishes our random
projections as a JL map for graphs. ∎