# Compact fermionic quantum state preparation
with a natural-orbitalizing variational quantum eigensolving scheme

P. Besserve  Eviden Quantum Laboratory, 78340 Les Clayes-sous-Bois, France
CPHT, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91128
Palaiseau, France  M. Ferrero  CPHT, CNRS, Ecole Polytechnique, Institut
Polytechnique de Paris, 91128 Palaiseau, France  Collège de France, 11 place
Marcelin Berthelot, 75005 Paris, France  T. Ayral  Eviden Quantum Laboratory,
78340 Les Clayes-sous-Bois, France

###### Abstract

Assemblies of strongly interacting fermions, whether in a condensed-matter or
a quantum chemistry context, range amongst the most promising candidate
systems for which quantum computing platforms could provide an advantage.
Near-term quantum state preparation is typically realized by means of the
variational quantum eigensolver (VQE) algorithm. One of the main challenges to
a successful implementation of VQE lies in the sensitivity to noise exhibited
by deep variational circuits. On the other hand, sufficient depth must be
allowed to be able to reach a good approximation to the target state. In this
work, we present a refined VQE scheme that consists in topping VQE with state-
informed updates of the elementary fermionic modes (spin-orbitals). These
updates consist in moving to the natural-orbital basis of the current,
converged variational state, a basis we argue eases the task of state
preparation. We test the method on the Hubbard model in the presence of
experimentally relevant noise levels. For a fixed circuit structure, the
method is shown to enhance the capabilities of the circuit to reach a state
close to the target state without incurring too much overhead from shot noise.
Moreover, coupled with an adaptive VQE scheme that constructs the circuit on
the fly, we evidence reduced requirements on the depth of the circuit as the
orbitals get updated.

Many-body systems such as strongly interacting condensed-matter systems and
chemical compounds are expected to range amongst the first systems whose
understanding could benefit from quantum computational power [1]. Their
description with classical computers indeed always faces an exponential wall:
it may stem from the exponential size of the Hilbert space with the number
$N_{e}$ of particles, which strongly limits the scope of exact diagonalization
techniques (also known as full configuration interaction in chemistry), or an
exponentially growing variance in Monte Carlo methods due to the fermionic
sign problem, or an exponential bond dimension in tensor-network methods in
volume-law regimes.

The first prototypes of quantum computers could help overcome this exponential
barrier. Current prototypes, dubbed noisy, intermediate scale quantum (NISQ,
[2]) computers, feature high noise levels and low qubit counts. In the absence
of quantum error correction, the fidelity of quantum states decreases
exponentially with the circuit depth, placing very strong limitations on the
number of gates that can be executed reliably [3]. This precludes the use of
”textbook” quantum algorithms with provable speedups such as quantum phase
estimation, and ask for the development of algorithms that require only very
shallow circuits. This is typically achieved by doing most of the heavy-
lifting with classical computers (CPUs) and only spoon-feeding very specific
tasks to the quantum processing unit (QPU).

The epitome of such hybrid algorithms is the variational quantum eigensolver
(VQE [4], see [5, 6] for a review): one defines a parameterized circuit (also
called ansatz) and tunes its parameters by minimizing the energy, a quantum
observable over the circuit instance, with a classical optimizer. The success
of the (classical) parameter optimization is by no means guaranteed due to the
barren plateau problem ([7], see [8] for a review), which colloquially stems,
as the number of qubits gets large, from (i) a variational circuit which is
too expressive, or (ii) a level of noise which is too high, or (iii) an
observable which is too global, or (iv) an initial state that is too entangled
[9].

In this work, we address issues (i) and (ii) by focusing on a degree of
freedom specific to fermionic problems, namely the single-particle orbital
basis from which fermionic many-body states are constructed. There are two
ways to perform a single-particle orbital basis rotation
$\mathcal{U}^{\dagger}H\mathcal{U}$ on a Hamiltonian $H$, illustrated in Fig.
1. On the one hand, one can rotate, using so-called Givens rotations [10],
quantum states using a quantum circuit _on a quantum computer_. This approach
has been carried out e.g in [11] and is the inspiration of some fermionic
ansätze [12]. On the other hand, one can change the representation of the
Hamiltonian by performing $\mathcal{U}^{\dagger}H\mathcal{U}$ _on a classical
computer_ : since $\mathcal{U}$ corresponds to a single-particle basis change,
the corresponding computation is classically efficient (i.e polynomial in the
number of orbitals). We will dub this approach the ”off-chip” approach as it
is not performed on the quantum chip. Deciding on one or the other strategy
results from the following tradeoff: performing the basis rotations on the
quantum computer will render the computation more sensitive to quantum noise,
but can reduce the number of terms in the Hamiltonian and thus the measurement
overhead [13, 14]. Conversely, performing the basis change classically at the
level of the Hamiltonian could increase the number of terms and thus the
measurement overhead, but also help reduce the circuit depth and thus make it
more noise-robust.

Here, we focus on the second, ”off-chip” approach because we are focusing on
current (noisy) devices. Several such orbital optimization schemes have been
proposed in the past few years in various VQE or broader variational contexts:
mainly quantum chemistry [15, 16, 17, 18, 19, 20], but also nuclear physics
[21], tensor networks [22, 23] as well as neural networks [24]. In these
schemes, $\mathcal{U}$ is optimized through parameters that are put on an
equal footing with that of the variational circuit. These approaches either
add a large number of parameters to optimize, and/or a large number of
gradients to be evaluated in the case of adaptive ansatz construction [25].

Here, we present a lighter orbital-adapting technique consisting in looking
for a physically-motivated orbital basis instead of brute-force search of an
optimal basis. More specifically, we propose to leverage the so-called
natural-orbital (NO) basis, which has been extensively used for quantum
chemistry in a classical computing context.

We propose two different versions of our natural-orbitalization strategy. The
first version (already introduced by some of us in the context of Anderson
impurity models, [26]) consists in picking an ansatz circuit and simply
interlacing plain-vanilla VQE runs with rotations to the optimal state’s NO
basis. Doing so, we observe, for the Hubbard model, considerable increases in
the expressiveness of the ansatz at hand, rendering very shallow—and thus,
noise-resilient—ansätze capable of preparing the target state. For deep
circuits, noisy performances are also greatly improved. Interestingly, we find
that the increase in the number of Hamiltonian Pauli terms incurred by
rotating the single-particle basis seems to have a limited impact on the
measurement overhead due to the concentration of Pauli weights around a few
terms only.

The second version consists in combining the natural-orbitalization framework
with the ADAPT-VQE method [27] which constructs the circuit iteratively. This
allows to both ’learn’ a circuit and a representation. In this framework, we
observe that very shallow circuit ansätze can be discovered even in the
presence of noise.

Figure 1: ”On-chip” (purely quantum, left) vs ”off-chip” (classical, right)
strategy to rotate the single-particle basis of a variational ansatz
$\mathcal{W}(\bm{\theta})$.

##  I Variational quantum eigensolving with fermions: the role of the orbital
basis

###  I.1 Fermions and the orbital basis

We are interested finding the ground-state energy of a fermionic Hamiltonian
$H$ given in its second-quantized form

$H=\sum\limits_{p,q=1}^{M}h_{pq}c^{\dagger}_{p}c_{q}+\frac{1}{2}\sum\limits_{p,q,r,s=1}^{M}h_{pqrs}c^{\dagger}_{p}c^{\dagger}_{q}c_{r}c_{s},$
(1)

with $h_{pq}$ and $h_{pqrs}$ the kinetic and interaction tensors, and
$c^{\dagger}_{p},c_{p}$ the creation and annihilation operators for the $p$-th
spin-orbital ($p$ encompasses both a spatial and a spin degree of freedom).

Any many-body state $\ket{\psi}$ (including the ground state of $H$) can be
written as a linear combination

$\ket{\psi}=\sum_{n_{1}n_{2}\dots n_{M}\in\\{0,1\\}^{M}}\psi_{n_{1}n_{2}\dots
n_{M}}\ket{\phi_{n_{1},\dots,n_{M}}}$ (2)

of an exponential number of Fock states

$\ket{\phi_{n_{1},\dots,n_{M}}}=(c^{\dagger}_{1})^{n_{1}}\cdots(c^{\dagger}_{M})^{n_{M}}\ket{\bar{0}},$
(3)

which are multi-fermion states corresponding to the excitation of some of the
fermionic modes through creation operators $c^{\dagger}_{p}$ over the fermion
vacuum state $\ket{\bar{0}}$. There exists a freedom in the choice of this
basis of single-particle fermionic modes $c^{\dagger}_{p}$, namely one could
also write

$\ket{\psi}=\sum_{n_{1}n_{2}\dots
n_{M}\in\\{0,1\\}^{M}}\tilde{\psi}_{n_{1}n_{2}\dots
n_{M}}\ket{\tilde{\phi}_{n_{1},\dots,n_{M}}}$ (4)

with Fock states
$\ket{\tilde{\phi}_{n_{1},\dots,n_{M}}}=(\tilde{c}^{\dagger}_{1})^{n_{1}}\cdots(\tilde{c}^{\dagger}_{M})^{n_{M}}\ket{\bar{0}}$
built out of rotated fermionic modes:

$\tilde{c}^{\dagger}_{i}=\sum_{p}R_{pi}c^{\dagger}_{p},\;\;\;\tilde{c}_{i}=\sum_{p}R^{\dagger}_{ip}c_{p},$
(5)

with $R\in U(M)$ to ensure the preservation of the anticommutation rules.

The number of nonzero coefficients $\psi_{\alpha}$ (with
$\alpha=n_{1},\dots,n_{M}$) in the first representation is in general
different from the number of nonzero coefficients $\tilde{\psi}_{\beta}$. An
extreme case is the case when $|\psi_{0}\rangle$ can be represented as a
single Fock state in a basis (e.g the $c^{\dagger}_{p},c_{p}$ basis): only one
$\psi_{\alpha}$ coefficient is nonzero. The representation of such a state,
called a Slater determinant, can have many nonzero $\tilde{\psi}_{\beta}$
coefficients in another basis.

The rotation $R$ will lead to the following equivalent expression for the
Hamiltonian:

$H=\sum\limits_{i,j=1}^{M}\tilde{h}_{ij}\tilde{c}^{\dagger}_{i}\tilde{c}_{j}+\frac{1}{2}\sum\limits_{i,j,k,l=1}^{M}\tilde{h}_{ijkl}\tilde{c}^{\dagger}_{i}\tilde{c}^{\dagger}_{j}\tilde{c}_{k}\tilde{c}_{l},$
(6)

with modified tensors

$\displaystyle\tilde{h}_{ij}$
$\displaystyle=\sum_{pq}R^{\dagger}_{ip}h_{pq}R_{qj},$ (7a)
$\displaystyle\tilde{h}_{ijkl}$
$\displaystyle=\sum_{pqrs}R^{\dagger}_{ip}R^{\dagger}_{jq}h_{pqrs}R_{rk}R_{sl}.$
(7b)

One can thus seek to optimize the rotation basis $R$. For instance, one may
want to reduce the number of nonzero coefficients of the wavefunction, or,
alternatively, reduce the number of terms in the Hamiltonian (Eq. (6)).

Before we turn to the optimization per se, let us note that a change in $R$
will effect a change of coordinates

$\tilde{\psi}_{\beta}=\sum_{\alpha}\mathcal{U}_{\beta\alpha}\psi_{\alpha},$
(8)

with the expression of the many-body unitary transformation $\mathcal{U}\in
U(2^{M})$ as a function of $R$ given by Thouless’s theorem [28]

$\mathcal{U}=e^{-\sum_{pq}\left(\log R\right)_{pq}c^{\dagger}_{p}c_{q}}.$ (9)

In other words, an orbital change of basis is a unitary operation
$\mathcal{U}$ on the many-body wavefunction with the specific form given in
Eq. (9).

###  I.2 Optimizing the single-particle basis

In a standard VQE workflow, one seeks to optimize the variational energy

$E(\bm{\theta})=\langle\bar{0}|\mathcal{W}(\bm{\theta})^{\dagger}H\mathcal{W}(\bm{\theta})|\bar{0}\rangle,$
(10)

with $\mathcal{W}(\bm{\theta})$ a parameterized unitary implemented with
quantum gates. Adding an orbital optimization amounts to replacing the energy
function by

$E(\bm{\theta},\bm{\theta^{\prime}})=\langle\bar{0}|\mathcal{W}(\bm{\theta})^{\dagger}\mathcal{U}(\bm{\theta^{\prime}})^{\dagger}H\mathcal{U}(\bm{\theta^{\prime}})\mathcal{W}(\bm{\theta})|\bar{0}\rangle,$
(11)

with $\bm{\theta^{\prime}}$ a vector of the $M(M-1)/2$ parameters defining the
anti-hermitian matrix $\log R$ introduced in Eq. (9).

Many VQE schemes exploiting this freedom in the choice of single-particle
modes have been proposed in the past few years. They differ in (i) whether the
rotation is done classically (”off-chip”, at the level of the Hamiltonian) or
quantumly (”on-chip”, at the level of the state), as illustrated in Fig. 1,
and (ii) what they attempt to optimize with this new degree of freedom.

”Off-chip” schemes [17, 16, 18, 24, 20, 19, 21] consist in dressing the
Hamiltonian on a classical computer

$\tilde{H}(\bm{\theta^{\prime}})=\mathcal{U}(\bm{\theta^{\prime}})^{\dagger}H\mathcal{U}(\bm{\theta^{\prime}}).$
(12)

and preparing the wavefunction $\mathcal{W}(\bm{\theta})|0\rangle$ on a
quantum computer. This is a classically efficient (polynomial) operation
thanks to Thouless’s theorem.

One optimization criterion is to minimize of the one-norm of the Hamiltonian
[15]: after decomposing $\tilde{H}(\bm{\theta^{\prime}})$ in terms of Pauli
words $\tilde{H}(\bm{\theta^{\prime}})=\sum_{d}\lambda_{d}H_{d}$ given an
encoding, its one-norm is defined as
$\left\|\tilde{H}(\bm{\theta^{\prime}})\right\|_{1}=\sum_{d}|\lambda_{d}|$.
The standard error $\Delta E$ stemming from a limited number
$n_{\mathrm{shots}}$ of shots is upper-bounded by
$\left\|\tilde{H}(\bm{\theta^{\prime}})\right\|_{1}/\sqrt{n_{\mathrm{shots}}}$
[29], so that minimizing the one-norm should decrease the measurement overhead
for precise evaluation of the energy in quantum chemistry.

Another more widespread optimization criterion is to minimize the energy with
the help of these additional parameters $\bm{\theta^{\prime}}$. In [17, 24,
21], the optimization on $\bm{\theta}$ and $\bm{\theta^{\prime}}$ is done
jointly, while it is done in an alternate fashion in [18].

The basis in which the Hamiltonian is written—namely the very form of
$\tilde{H}(\bm{\theta^{\prime}})$—generally impacts the behavior of the
optimization. For instance, in [18], greater correlation energy was observed
for molecular systems such as $LiH$ with shallow empirical ansätze when
carrying out orbital optimization along VQE. Along with boosted
expressiveness, this could be due to a favorable effect of the variational
unitary dressing onto the optimization landscape [24].

An extreme example of the off-chip vs on-chip tradeoff and of the off-loading
of the optimization difficulty to classical vs quantum processors is supplied
by the implementation of the Hartree-Fock method presented in [11]. In this
work, the Hartree-Fock method, which consists in finding $R$, or equivalently
$\mathcal{U}(\bm{\theta^{\prime}})$, that minimizes the energy (i.e here
$\mathcal{W}(\bm{\theta})=I$), is implemented entirely ”on-chip”: states
$\mathcal{U}(\bm{\theta^{\prime}})|0\dots 0,1\dots 1\rangle$ are prepared on a
superconducting processor with so-called Givens rotation and their energy
measured. By contrast, the usual implementation of the Hartree-Fock method can
be regarded as performing the orbital rotation ”off-chip”: in this case, the
”quantum” part, namely the state preparation, is trivial as it consists in
producing states $|0\dots 0,1\dots 1\rangle$. This extreme example illustrates
the fact that if the Hamiltonian is adequately (classically) preprocessed, the
depth required to reach its ground state may be reduced.

One major drawback of the above methods is that they incur a rather large cost
in terms of classical computation: the parameter space is enlarged, and
gradient-based minimizations will require more circuit samples than the mere
circuit parameter update does. Additionally, more complex landscapes arising
from this added variational freedom have been reported on for chemical systems
studied with the unitary coupled-cluster (UCC) ansatz [16].

This motivates our choice to use a physically-motivated orbital change instead
of optimizing a parameterized orbital change. Here, our target orbital change
will be the so-called natural orbitals (NO).

###  I.3 Natural orbitals

Natural orbitals [30] are defined as the single-particle orbitals that
diagonalize the one-particle reduced density matrix (1-RDM) $D$, defined by
its matrix elements

$D_{pq}\equiv\langle\psi|c^{\dagger}_{p}c_{q}|\psi\rangle,$ (13)

of a given $M$-qubit fermionic quantum state $\ket{\psi}$. Namely, if $D$, a
positive semidefinite matrix, is diagonalized by the unitary transformation
$V$, $D=V\mathrm{diag}(n_{1}\dots n_{M})V^{\dagger}$, then the
creation/annihilation operators

$\tilde{c}^{\dagger}_{i}=\sum_{p}V^{\dagger}_{ip}c^{\dagger}_{p},\;\;\;\tilde{c}_{i}=\sum_{p}V_{pi}c_{p},$
(14)

create/annihilate electrons in the $i$-th natural orbital. The diagonal
elements $n_{1}\dots n_{M}$ are called the natural orbital occupation numbers
(NOONs) and are such that $0\leq n_{k}\leq 1$, 0 meaning an empty orbital and
1 an occupied orbital.

The natural-orbital basis is colloquially said [30] to be the basis that
minimizes the number of Fock states (or Slater determinants) in the
decomposition of a given state (or equivalently minimizes the number of
nonzero coefficients in Eq. (2)). A more rigorous statement (proven in
Appendix A) is that the NO basis minimizes the correlation entropy

$S_{\mathrm{corr}}=-\sum_{p}D_{pp}\log{D_{pp}}$ (15)

where $D_{pp}$ is the occupancy of the $p$-th single-particle orbital. Thus,
put informally, the NO basis has the least partially filled orbitals, which
explains the original colloquial statement (since states with no partial
occupations are single Slater determinant, aka single-reference states, states
with the least partial occupations should be representable by the fewest
Slater determinants).

This property has prompted the widespread use of NOs in classical
computational methods. In quantum chemistry, configuration interaction (CI)
expansions based on natural orbitals were found to range amongst the fastest
and most accurate (see e.g. [30, 31]). NOs have also been used to tackle
Anderson impurity models [32, 33, 34, 35]. There, they provide a way to
alleviate the effect of the harsh bath truncations that the exponential
scaling of the Hilbert space with the number of bath sites calls for.

We claim that the properties of the natural orbitals can also be beneficial
for the preparation of quantum states with quantum processors: intuitively,
states that can be represented as a linear combination of few Slater
determinants will require less complex circuits than states requiring more
Slater determinants. In a way, the NO basis should be the basis such that
single-particle orbital rotations (with generators quadratic in the creation
and annihilation operators, such as in Eq. (9)) have been stripped away from
the circuit and put inside the representation of the Hamiltonian, so that the
quantum processor only executes ”nontrivial” unitary operations.

Of course, we have overlooked one major hurdle so far: computing the NO
transformation requires the knowledge of the 1-RDM $D$ and thus of the state
$\ket{\psi}$ it is defined upon… Yet, this state is also the ground state we
are looking for.

In the next section, we will present an iterative method to solve this
conundrum.

Before this, let us mention recent work [36] implementing the kind of orbital
optimization we mentioned in subsection I.2 and that a posteriori backs up our
choice: in this work, using heuristic ladder ansätze, the orbital optimization
scheme was observed to converge towards orbitals very close to the NOs for
most of the chemical compounds investigated ($LiH$, $H_{2}$, $H_{2}O$, $HF$,
$NH_{3}$).

##  II Methods

Figure 2: _Flow diagram of the natural-orbitalizing VQE methods_. The
algorithm receives a Hamiltonian $\hat{H}$, defined by the tensors $h_{pq}$
and $h_{pqrs}$, whose ground state is to be determined. Repeating interleaved
sequences of (possibly adaptive) VQE runs and single-particle basis updates
until convergence, the algorithm outputs a final variational approximation
$\ket{\psi^{(K)}}=\mathcal{W}^{(K)}(\bm{\theta}^{*})\ket{\bar{0}}$ to the
target ground state, expressed in the rotated basis associated with the
transformed Hamiltonian $\hat{H}^{(K)}=(h^{(K)}_{pq},h^{(K)}_{pqrs})$.

As sketched in the previous section, the natural orbital basis we are aiming
to rotate our representation into requires the prior knowledge of the very
ground state we are looking for. In this section, we describe how we solve
this problem. Essentially, we gather information about the true NO basis by
computing the NO basis corresponding to the current converged VQE state, in
the hope that as the basis gets rotated to a basis closer to the true NO
basis, the quality of VQE, and thus of the converged state, improves.

In other words, we propose to interlace VQE runs—usual ”circuit learning”
steps—with rotations onto the natural-orbitals basis of the current
variational state—adding a representation learning aspect, as illustrated on
Fig. 2. This amounts to pushing the QPU-CPU hybridization further. We dub this
method ’natural-orbitalization’ and abbreviate it as NOization [26].

###  II.1 Representation learning through natural-orbitalization (NOization)

In the first variant of the algorithm, introduced by some of us to solve
Anderson impurity models [26], a variational circuit is picked once and for
all, as in usual VQE. A first VQE optimization over the circuit parameters in
the original basis selected to write down the Hamiltonian initially
(typically, in a condensed-matter context, the site-spin basis) is run in
order to provide a first approximation to the ground state. The 1-RDM of the
attained state is then diagonalized to provide a new set of orbitals to work
with. The VQE procedure is then repeated, and the orbitals get re-updated
according to the new converged state. The number of orbital rotations is a
hyperparameter of the method we note $K$. We will see that a few steps
($K\simeq 3$) are enough to drastically increase the expressiveness of the
circuit.

Let us now spell out how to perform the orbital rotation. At step $k\in 0\dots
K-1$ of the method, one aims at preparing the ground state of the current
Hamiltonian

$H^{(k)}=\sum\limits_{pq}h^{(k)}_{pq}c^{\dagger}_{p}c_{q}+\frac{1}{2}\sum\limits_{pqrs}h^{(k)}_{pqrs}c^{\dagger}_{p}c^{\dagger}_{q}c_{r}c_{s}.$
(16)

The tensors $h^{(k)}_{pq}$ and $h^{(k)}_{pqrs}$ are initialized e.g. as the
usual $h_{pq}$ and $h_{pqrs}$ tensors of the Hamiltonian in the site-spin
basis. They are the objects which will be rewritten to reflect the single-
particle orbital rotation.

The key step of the algorithm goes as follows. Provided the optimal VQE state
$\ket{\psi(\bm{\theta}^{*(k)})}$ returned by VQE at NOization step $k$ of the
procedure—namely minimizing the expectation value of $H^{(k)}$—we use the
quantum computer to compute the matrix elements of the 1-RDM

$D^{(k)}_{pq}=\langle\psi(\bm{\theta}^{*(k)})|c^{\dagger}_{p}c_{q}|\psi(\bm{\theta}^{*(k)})\rangle.$
(17)

In the presence of noise, the VQE optimization yields a density matrix
$\rho(\bm{\theta}^{*(k)})$ instead of a pure state
$\ket{\psi(\bm{\theta}^{*(k)})}$ and, accordingly, we actually measure
expectation values that read

$D^{(k)}_{pq}=\mathrm{Tr}(\rho(\bm{\theta}^{*(k)})c^{\dagger}_{p}c_{q}).$ (18)

Note that this quantity can actually be reconstructed at no additional cost by
aggregating the corresponding terms measured to compute the energy of the
variational state $|\psi(\bm{\theta}^{*(k)})\rangle$ (resp.
$\rho(\bm{\theta}^{*(k)})$). Then, we classically compute the transformation
$V^{(k)}$ that diagonalizes $D^{(k)}$:

$D^{(k)}_{pq}=\sum\limits_{\alpha}V^{(k)}_{p\alpha}n_{\alpha}V^{(k)\dagger}_{\alpha
q}.$ (19)

This matrix is used to update the orbital basis as:

$c_{\alpha}\rightarrow\sum_{q}V^{(k)}_{q\alpha}c_{q}$ (20)

which effectively corresponds to the following transformation being applied
onto the fermionic Hamiltonian’s coefficients:

$\displaystyle
h^{(k+1)}_{pq}=\sum\limits_{p^{\prime}q^{\prime}}V^{(k)}_{p^{\prime}p}h^{(k)}_{p^{\prime}q^{\prime}}(V^{(k)\dagger})_{qq^{\prime}}$
(21a) $\displaystyle
h^{(k+1)}_{pqrs}=\sum\limits_{p^{\prime}q^{\prime}r^{\prime}s^{\prime}}V^{(k)}_{p^{\prime}p}V^{(k)}_{q^{\prime}q}h^{(k)}_{p^{\prime}q^{\prime}r^{\prime}s^{\prime}}(V^{(k)\dagger})_{rr^{\prime}}(V^{(k)\dagger})_{ss^{\prime}}.$
(21b)

The procedure is repeated $K$ times until a convergence criterion is met or an
update budget is reached. As will be evidenced in the Results section, the
effect is to boost the expressiveness of the quantum circuit.

The strategy described above is similar in spirit to that implemented in the
so-called Perm-VQE method [37], in which qubit ordering (as opposed to the
full single-particle orbital basis) is updated iteratively, so that leading
correlations are displayed by neighboring qubits.

###  II.2 Representation and circuit learning: natural-orbitalizing adaptive
VQE (NOA-VQE) scheme

The previous method works with a fixed circuit structure, and aims at
enhancing the expressiveness of that circuit.

Having a fixed ansatz imposes some limitations. For instance, in the
noninteracting limit, $h_{pqrs}=0$, we know that the ground state is a single
Fock state in the NO basis. Therefore, in the NO basis, a circuit with only
single-qubit rotations will capture the solution. But in order to reach this
NO basis via NOization, we need to start from the original basis, in which the
ground state is not a single Fock state… and therefore requires entangling
circuits. We thus need to allow for an adaptive ansatz to capture this limit.
In fact, this strategy is probably useful in other limits: not only should the
ansatz adapt to the current single-orbital basis, but also to the noise level.

We propose to resort to ADAPT-VQE [25] (in its ”QUBIT” version [38]) to
construct the circuit at step $k$ of the NOization procedure. This variant
simply consists in using ADAPT-VQE in the VQE box of Fig. 2 instead of plain
VQE. By allowing the circuit to be constructed on the fly, one can use the
NOization procedure to reduce the depth of the variational circuit as the
orbitals get rotated.

The specifics of our implementation are the following:

(1) We pick, as a set of generators $P$ to construct our gates $e^{i\theta P}$
from, the particle-number-preserving two-qubit operators
$X_{i}X_{j}+Y_{i}Y_{j}$, the non particle-number-preserving two-qubit
operators $X_{i}Y_{j}$, the density-density-like operators $Z_{i}Z_{j}$, and
the single-qubit generators $X_{i},Y_{i}$ and $Z_{i}$. All qubit pairs $(i,j)$
are allowed.

(2) We get a reference state $|\psi_{\mathrm{ref}}\rangle$ by running VQE on a
product circuit $\otimes RY(\theta_{i},q_{i})$ (rotation of qubit $q_{i}$ with
an angle $\theta_{i}$ around the $y$ axis), as in [39]. We do this rather than
starting from some product state $\otimes X$ because the latter has vanishing
gradients for the operators of the pool considered here.

(3) We run qubit-ADAPT VQE with our previously-defined pool, with a budget of
10 circuit growth steps. At step $n$, the operator $e^{i\theta_{n}P_{n}}$
yielding the highest-magnitude gradient at $\theta_{n}=0$ over the variational
state

$|\psi(\bm{\theta_{n-1}^{*}},0)\rangle=\prod_{m=1}^{n-1}P_{m}(\theta_{m}^{*})|\psi_{\mathrm{ref}}\rangle$
(22)

is added to the circuit. The magnitude around $\theta_{k}=0$ of the gradient
associated with generator $P_{k}$ can be shown [25] to be expressible as

$G_{k}=\left|\langle\psi(\bm{\theta_{n-1}^{*}},0)|[H,P_{k}]|\psi(\bm{\theta_{n-1}^{*}},0)\rangle\right|.$
(23)

All of the parameters of the circuit are then reoptimized within VQE with the
COBYLA optimizer, except for the parameters defining the reference state.

(4) We perform orbital rotation with NOization on the last converged
variational state.

(5) We repeat steps 2 to 4 several times, either until a budget of orbital
rotations has been exhausted or until no further improvement is brought about
by the orbital updates. Here we repeated the adaptive procedure five times.

In both schemes, illustrated on Fig. 2, the orbital rotation at step $k$ is
determined by the 1-RDM $D^{(k)}$ of the current approximate ground state
$\ket{\psi^{(k)}(\vec{\theta^{*}})}$: it is the basis change which
diagonalizes this matrix. The VQE runs consist in ’circuit learning’ in that
they either simply optimize the variational parameters of a fixed ansatz, or
determine both the ansatz circuit and its optimal parameters. The orbital
rotation steps can be understood as ’representation learning’ steps.

We note that Ref. [20] uses a similar combination of adaptive VQE with orbital
optimization, where the latter is however guided by the two-particle reduced
density matrix (2-RDM).

In the next section, we present results with and without noise for the Hubbard
model with 2 and 4 sites, with onsite interactions $U=0$ and $U=1$, namely in
an uncorrelated regime ($U=0$), where the ground state is a Slater
determinant, but also in a correlated regime ($U=1$). We first document the
expressiveness enhancement allowed by the NOization procedure when using a
fixed circuit structure.

Figure 3: _Ansatz circuits_ , illustrated in the 4-qubit case (Hubbard dimer,
$N=2$).

##  III Results

Figure 4: Converged VQE energies for the half-filled Hubbard model at each
step of the natural-orbitalizing scheme for the three ansatz circuits of Fig.
3. $N=2$ (left panel) and $N=4$ (right panel) sites, with interaction $U=0$
(top row) and $U=1$ (bottom row). Noiseless and noisy (depolarizing noise)
simulations.

In this section, we apply the NOization method, in its fixed-ansatz and
adaptive versions, to the Hubbard model, defined by the Hamiltonian

$H=-t\sum\limits_{\langle
i,j\rangle,\sigma}(c^{\dagger}_{i\sigma}c_{j\sigma}+\mathrm{h.c.})+U\sum\limits_{i}n_{i\uparrow}n_{i\downarrow}-\mu\sum\limits_{i}(n_{i\uparrow}+n_{i\downarrow})$
(24)

where $t$ and $U$ denote respectively the nearest-neighbor ($\langle
i,j\rangle$) tunneling amplitude and on-site interaction of fermions of spin
$\sigma=\uparrow,\downarrow$ on a lattice with sites $i=1\dots N$. We set the
chemical potential $\mu=U/2$ to enforce half-filling of the ground state and
fix the reference energy scale through $t=1$. Values of $U$ must thus be
understood as values of the ratio $U/t$. We will investigate the non-
interacting regime $U=0$ and a correlated regime, $U=1$. The number of sites
considered will be $N=2$ (Hubbard dimer) as well as $N=4$ in a square lattice
geometry (Hubbard plaquette).

To transform $H$ into a sum of Pauli terms, we use the Jordan-Wigner
transformation [40]:

$\displaystyle
c^{\dagger}_{p}\rightarrow\hat{Z}_{1}\otimes\hat{Z}_{2}\otimes...\otimes\hat{Z}_{p-1}\otimes\frac{1}{2}(\hat{X}_{p}-i\hat{Y}_{p}),$
(25a) $\displaystyle
c_{p}\rightarrow\hat{Z}_{1}\otimes\hat{Z}_{2}\otimes...\otimes\hat{Z}_{p-1}\otimes\frac{1}{2}(\hat{X}_{p}+i\hat{Y}_{p}),$
(25b)

which maps the $M=2N$ spin-orbitals $p=(i,\sigma)$ onto qubits.

Whenever considering the effect of decoherence, we add a depolarizing noise
channel after each one-qubit gate:

$\mathcal{E}^{(1)}_{p_{1}}(\rho)=(1-p_{1})\rho+\frac{p_{1}}{3}\left(X\rho
X+Y\rho Y+Z\rho Z\right).$ (26)

After each two-qubit gate, likewise, we add a two-qubit channel
$\mathcal{E}^{(2)}_{p_{2}}(\rho)=(\mathcal{E}_{p_{2}}^{(1)}\otimes\mathcal{E}_{p_{2}}^{(1)})(\rho)$.
We adjust the depolarizing probabilities $p_{1}$ and $p_{2}$ of the single and
two-qubit depolarizing channels to match the error rates measured in
randomized benchmarking experiments on current NISQ processors, using the
relation [41] $p=1-\left(1+\frac{1}{d}\right)\epsilon_{\mathrm{RB}}$, with
$d=2^{n}$ the dimension of the subspace that is acted on by the channel. For
one and two-qubit depolarizing noise, this gives
$p_{1}=\frac{3}{2}\epsilon_{\mathrm{RB}}^{(1)}$ and
$p_{2}=1-\sqrt{1-\frac{5}{4}\epsilon_{\mathrm{RB}}^{(2)}}$.

The specific RB values we choose are those referenced by Google for the
Sycamore chip as of 2019 [42]: $\epsilon_{\mathrm{RB}}^{(1)}=0.16\%$ and
$\epsilon_{\mathrm{RB}}^{(2)}=0.6\%$.

Unless stated otherwise, VQE runs are executed with a limit of 1000 COBYLA
optimization steps. We consider several random initializations to avoid local
minima (10 random initializations for subsection III.1, 5 for subsection
III.2).

###  III.1 Expressiveness enhancement via natural orbitalization

We investigate the NOization method for three fixed ansätze documented on Fig.
3:

(i) “product” consists in a product of $RY$ parametrized rotations. It
produces only product (factorized) states, and is adapted only to non-
interacting Hamiltonians expressed in their diagonal basis;

(ii) “fSim” consists in $X$ gates on half of the qubits to fix the number of
excitations, followed by a routine made of excitation-number preserving
Google’s fSim gates [43, 11] (with one layer here), see also [44];

(iii) “LDCA” is the low-depth circuit ansatz [45] (with one cycle here).
Despite its name, it contains a relatively large number of gates, which gives
it a high expressiveness, but makes it very sensitive to decoherence.

####  III.1.1 Results for 2 sites

Figure 4(a) illustrates the behavior of the NOization procedure for the three
above ansätze, different correlation regimes and noise levels.

For $U=0$ and in the absence of noise, all circuits, including the “product”
ansatz, are able to recover the exact energy after a few steps (in the NO
basis, for $U=0$, the GS is a product state). Only the LDCA circuit is capable
of reaching the target state directly in the original basis due to its very
high expressiveness. It thus does not benefit from the NOization procedure.
When noise is included, the NOization procedure for the product circuit is
slightly impacted. Looking at the results for the fSim circuit, which is also
very shallow but exhibits noisy results similar to that of perfect results, it
is probably an artefact of the optimization and should not be observed as runs
are repeated. For LDCA, we observe that the converged energy is greatly
degraded compared with noise-free performances. This is attributable to a
large gate count which increases its sensitivity to hardware noise.
Nevertheless, NOization interestingly has the effect of slightly improving the
LDCA energy in the noisy setting.

Figure 5: Evolution of the spin-orbitals along the NOization procedure for the
Hubbard dimer ($N=2$), $U=1$, product ansatz circuit.

For $U=1$, the “product” ansatz fails to reach the exact GS energy because the
true GS is entangled due to interactions. However, the energy gets lower as
one rotates to a NO basis. For illustration purposes we display on Fig. 5 the
evolution of the spin-orbitals along the NOization procedure in this case.
Originally, we start in the site-spin basis, which means we consider the
single-fermion modes to be spanned by one Wannier orbital centered around each
site per spin species. The natural-orbital basis associated with a two-fermion
system is the so-called bonding-antibonding basis, which consists in bonding
modes $\frac{1}{\sqrt{2}}(c^{\dagger}_{\sigma 1}+c^{\dagger}_{\sigma 2})$ and
antibonding modes $\frac{1}{\sqrt{2}}(c^{\dagger}_{\sigma
1}-c^{\dagger}_{\sigma 2})$. We observe that the procedure yields orbitals
very close to this basis set, with a few spurious but small peaks. This means
that the energy bias is only attributable to the lack of expressiveness of the
circuit even in NOs. Note that this corresponds to the very particular case of
a two-fermion system, where the set of ground state natural orbitals is
independent of $U$. Similarly the fSim circuit does not reach exactly the GS
energy but gets very close as the orbitals get rotated. LDCA retains its
expressiveness and displays a behaviour very similar to the $U=0$ case.

The key point to be observed is that the NOization procedure is always found
to yield improved performances: the converged VQE energy at NOization step
$k+1$ is systematically observed to be lower than or equal to the energy at
previous step $k$. The improvement goes as high as $75\%$ improvement for the
product circuit at $U=0$. This substantiates the claim that the NOization
procedure enables to break an expressiveness wall due to unadapted
representation. It is remarkable that states which are greatly off-the-point
such as the first converged VQE states for the product and fSim circuits
provide orbital updates leading in fine to the VQE energy they yield in the NO
basis associated with the target ground state (materialized with brown lines
on Fig. 4(a)).

Raw VQE data is reported in Appendix B, in Fig. 1(a) and Fig. 1(b).

####  III.1.2 Results for 4 sites

We now turn to the Hubbard plaquette ($N=4$) in order to make sure the good
performances of the NOization procedure are not attributable to the special
character of the two-fermion case. We now have four fermions and circuits with
eight qubits.

The results are summarized in Fig. 4 (b). We can essentially make the same
observations as for 2 sites. There is however a caveat: for the LDCA ansatz,
convergence was not achieved after 1000 steps (see Fig. 2(a) and 2(b)); hence,
probably, the nonmonotonic shape of LDCA curves in 4 (b).

Raw VQE results without and with noise are given in Fig. 2(a) and 2(b) of
Appendix B.

###  III.2 Trade-off between circuit sampling and gate noise

Figure 6: _Gate noise vs shot noise trade-off_. Fixed-ansatz natural-
orbitalization strategy with the fSim circuit (dashed lines) compared to
standard VQE with the LDCA ansatz (solid lines), $U=1$. _Left_ : total energy
estimation error (bias and variance) as a function of the noise ratio $r$ for
different total shot budgets $N_{\mathrm{shots}}$, averaged over 10
realizations. _Right_ : distribution of the Pauli coefficients $|\lambda_{i}|$
for the various NOization steps at $r=0.01$. Inset: one-norm as a function of
the NOization step.

In the previous sections, we focused on remedying the problem of quantum gates
currently being considerably faulty and put forth an orbital-enhancing scheme
based on natural orbitals. As mentioned earlier, the price to pay for reduced
requirements over circuit depth that come with this scheme is that the number
of terms in the Hamiltonian may increase. As a consequence, given a fixed
number of shots for each energy evaluation, the statistical variance (or ”shot
noise”) in the energy measurement may increase because each individual term
may receive a smallest share of the shot budget. Alternatively, if one allots
a number of shots proportional to the relative strength of each term [29], the
statistical error is bounded by the one-norm of the Hamiltonian. Changing the
orbital basis likely changes this one-norm, thereby possibly increasing shot
noise.

In this section, we carry out a study on the specific case of the half-filled
Hubbard dimer at $U/t=1$ to investigate the trade-off between gate noise and
shot noise. We compare the effect of shot noise on standard VQE (with the LDCA
ansatz) and on NOization strategy with the fSim ansatz (Fig. 3(b)). Note that
unlike the LDCA ansatz which does provide a means to reach the ground state
energy $E_{0}=-2.56$ in the absence of noise, this ansatz was not found to
reach energies lower than $E_{\mathrm{min}}^{\mathrm{(fSim)}}=-2.50$ even in
noise-free simulation of VQE in natural orbitals. This means that below a
certain level of noise, the LDCA ansatz will always provide better
performances than the NOized fSim ansatz. We try here to determine what is
this level of noise, and how it changes with the number of shots.

The depolarizing noise levels we consider are chosen relatively to current
NISQ noise levels. As mentioned in the introduction to this section, we set
the maximum noise level to correspond to that of Google’s Sycamore chip [42].
We investigate fractions $r$ of this noise level in order to take into account
conceivable future hardware improvements, namely we consider depolarizing QPUs
with

$\epsilon_{\mathrm{RB}}^{(1)}(r)=r\times\epsilon_{\mathrm{RB}}^{(1)},\;\;\;\epsilon_{\mathrm{RB}}^{(2)}(r)=r\times\epsilon_{\mathrm{RB}}^{(2)}.$
(27)

To observe the effect of shot noise, we consider different budgets for the
total number of shots $n^{\mathrm{tot}}_{\mathrm{shots}}$ allocated to the
task of estimating the ground state energy. At the most atomistic level this
total budget is turned into a budget $n_{\mathrm{shots}}$ to evaluate the
expectation value over the current variational circuit. These
$n_{\mathrm{shots}}$ shots are allocated proportionally to the absolute value
of the term’s coefficient, $n_{\mathrm{shots}}^{(i)}\propto|\lambda_{i}|$, a
strategy that yields an upper bound
$\left\|\tilde{H}(\bm{\theta^{\prime}})\right\|_{1}/\sqrt{n_{\mathrm{shots}}}$
[46, 47] to the standard error on the mean. For a single VQE procedure, we
allow for up to $n_{\mathrm{iter}}=1000$ optimization steps (here, with the
COBYLA minimizer). When NOizing, the VQE procedure is repeated a number $K$ of
times. Also, to account for sensitivity of the VQE result to the initial
parameters we repeat each VQE step a number $n_{\mathrm{repeats}}=5$ of times
with random initialization, and keep the lowest-energy result. All in all,
each energy measurement of in the NOizing scheme is given a budget

$n_{\mathrm{shots}}^{\mathrm{(NOizing)}}=\frac{n^{\mathrm{tot}}_{\mathrm{shots}}}{n_{\mathrm{iter}}\times
n_{\mathrm{repeats}}\times K}$ (28)

whereas in the direct LDCA scheme we allocate a number

$n_{\mathrm{shots}}^{\mathrm{(direct)}}=\frac{n^{\mathrm{tot}}_{\mathrm{shots}}}{n_{\mathrm{iter}}\times
n_{\mathrm{repeats}}}.$ (29)

Note that in the NOization procedure, compared to the straightforward one, we
have additional 1-RDM measurements. We neglect their shot budget consumption
(here, it represents $16\times n_{\mathrm{shots}}^{\mathrm{(NOizing)}}$
additional shots per NOization step whereas a single VQE run consumes
$n_{\mathrm{iter}}\times n_{\mathrm{shots}}^{\mathrm{(NOizing)}}$ with
$n_{\mathrm{iter}}=1000$).

The final energy error stems from a bias and a variance term

$\displaystyle\mathrm{err}$
$\displaystyle=\left\langle\left|\frac{E-E_{0}}{E_{0}}\right|\right\rangle,$
(30) $\displaystyle\mathrm{Var}(E)$
$\displaystyle=\left\langle\left(\frac{E-\langle
E\rangle}{E_{0}}\right)^{2}\right\rangle,$ (31)

with $E_{0}$ the exact ground state energy. The average here corresponds to
repeating the whole VQE+NOization procedure 10 times. We will consider the
mixed figure of merit $\sqrt{\mathrm{err}^{2}+\mathrm{Var}(E)}$.

The results are summarized on Fig. 6. On the left panel, first of all, we
observe that both the standard VQE (LDCA) and the NOization (fSim) approaches
are impacted by shot noise in a similar fashion. As for the effect of
depolarizing noise, we observe a roughly linear increase in
$\sqrt{\mathrm{err}^{2}+\mathrm{Var}(E)}$ with $r$ for the standard VQE (LDCA)
approach, as opposed to virtually no impact on the NOization (fSim) method.

To explain why shot noise does not impact NOization substantially more than it
does VQE, we study the distribution of the Pauli weights for different NO
steps in the right panel of Fig. 6. In the original (site-spin) basis, the
number of terms is very small (6 terms) with substantial weight for all terms.
As NOization proceeds, the number of terms increases (reaching a hundred
terms). Yet, the distribution of the coefficients remains concentrated around
a few sizable contributions, with a lot of small coefficients. Because we
allocate the number of shots proportionally to the weight, rotating to the NO
basis thus does not substantially change the shot noise. This can be checked
more quantitatively by looking at the evolution of the one-norm (see inset),
which only doubles when going from the original basis to the NO basis.

Another striking observation is that for the smallest shot count ($5\cdot
10^{7}$ shots, corresponding to $10,000$ shots per energy evaluation),
NOization seems to outperform standard VQE for all decoherence levels.

###  III.3 Circuit shortening with the natural orbitalizating adaptive VQE
(NOA-VQE) approach

Figure 7: Adaptive natural orbitalization for the Hubbard plaquette ($N=4$).
_Upper panels:_ $U=0$. _Lower panels:_ $U=1$. _Left:_ noiseless simulations.
_Right_ : noisy simulations (depolarizing noise).

We now consider the adaptive version of the algorithm. Results are presented
for the more challenging four-site Hubbard model. We still consider the $U=0$
and $U=1$ cases, with and without depolarizing noise. Unlike the previous two
sections, we do not decide on a circuit but rather construct it on the fly by
selecting operators from a pool of our choice based on a gradient criterion,
as detailed in Section II.2. The pool corresponds to a QUBIT-ADAPT
implementation [38] with operators selected for their relevance regarding the
physics of correlated fermionic states. We consider a small budget of 10
operator additions so as to focus the study on the effect of orbital updates.

Results are presented on Fig. 7. We observe a systematic improvement in terms
of energy from one run to another, both in terms of the energy of the initial,
reference state (as wall of $RY$ gates) as well as in terms of the final
energy after the final ten-operator circuit is optimized. In all settings, we
observe that the first convergence in terms of NO steps is very quick: at the
second step, already, the method reaches energies very close to the exact
ground state energy. In the $U=0$, noiseless case, the convergence to the
exact ground state energy is perfect. Going to a correlated case ($U=1$), we
get a converged energy that is close, although not exactly equal to the exact
ground state energy. We ascribe this small deviation to the quite simple pool
we picked, the lack of a smart initialization, and the small number of ADAPT-
VQE steps we chose.

In the presence of depolarizing noise, results are similar to the noiseless
case, namely an improvement with NOization, with the important difference that
adding more operators in the adaptive ansatz is seen to degrade the converged
energy in some cases. This is because longer circuits, although possibly more
expressive, are more sensitive to noise.

##  IV Conclusion

In this work, we proposed a method to reduce the required gate count for the
task of variational state preparation of fermionic ground states.

More precisely, we put forth an orbital update technique that allows to tailor
the single-particle basis to the target state and thus incorporate part of the
correlations of the target state at the level of the Hamiltonian. This
contrasts with the usual paradigm in which one relies on the variational
circuit to generate these correlations, as well as with recent proposals which
carry out full-fledged orbital optimization. The method consists in
interleaving VQE runs with rotation of the single-particle modes to the
natural-orbital basis of the current optimized VQE state.

It was tested both within perfect state simulation (assuming a perfect QPU)
and in a noisy setting, with a depolarizing model representative of recent
devices. The systems to which the method was applied are instances of the
Hubbard model with two and four sites, without ($U/t=0$) and in the presence
of ($U/t=1$) correlations.

Within usual VQE—that is, with a fixed circuit structure—we evidence increased
expressiveness in the sense that the optimized energy gets lowered as the
orbitals get updated in both cases. Within an adaptive scheme, where the
circuit is constructed on the fly, we observe that reduced depths are
necessary to go below a given energy threshold as the orbitals get updated.
The method is relevant in the presence of noise as it provides a way to resort
to small variational circuits that are less sensitive to noise. This contrasts
with deep circuits that perform well without noise such as the low-depth
circuit ansatz or the unitary coupled cluster ansatz, as these are not
compatible with current noise levels.

One main drawback of the method is that orbital rotations do not preserve the
number of terms of the Hamiltonian or its one-norm. As a consequence, the
Hamiltonian observable to be measured onto circuits has more terms. This means
that a careful analysis of resources, both in terms of qubit quality and in
terms of shot budget, must be performed in order to assess the relevance of
our method. We show that at current and improved noise levels, for the ansätze
we studied, NOization is still preferable over standard VQE.

Another drawback of our method concerns asymmetric models such as impurity
models in which part of the qubits represent correlated impurities whereas the
others represent uncorrelated ’bath’ modes. In this context, if the different
kinds of modes do mix—as is allowed in the absence of further restrictions to
the method—then the correlation term does not concern the few impurity modes
anymore, but is being spread onto all of the modes by orbital rotation. This
prevents leveraging the different kind of physics different parts of the qubit
register must display. A way to avoid this would be to restrict the modes
rotations to the bath qubits, as is done in classical settings, for instance
in [48].

Finally, one could think about other kinds of orbital updates. For instance,
the basis that minimizes the seniority of the state (that is, the number of
unpaired electrons) is associated with a lower informational content (the
entropy of the probability distribution of the state across the different
computational basis states) [49]. From a physics point of view, the set of
natural orbitals is only concerned with single-particle degrees of freedom, as
it is defined from the 1-RDM. It would be interesting to identify a set of
spin-orbitals which are tailored to a target state according to, also, 2-RDM
information.

###### Acknowledgements.

We acknowledge useful discussions with Gaurav Gyawali. This work is part of
HQI initiative (www.hqi.fr) and is supported by France 2030 under the French
National Research Agency award number “ANR-22-PNCQ-0002”. The quantum circuit
emulations were performed on the Eviden Qaptiva platform.

##  Appendix A Proof of the correlation-entropy minimization property of
natural orbitals

Let’s introduce the so-called correlation entropy, defined for any matrix $A$
with non-negative diagonal entries as

$S_{\mathrm{corr}}(A)\equiv-\sum_{i}A_{ii}\log{A_{ii}}$ (32)

From definition (13), such a quantity can be defined for $D$ since diagonal
elements are occupation numbers associated to the single-particle orbitals
defined by the set of creation operators $\\{c^{\dagger}_{i}\\}$. It provides
a measure of the spread of the distribution onto the different orbitals of the
$\mathrm{Tr}(D)\equiv n_{\mathrm{part}}$ particles $\ket{\psi}$ describes the
state of.

We aim at showing that the transformation to NO

$D=UnU^{\dagger}$ (33)

yields a minimal correlation entropy 1-RDM (and thus, a maximally localized
occupation numbers distribution in a sense that will be made clearer later):

$U=\operatorname*{arg\,min}_{V\in U(M)}(S_{\mathrm{corr}}(V^{\dagger}DV))$
(34)

or, in other words,

$S_{\mathrm{corr}}(D)\geq S_{\mathrm{corr}}(n)$ (35)

Since $n$ is diagonal we have that $S_{\mathrm{corr}}(n)=-\mathrm{Tr}(n\log
n)\equiv S(n)$ where $S$ denotes the von-Neumann entropy. Since the von-
Neumann entropy is invariant under unitary transformation, we also have
$S(n)=S(D)$. Conversely, $S_{\mathrm{corr}}(D)=S(\mathrm{diag}(D))$ with
$\mathrm{diag}(D)$ referring to the diagonal matrix obtained from $D$ by
setting its off-diagonal elements to 0. $S(D)$ and $S(\mathrm{diag}(D))$ are
related by the so-called relative entropy of coherence [50]
$\mathcal{C}_{\mathrm{RE}}(D)\equiv
S(\mathrm{diag}(D))-S(D)=S_{\mathrm{corr}}(D)-S_{\mathrm{corr}}(n)$. From Ref.
[50] we can rewrite $\mathcal{C}_{\mathrm{RE}}(D^{1})$ as the minimum of the
relative entropy $S(D||\delta)\equiv\mathrm{Tr}(D\log
D)-\mathrm{Tr}(D\log\delta)$ with regards to incoherent (diagonal) matrices
$\delta$. Since this quantity is known to be always positive (see eg [51]), we
have the result.

##  Appendix B Raw VQE traces

((a))

((b))

Figure A.1: Raw VQE optimization traces at each step of the natural-
orbitalizing scheme for the three ansatz circuits of Fig. 3. The method is
applied here to the Hubbard model with $N=2$ sites, at half-filling and with
Coulomb interaction $U=0$ and $U=1$. We consider both noise-free state
preparation assuming a perfect QPU (upper panel 1(a)) and noisy preparation
with a depolarizing QPU model (lower panel 1(b)) with noise rates
representative of current devices.

((a))

((b))

Figure A.2: Same as Fig. A.1, with $N=4$ sites.

We report on Fig. A.1 and Fig. A.2 raw data for VQE results corresponding the
different VQE optimizations along the NOization procedure for the Hubbard
model with respectively two and four sites, at $U=0$ and $U=1$, with and
without depolarizing noise.

## References

  * Ayral  _et al._ [2023] T. Ayral, P. Besserve, D. Lacroix, and E. A. Ruiz Guzman, The European Physical Journal A 59, 10.1140/epja/s10050-023-01141-1 (2023).
  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018), arXiv: 1801.00862.
  * Louvet  _et al._ [2023] T. Louvet, T. Ayral, and X. Waintal, (2023), arXiv:2306.02620v1 .
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014), number: 1 Publisher: Nature Publishing Group.
  * Cerezo  _et al._ [2020] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, arXiv:2012.09265 [quant-ph, stat] (2020), arXiv: 2012.09265.
  * Tilly  _et al._ [2021] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, arXiv:2111.05176 [quant-ph] (2021), arXiv: 2111.05176.
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 4812 (2018), arXiv: 1803.11173.
  * Larocca  _et al._ [2024] M. Larocca, S. Thanasilp, S. Wang, K. Sharma, J. Biamonte, P. J. Coles, L. Cincio, J. R. McClean, Z. Holmes, and M. Cerezo,  , 1 (2024), arXiv:2405.00781 .
  * Ragone  _et al._ [2023] M. Ragone, B. N. Bakalov, F. Sauvage, A. F. Kemper, C. O. Marrero, M. Larocca, and M. Cerezo,  , 33 (2023), arXiv:2309.09342 .
  * Wecker  _et al._ [2015a] D. Wecker, M. B. Hastings, N. Wiebe, B. K. Clark, C. Nayak, and M. Troyer, Physical Review A 92, 062318 (2015a), arXiv: 1506.05135.
  * Arute  _et al._ [2020] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. J. Huggins, L. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Y. Niu, T. E. O’Brien, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, M. Szalay, T. Y. Takeshita, A. Vainsencher, T. White, N. Wiebe, Z. J. Yao, P. Yeh, and A. Zalcman, arXiv:2004.04174 [physics, physics:quant-ph] (2020), arXiv: 2004.04174.
  * Dallaire-Demers  _et al._ [2019] P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik, Quantum Science and Technology 4, 045005 (2019), publisher: IOP Publishing.
  * Huggins  _et al._ [2021] W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, npj Quantum Information 7, 23 (2021), arXiv:1907.13117 .
  * Yen and Izmaylov [2021] T.-c. Yen and A. F. Izmaylov, PRX Quantum 2, 040320 (2021), arXiv:2007.01234 .
  * Koridon  _et al._ [2021] E. Koridon, S. Yalouz, B. Senjean, F. Buda, T. E. O’Brien, and L. Visscher, arXiv:2103.14753 [quant-ph] 10.1103/PhysRevResearch.3.033127 (2021), arXiv: 2103.14753.
  * Sokolov  _et al._ [2020] I. O. Sokolov, P. K. Barkoutsos, P. J. Ollitrault, D. Greenberg, J. Rice, M. Pistoia, and I. Tavernelli, The Journal of Chemical Physics 152, 124107 (2020), arXiv: 1911.10864.
  * Mizukami  _et al._ [2020] W. Mizukami, K. Mitarai, Y. O. Nakagawa, T. Yamamoto, T. Yan, and Y.-y. Ohnishi, Physical Review Research 2, 033421 (2020), arXiv: 1910.11526.
  * Ratini  _et al._ [2022] L. Ratini, C. Capecci, F. Benfenati, and L. Guidoni, Journal of Chemical Theory and Computation 18, 899 (2022), publisher: American Chemical Society.
  * Bierman  _et al._ [2022] J. Bierman, Y. Li, and J. Lu,  (2022), arXiv:2208.14431 .
  * Fitzpatrick  _et al._ [2022] A. Fitzpatrick, A. Nyk, N. W. Talarico, A. Lunghi, and S. Maniscalco, , 1 (2022), arXiv:2212.11405v1 .
  * Robin and Savage [2023] C. E. P. Robin and M. J. Savage,  (2023), arXiv:2301.05976 [nucl-th, physics:quant-ph].
  * Zgid and Nooijen [2008] D. Zgid and M. Nooijen, Journal of Chemical Physics 128, 10.1063/1.2883981 (2008).
  * Ghosh  _et al._ [2008] D. Ghosh, J. Hachmann, T. Yanai, and G. K.-L. Chan, The Journal of Chemical Physics 128, 10.1063/1.2883976 (2008), arXiv:0712.2475 .
  * Moreno  _et al._ [2023] J. R. Moreno, J. Cohn, D. Sels, and M. Motta,  _Enhancing the Expressivity of Variational Neural, and Hardware-Efficient Quantum States Through Orbital Rotations_, Tech. Rep. arXiv:2302.11588 (arXiv, 2023) arXiv:2302.11588 [cond-mat, physics:quant-ph].
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature Communications 10, 3007 (2019), arXiv: 1812.11173.
  * Besserve and Ayral [2022] P. Besserve and T. Ayral, Physical Review B 105, 115108 (2022), publisher: American Physical Society.
  * Grimsley  _et al._ [2022] H. R. Grimsley, G. S. Barron, E. Barnes, S. E. Economou, and N. J. Mayhall,  _ADAPT-VQE is insensitive to rough parameter landscapes and barren plateaus_, Tech. Rep. arXiv:2204.07179 (arXiv, 2022) arXiv:2204.07179 [physics, physics:quant-ph].
  * Thouless [1960] D. J. Thouless, Nuclear Physics 21, 225 (1960).
  * Rubin  _et al._ [2018a] N. C. Rubin, R. Babbush, and J. McClean, New Journal of Physics 20, 053020 (2018a), arXiv:1801.03524 .
  * Löwdin [1955] P.-O. Löwdin, Physical Review 97, 1474 (1955).
  * Davidson [1972] E. R. Davidson, in  _Advances in Quantum Chemistry_, Vol. 6 (Elsevier, 1972) pp. 235–266.
  * Lin and Demkov [2013] C. Lin and A. A. Demkov, Physical Review B 88, 035123 (2013), arXiv: 1307.4982.
  * Lu and Haverkort [2017] Y. Lu and M. W. Haverkort, The European Physical Journal Special Topics 226, 2549 (2017).
  * Lu  _et al._ [2019] Y. Lu, X. Cao, P. Hansmann, and M. W. Haverkort, Physical Review B 100, 115134 (2019), arXiv: 1909.02757.
  * Bi  _et al._ [2019] S. Bi, L. Huang, and N.-H. Tong, Computer Physics Communications 235, 196 (2019), arXiv: 1802.02340.
  * Ratini  _et al._ [2023] L. Ratini, C. Capecci, and L. Guidoni,  (2023), arXiv:2308.08056 .
  * Tkachenko  _et al._ [2021] N. V. Tkachenko, J. Sud, Y. Zhang, S. Tretiak, P. M. Anisimov, A. T. Arrasmith, P. J. Coles, L. Cincio, and P. A. Dub, PRX Quantum 2, 020337 (2021), publisher: American Physical Society.
  * Tang  _et al._ [2021] H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021).
  * Gyawali and Lawler [2022] G. Gyawali and M. J. Lawler, Physical Review A 105, 012413 (2022), arXiv:2109.12126 .
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Z. Physik 47, 631 (1928).
  * Magesan  _et al._ [2012] E. Magesan, J. M. Gambetta, and J. Emerson, Physical Review A 85, 042311 (2012).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019), number: 7779 Publisher: Nature Publishing Group.
  * Foxen  _et al._ [2020] B. Foxen, C. Neill, A. Dunsworth, P. Roushan, B. Chiaro, A. Megrant, J. Kelly, Z. Chen, K. Satzinger, R. Barends, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, S. Boixo, D. Buell, B. Burkett, Y. Chen, R. Collins, E. Farhi, A. Fowler, C. Gidney, M. Giustina, R. Graff, M. Harrigan, T. Huang, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, P. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, E. Lucero, J. Mcclean, M. Mcewen, X. Mi, M. Mohseni, J. Y. Mutus, O. Naaman, M. Neeley, M. Niu, A. Petukhov, C. Quintana, N. Rubin, D. Sank, V. Smelyanskiy, A. Vainsencher, T. C. White, Z. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Physical Review Letters 125, 120504 (2020), arXiv:2001.08343 .
  * Cade  _et al._ [2019] C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, arXiv:1912.06007 [quant-ph] (2019), arXiv: 1912.06007.
  * Dallaire-Demers  _et al._ [2018] P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik,  (2018), arXiv:1801.01053 .
  * Wecker  _et al._ [2015b] D. Wecker, M. B. Hastings, and M. Troyer, Physical Review A 92, 042303 (2015b), arXiv: 1507.08969.
  * Rubin  _et al._ [2018b] N. C. Rubin, R. Babbush, and J. McClean, New Journal of Physics 20, 053020 (2018b).
  * Debertolis  _et al._ [2021] M. Debertolis, S. Florens, and I. Snyman, Physical Review B 103, 235166 (2021), publisher: American Physical Society.
  * Alcoba  _et al._ [2014] D. R. Alcoba, A. Torre, L. Lain, G. E. Massaccesi, and O. B. Oña, The Journal of Chemical Physics 140, 234103 (2014).
  * Baumgratz  _et al._ [2014] T. Baumgratz, M. Cramer, and M. B. Plenio, Physical Review Letters 113, 140401 (2014), arXiv: 1311.0275.
  * Wehrl [1978] A. Wehrl, Reviews of Modern Physics 50, 221 (1978).


# K-Core Maximization through Edge Additions

Zhongxin Zhou1,2,4 and Fan Zhang1 and Xuemin Lin2,3,4 and Wenjie Zhang3&Chen
Chen2
1Guangzhou University, Guangzhou, China
2East China Normal University, Shanghai, China
3University of New South Wales, Sydney, Australia
4Zhejiang Lab, Hangzhou, China
{zzxecnu<EMAIL_ADDRESS>{lxue<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

A popular model to measure the stability of a network is $k$-core - the
maximal induced subgraph in which every vertex has at least $k$ neighbors.
Many studies maximize the number of vertices in $k$-core to improve the
stability of a network. In this paper, we study the edge $k$-core problem:
Given a graph $G$, an integer $k$ and a budget $b$, add $b$ edges to non-
adjacent vertex pairs in $G$ such that the $k$-core is maximized. We prove the
problem is NP-hard and APX-hard. A heuristic algorithm is proposed on general
graphs with effective optimization techniques. Comprehensive experiments on 9
real-life datasets demonstrate the effectiveness and the efficiency of our
proposed methods.

††Zhongxin Zhou and Fan Zhang are the joint first authors. Fan Zhang is the
corresponding author.

##  1 Introduction

Graphs are widely used to model networks, where each vertex represents a user
and each edge represents a connection between two users. The cohesive subgraph
model of $k$-core, introduced by Seidman (1983), is defined as the maximal
induced subgraph in which every vertex has at least $k$ neighbors (adjacent
vertices) in the subgraph. The $k$-core of a network corresponds to the
natural equilibrium of a user engagement model: each user incurs a cost (e.g.,
$k$) to remain engaged but receives a benefit proportional to (e.g., equal to)
the number of engaged neighbors. Since the number of vertices (the size) of
$k$-core reflects the stability of a network, it is widely adopted in the
study of network engagement (stability), e.g., Bhawalkar et al. (2015);
Malliaros and Vazirgiannis (2013); Wu et al. (2013).

To prevent network unraveling, Bhawalkar and Kleinberg et al. propose the
anchored $k$-core problem which maximizes the $k$-core by anchoring $b$
vertices Bhawalkar et al. (2015), where the degree of an anchor is considered
as infinitely large. There are a series of following work to maximize the
$k$-core, e.g., Zhang et al. (2018b, 2017a); Chitnis et al. (2013). In order
to improve network stability, another basic graph operation is edge addition,
which can also be applied to $k$-core maximization. Thus, the edge $k$-core
problem is proposed Chitnis and Talmon (2018): Given a graph $G$, an integer
$k$ and a budget $b$, add $b$ edges to non-adjacent vertex pairs in $G$ such
that the $k$-core is the largest.

Figure 1: An Example of $k$-Core and Anchoring, $k=3$

######  Example 1.

Figure 1 depicts a social group $G$ with 9 users and their connections. The
willingness of a user to keep engaged is influenced by the number of her
friends (neighbors) in this group. According to the $k$-core model, suppose
$k=3$, $v_{4},v_{6}$, and $v_{7}$ firstly drop out. Their departure leads to
the leave of $v_{5}$ and $v_{8}$, as their degrees decrease to 1 which is less
than $k$. To improve network stability, we can anchor $v_{6}$ and $v_{7}$
based on anchored $k$-core model, or add an edge between $v_{6}$ and $v_{7}$
based on edge $k$-core model. Both solutions lead to a larger $k$-core induced
by the vertices in $G$ except $v_{4}$.

The edge $k$-core problem can find many applications on real-life networks:
friend recommendation in social networks, connection construction in telecom
networks, etc. For instance, in a P2P network, any user benefiting from the
network should be connected to at least $k$ other users, to exchange
resources. The holder of a P2P network can use the edge $k$-core model to find
which connections should be added between users so that a large number of
users can successfully use the P2P network Chitnis and Talmon (2018).

#### Challenges and Contributions.

In this paper, we propose a concise reduction to prove that the problem is NP-
hard and APX-hard. The only existing solution is proposed for graphs with
bounded tree-width Chitnis and Talmon (2018). However, this assumption usually
does not hold in real-life graphs, and their techniques cannot be extended to
handle general graphs. Due to the hardness of the problem, we propose a
heuristic algorithm with effective pruning techniques. The experiments are
conducted on 9 real networks to demonstrate the effectiveness and the
efficiency of the proposed methods.

##  2 Related Work

Graph processing on large data may require higher computation efficiency than
traditional queries Luo et al. (2008); Cheema et al. (2010); Luo et al.
(2011). Cohesive subgraph mining is a fundamental graph problem, with various
models such as clique Luce and Perry (1949), $k$-core Seidman (1983), $k$-fami
Zhang et al. (2018a), etc. Among the models, $k$-core is the only one known to
have a linear time algorithm Batagelj and Zaversnik (2003). The $k$-core has a
wide range of applications such as social contagion Ugander et al. (2012),
influential spreader identification Kitsak et al. (2010), collapse prediction
Morone et al. (2019), user engagement study Malliaros and Vazirgiannis (2013),
etc.

There is an efficient heuristic algorithm for the anchored $k$-core problem
Zhang et al. (2017a), while it cannot be simply applied to solve the edge
$k$-core problem. One major reason is that the anchored $k$-core model does
not change the topology of the graph while the edge $k$-core model needs to
add new edges. Besides the $k$-core maximization work introduced in Section 1,
there are some studies on $k$-core minimization under the view of against
attack Zhang et al. (2017b); Zhu et al. (2018); Medya et al. (2019).

The edge addition has been studied in different topics. Natanzon et al. (2001)
studies the hardness of edge modification problems on some classes of graphs.
Suady and Najim (2014) aims to reduce the diameter of a graph by adding edges.
Lai et al. (2005) aims to add a small number of edges in a graph to enlarge
the bandwidth. Kapron et al. (2011) aims to anonymize a given vertex set by
adding fewest edges.

##  3 Preliminaries

We consider a simple, undirected and unweighted graph $G=(V,E)$, where $V$ is
a set of vertices and and $E$ is a set of edges. We denote $n=|V|$, $m=|E|$
and assume $m>n$. Let $S=(V^{\prime},E^{\prime})$ be an induced subgraph of
$G$, where $V^{\prime}\subseteq V$ and $E^{\prime}\subseteq E$. The notations
are summarized in Table 1.

######  Definition 1.

$k$-core. Given a graph $G$, a subgraph $S$ is the $k$-core of $G$, denoted by
$C_{k}(G)$, if ($i$) $deg(u,S)\geq k$ for every vertice $u\in S$; ($ii$) $S$
is maximal, i.e., any subgraph $S^{\prime}\supset S$ is not a $k$-core.

The $k$-core of a graph $G$ can be obtained by recursively removing every
vertex $u$ and its incident edges in $G$ if $deg(u,G)<k$, with a time
complexity of $O(m)$. The $k$-cores of $G$ with different inputs of $k$
constitute a hierarchical structure of $G$, i.e., $C_{k+1}(G)\subseteq
C_{k}(G)$ for every value of $k$. The definition of $k$-shell is then derived.

######  Definition 2.

$k$-shell. Given a graph G, the $k$-shell of G, denoted by $H_{k}(G)$, is the
set of vertices in $k$-core but not in ($k$+1)-core, i.e.,
$H_{k}(G)=V(C_{k}(G)-C_{k+1}(G))$.

If we add some new edges among the vertices which are not adjacent, the
$k$-core of the graph may contain more vertices, which is named the edge
$k$-core. The added new edges are called anchors or anchor edges. In this
paper, we say anchor, add, or insert an edge interchangeably, e.g., an
inserted edge is also called an anchored edge. The edges, which may be
inserted to the graph, are called candidate anchors or candidate edges.

######  Definition 3.

edge $k$-core. Given a graph $G$ and a set of anchor edges
$A\subseteq(\tbinom{V}{2}\setminus E)$, the edge $k$-core, denoted by
$C_{k}(G+A)$, is the $k$-core of the graph $G^{\prime}=(V,E\cup A)$.

Due to the addition of anchor edges ($A$), more vertices might be retained in
$C_{k}(G+{A})$, in addition to the vertices in $C_{k}(G)$. Note that the
vertices not incident to the anchor edges may also be retained, according to
the contagious nature of $k$-core computation. The vertices following the
anchor edges $A$ to engage in $k$-core are named the followers of $A$, denoted
by $\mathcal{F}(A,G)$. Formally, $\mathcal{F}(A,G)$ is the set of vertices in
$C_{k}(G+A)\setminus C_{k}(G)$. The number of the followers reflects the
importance of the corresponding anchor edges.

#### Problem Statement.

Given a graph $G$, a degree constraint $k$ and a budget $b$, the edge $k$-core
problem aims to find a set $A$ of $b$ edges in ${V\choose 2}\setminus E$ such
that the number of followers of $A$ is maximized, i.e., $\mathcal{F}(A,G)$ is
maximized.

Notation |  Definition
---|---
$G$ |  an unweighted and undirected graph
$u$, $v$; $e$, $(u,v)$ |  a vertex in $G$; an edge in $G$
$m$; $n$ |  the number of edges in $G$; the number of vertices in $G$
$N(u,G)$ |  the set of adjacent vertices (neighbors) of $u$ in $G$
$deg(u,G)$ |  the number of adjacent vertices of $u$ in $G$
$S$ |  a subgraph of $G$
$V(S)$; $E(S)$ |  the vertex set of $S$; the edge set of $S$
$G[X]$ |  induced subgraph of the vertex set $X$ in $G$
$C_{k}(G)$; $H_{k}(G)$ |  the $k$-core of $G$; the $k$-shell of $G$
$k$; $b$ |  the degree constraint; the anchor budget
$A$ |  a set of anchor edges
$G_{A}$; $G_{e}$ |  the graph $G+A$; the graph $G+\\{e\\}$
$\mathcal{F}(A,G)$ |  the followers of the anchor set $A$ in $G$
$\mathcal{L}$; $L_{i}$ |  the onion layers of $G$; the $i$-th layer of $\mathcal{L}$
$l(u)$ |  layer index of $u$ in $\mathcal{L}$
$d^{*}(u)$ |  number of neighbors of $u$ in its higher layers
Table 1: Summary of Notations

##  4 Complexity

######  Theorem 1.

Edge $k$-core problem is NP-hard when $k\geq 3$.

###### Proof.

We reduce the edge $k$-core problem from the maximum coverage (MC) problem
Karp (1972) which is NP-hard. The MC problem is to find at most $b$ sets to
cover the largest number of elements, where $b$ is a given budget. We consider
an arbitrary instance of MC with $c$ sets $T_{1},..,T_{c}$ and $d$ elements
$\\{e_{1},..,e_{d}\\}=\cup_{1\leq i\leq c}T_{i}$. We suppose $c>k$, $c>b$, and
each of the resulting $b$ sets contains at least 2 elements, without loss of
generality. Then we construct a corresponding instance of the edge $k$-core
problem on a graph $G$. Figure 2 shows a construction example from 3 sets and
4 elements when $k=3$.

The set of vertices in $G$ consists of three parts: $M$, $N$ and $Q$. The part
$M$ contains $c$ set of vertices where each set has $d+3$ vertices, i.e.,
$M=\cup_{1\leq i\leq c}M_{i}$ where $M_{i}=\cup_{1\leq j\leq d+3}u^{i}_{j}$.
The part $N$ contains $d$ vertices. The part $Q$ is a ($k$+1)-clique where
every two vertices of the $k+1$ vertices are adjacent. For every $i$ and $j$,
if $e_{i}\in T_{j}$ in the MC instance, we add an edge between $v_{i}$ and
$u^{j}_{i}$. In Figure 2, these edges are marked in bold. For every $M_{i}$ in
$M$, we connect $u^{i}_{j}$ and $u^{i}_{j+1}$ by an edge for every
$j\in[1,d+2]$, and we also connect $u^{i}_{1}$ and $u^{i}_{d+3}$. For every
$M_{i}$, we add edges between every vertex in $M_{i}$ and the vertices in $Q$
so that every vertex in $M_{i}\setminus\\{u^{i}_{d+1},u^{i}_{d+3}\\}$ has a
degree of $k$ and every vertex in $\\{u^{i}_{d+1},u^{i}_{d+3}\\}$ has a degree
of $k-1$. We add $k-1$ edges between $v_{i}$ and the vertices in $Q$ for every
$i\in[1,d]$. The construction of $G$ is completed.

The $k$-core computation will delete all the vertices in $M$ and $N$. Thus,
the $k$-core of $G$ is $Q$. To enlarge the $k$-core, the solution (A) is to
add edges between $u^{i}_{d+1}$ and $u^{i}_{d+3}$, which is most cost-
effective. Another solution is to add edges between $M_{i}$ and $M_{j}$ where
$i\neq j$ may also enlarge the $k$-core, while this solution can always be
replaced by solution (A) with at most the same number of followers. Adding an
edge to a vertex in $N$ is not worthwhile, since we suppose each of the
resulting $b$ sets contains at least 2 elements. Thus, the edge $k$-core
problem always chooses $b$ of the $M_{i}$ in $G$ which corresponds to $b$ sets
in the MC problem. If there is a polynomial time solution for the edge
$k$-core problem, the MC problem will be solved in polynomial time. ∎

Figure 2: Complexity Reduction, $k=3$

######  Theorem 2.

For $k\geq 3$ and any $\epsilon>0$, the edge $k$-core problem cannot be
approximated in polynomial time within a ratio of $(1-1/e+\epsilon)$, unless
$P=NP$.

###### Proof.

We reduce from the MC problem using a reduction similar to that in the proof
of Theorem 1. For any $\epsilon>0$, the MC problem cannot be approximated in
polynomial time within a ratio of $(1-1/e+\epsilon)$, unless $P=NP$ Feige
(1998). Let $p$ be an arbitrarily large constant. There are two differences in
the construction of $G$: (i) $Q$ is a $p$-clique; and (ii) every $v_{i}$ is
attached by a loop of $p$ vertices where each vertex is connected to $Q$ by
$k-2$ edges except $v_{i}$. Let $\gamma>1-1/e$, if there is a solution with
$\gamma$-approximation on optimal follower number for the edge $k$-core
problem, there will be a $\lambda$-approximate solution on optimal element
number for MC, where $\lambda=\gamma+\frac{(\gamma-1)\times b(d+3)}{p\times
f}$ and $f$ is the number of followers of edge $k$-core problem. Thus, the
theorem is proved. The edge $k$-core problem is APX-hard. ∎

######  Theorem 3.

Let $f(A)=|\mathcal{F}(A)|$. We have that $f$ is not submodular for $k\geq 2$.

###### Proof.

For two arbitrary collapsers sets $A$ and $B$, if $f$ is submodular, it must
hold that $f(A)+f(B)\geq f(A\cup B)+f(A\cap B)$. Let $Q_{1}$ be a
($k$+1)-clique where vertices $u$ and $v$ are contained. Let $Q_{2}$ be
another ($k$+1)-clique. We create a vertex $w$ and connect it to the vertices
in $Q_{2}$ by $k-2$ edges. If $A=(u,w)$ and $B=(v,w)$, $f(A)+f(B)=0<f(A\cup
B)+f(A\cap B)=1$. ∎

##  5 Solution

Due to the NP-hardness and inapproximability of the problem, we resort to a
greedy heuristic which iteratively finds the best anchor, i.e., the edge with
the largest number of followers. The framework of the greedy algorithm is
shown in Algorithm 1. At Line 4, we compute the followers for each candidate
anchor edge in the complement graph of $G$, except the anchored edges in $A$
and the edges between the $k$-core vertices. Note that adding an edge between
two (non-adjacent) $k$-core vertices cannot enlarge the $k$-core, because all
the non-$k$-core vertices will still be deleted in $k$-core computation. After
the computation for every candidate anchor edge, the best anchor is chosen at
Line 6, and the $k$-core is updated by inserting the anchor edge. The time
complexity of Algorithm 1 is $O(b\times n^{2}\times m)$. As we do not
explicitly record the candidate edges, the space complexity of Algorithm 1 is
$O(m)$.

Algorithm 1 NaiveEKC

Input: $G$: a graph, $k$: degree constraint, $b$: budget
Output: a set $A$ of anchor edges

1: $A\leftarrow\emptyset$; $C_{k}\leftarrow C_{k}(G)$

2: for $i$ from $1$ to $b$ do

3: for each $e\in{V\choose 2}\setminus\\{E\cup A\cup{V(C_{k})\choose 2}\\}$ do

4: compute $\mathcal{F}(e,G+A)$;

5: end for

6: $e^{*}\leftarrow$ the edge with most followers;

7: $A\leftarrow A\cup e^{*}$; $C_{k}\leftarrow C_{k}(G+A)$;

8: end for

9: return $A$

###  5.1 Optimizations on Each Iteration

In this section, we introduce optimization techniques for the first iteration
of the greedy algorithm, i.e., the edge $k$-core problem with $b=1$. They can
be immediately applied to other iterations by replace the $k$-core of $G$ by
the edge $k$-core of $G+A$.

Basic Candidate Pruning

The following theorem locates the scope of valid candidate anchors where each
edge has at least one follower.

######  Theorem 4.

Given a graph $G$, if a candidate edge $e=(u,v)$ has at least one follower, we
have that $u\in C_{k-1}(G)$ and $v\in C_{k-1}(G)$, where at least one of $u$
and $v$ is in $H_{k-1}(G)$.

###### Proof.

Suppose $u\notin C_{k-1}(G)$. We have $e\in C_{k}(G_{e})$; otherwise, we have
$C_{k}(G_{e})=C_{k}(G)$, i.e., there is no follower of $e$. If we remove $e$
from $C_{k}(G_{e})$, then $deg(u,C_{k}(G_{e})\setminus\\{e\\})$ and
$deg(v,C_{k}(G_{e})\setminus\\{e\\})$ are at least $k-1$, because $(u,v)\in
C_{k}(G_{e})$ and their degrees decrease by 1. Thus,
$C_{k}(G_{e})\setminus\\{e\\}\subseteq C_{k-1}(G)$ which contradicts with
$u\notin C_{k-1}(G)$. Now we have proved that $u\in C_{k-1}(G)$ and $v\in
C_{k-1}(G)$. Suppose that $\\{u,v\\}\in C_{k}(G)$ and $C_{k}(G_{e})\setminus
C_{k}(G)=F\neq\emptyset$, we have that $C_{k}(G_{e})\setminus\\{e\\}$ belongs
to $C_{k}(G)$, which contradicts with $F\notin C_{k}(G)$. ∎

######  Example 2.

In Figure 3, when $k=3$, the 3-core $C_{3}(G)$ is induced by
$\\{v_{0},v_{1},v_{5},v_{6}\\}$, and the 2-shell
$H_{2}(G)=\\{v_{2},v_{3},v_{7},v_{8}\\}$. According to Theorem 4, we only need
to consider every new edge $(u,v)$ or $(v,u)$ with $u\in H_{2}(G)$ and $v\in
V(C_{2}(G))=V(C_{3}(G))\cup H_{2}(G)$, as candidates.

Onion Layer based Candidate Pruning

Given the ($k$-1)-core of $G$, the computation of $k$-core on $G$ recursively
deletes some vertices in the ($k$-1)-core. The vertices are deleted in batch
as their degrees are less than $k$ at a same time, like peeling an onion Zhang
et al. (2017c, a). The first layer of the onion, denoted by $L_{1}$, consists
of the vertices in ($k$-1)-core with degree less than $k$, i.e.,
$L_{1}=\\{u~{}|~{}deg(u,C_{k-1}(G))<k\\}$. The deletion of the first layer
vertices may decrease the degrees of other vertices, and produces the second
layer. Recursively, we have that $L_{i}=\\{u~{}|~{}deg(u,G_{i})<k\\}$ where
$G_{i}=C_{k-1}(G)-G[\cup_{1\leq j<i}L_{j}]$.

Let $l(u)$ denote the layer index of $u$, i.e., $u\in L_{l(u)}$. If
$l(u)<l(v)$, we say $u$ is at a lower layer of $v$, or $v$ is at a higher
layer of $u$. We use $\mathcal{L}$ to denote the union of layers, i.e.,
$\mathcal{L}=\cup_{1\leq i\leq s}L_{i}$ where $s=max(\\{l(u)~{}|~{}u\in
H_{k-1}(G)\\})$. Let $d^{*}(u)$ denote the number of adjacent neighbors of $u$
in higher layers and $k$-core, i.e., $d^{*}(u)=deg(u,G^{\prime})$ where
$G^{\prime}=C_{k-1}(G)-G[\\{\cup_{1\leq i\leq l(u)}L_{i}\\}\setminus u]$.

Benefit from $\mathcal{L}$, we propose an effective pruning technique.

######  Theorem 5.

Given a graph $G$, if a candidate edge $e=(u,v)$ has at least one follower, we
have that ($i$) $d^{*}(u)=k-1$ when $l(u)<l(v)$, ($ii$) $d^{*}(v)=k-1$ when
$l(v)<l(u)$, and ($iii$) $d^{*}(u)=d^{*}(v)=k-1$ when $l(u)=l(v)$;

###### Proof.

Let $O$ be the vertex deletion order of computing the $k$-core on the
($k$-1)-core. We have $d^{*}(u)\leq k-1$ and $d^{*}(v)\leq k-1$, since $u$ and
$v$ are deleted at their layers in $O$, respectively. When $l(u)<l(v)$, we
suppose $d^{*}(u)<k-1$. The anchoring of $(u,v)$ increases the degrees of $u$
and $v$ by 1, if we use the same order $O$ to compute the $k$-core again, then
$d^{*}(u)\leq k-1<k$, i.e., every vertex will be deleted at the same position
of the order $O$, including $u$ and $v$. It contradicts with $e=(u,v)$ has at
least one follower. Thus, case ($i$) and ($ii$) are proved. When $l(u)=l(v)$,
we suppose $d^{*}(u)<k-1$ without loss of generality. If we use the same order
$O$ to compute the $k$-core again with the anchored edge $e$, $u$ will be
deleted at the same position since $d^{*}(u)<k$. So $v$ can’t survive after
the next batch of deletions since $d^{*}(v)<k$. Finally, every vertex will be
deleted. ∎

Figure 3: Candidate Seletion, $k=3$

######  Example 3.

In Figure 3, when $k=3$, layer 1 contains $v_{3}$ and $v_{8}$, and layer 2
contains $v_{2}$ and $v_{7}$. By theorem 5, ($v_{2},v_{7}$) is not a promising
candidate edge because $d^{*}(v_{7})=1$, but ($v_{3},v_{8}$) is a proper
candidate edge.

Note that if $l(u)\leq l(v)$, in the $k$-core computation with anchor $e$, the
survive of $u$ may preserve some vertices before visiting $v$. Thus, in
Theorem 5, there is no degree requirement for $v$ to ensure that $e=(u,v)$ has
at least one follower.

Onion Layer based Follower Computation

A naive follower computation is to directly apply the $k$-core computation on
the graph with the existence of an anchor. An improved idea is to use the core
maintenance algorithms, which update the $k$-core for every $k$ with the
addition of an edge Zhang et al. (2017c). For the edge $k$-core algorithm, we
only need to update the $k$-core of a given $k$ with an anchor edge. Thus, we
adopt and refine the follower computation in the vertex-anchored $k$-core
algorithm (OLAK) which is shown to be more efficient than core maintenance for
$k$-core update with a fixed $k$ Zhang et al. (2017a).

In OLAK, anchoring a vertex means the degree of the vertex is infinitely
large. There is an observation that a vertex $u$ is the follower of the anchor
vertex $x$ if there is a path $x\leadsto u$ based on neighboring relations,
and $l(y)<l(z)$ for every two consecutive vertices $y$ and $z$ along the path.
This indicates that we do not need to consider the vertices without such paths
in the follower computation. Different from anchored $k$-core problem, edge
$k$-core problem requires us to add a new edge, where the vertex degree is
increased by exactly 1, and we have to consider two vertices rather than one.

In our algorithm, given an anchor edge $(u,v)$ with $l(u)<l(v)$, we generate
candidate followers layer-by-layer through activating the neighbors at higher
layers, starting from the layer of $u$. Here $u$ is the first activated vertex
and only the activated vertices can activate their neighbors at higher layers.
In this activation procedure, each activated vertex is assigned an upper bound
of its degree in the edge $k$-core. The degree upper bound of a vertex $u$ is
generated by counting its neighbors in higher layers and the activated
neighbors of $u$ at other layers. Once the degree upper bound of a vertex is
less than $k$, it will be deleted immediately and the upper bounds of its
neighbours will be decreased by 1. The follower computation is complete when
the layer-by-layer activation is finished. Note that the two vertices incident
to the anchor edge may also be deleted in this procedure. Once one of them is
deleted, there is no follower and the computation is returned.

Follower based Candidate Pruning

We can further prune some candidate edges by known results, when we retrieve
the followers of some candidate edges.

######  Theorem 6.

Given a candidate anchor edge $e_{1}=(u_{1},v_{1})$ and its follower set
$\mathcal{F}(e_{1})$, we have $\mathcal{F}(e)\subseteq\mathcal{F}(e_{1})$ for
every candidate edge
$e\in\\{(u,v)~{}|~{}u\in\mathcal{F}(e_{1})~{}$and$~{}v\in\mathcal{F}(e_{1})\cup
C_{k}(G)~{}$and$~{}u\neq v\\}$.

###### Proof.

Let $e=(u,v)$. If $\mathcal{F}(e)$ is not empty, then $C_{k}(G_{e})$ contains
$C_{k}(G)$, $u$, $v$ and $\mathcal{F}(e)$. Because $C_{k}(G_{e_{1}})$ contains
$u$ and $v$, i.e., more vertices may be added to $C_{k}(G_{e})$ after
anchoring $e_{1}$. So, $V(C_{k}(G_{e}))\subseteq V(C_{k}(G_{e_{1}}))$ and thus
$\mathcal{F}(e)\subseteq\mathcal{F}(e_{1})$. ∎

######  Example 4.

In Figure 3, when $k=3$, if we get that $v_{2}$ and $v_{7}$ are the followers
of edge ($v_{3},v_{8}$). According to Theorem 6, the followers of
$(v_{2},v_{7})$ are also the followers of ($v_{3},v_{8}$), thus
$(v_{2},v_{7})$ is not a promising candidate in current iteration.

Algorithm 2 EKC

Input: $G$: a graph, $k$: degree constraint, $b$: budget
Output: A set $A$ of anchor edges

1: $A\leftarrow\emptyset$;

2: for $i$ from $1$ to $b$ do

3: $N\leftarrow\\{(u,v)~{}|~{}u\in H_{k-1}(G+A),v\in C_{k-1}(G+A)\\}\setminus
E(C_{k-1}(G+A))$ (Theorem 4);

4: compute $\mathcal{L}$ and filter $N$ based on Theorem 5;

5: for each $e\in N$ do

6: $\mathcal{F}(e)\leftarrow$ FindFollower$(e,\mathcal{L})$ (Section 5.1);

7: update $N$ based on Theorem 6;

8: end for

9: $e^{*}\leftarrow$ the edge with most followers; $A\leftarrow A\cup e^{*}$;

10: update $N$ (Section 5.2);

11: end for

12: return $A$

###  5.2 Reusing Intermediate Results across Iterations

When one iteration of the greedy algorithm is completed, we get the best
anchor edge $A$ and the number of followers for every candidate edge in this
iteration. These results can be reused since some connected components in the
induced subgraph of ($k$-1)-shell may keep the same topology after anchoring
an edge. For each connected component $S$ where none of the vertices are
incident to the anchor $A$, we can record the largest number of followers of
one candidate anchor in $S$, for the later iterations. For the connected
component(s) $S$ where there is a vertex incident to the anchor $A$, it is
hard to reuse the results because the addition of edges may largely change the
vertex deletion order and the layer structure in $k$-core computation.

###  5.3 EKC Algorithm

We present the EKC algorithm in Algorithm 2, which optimizes the greedy
algorithm by adopting all the proposed techniques. At Line 3, the candidate
edge set is restricted to $N$ according to Theorem 4. In Line 4, we compute
the onion layers of $G$ and then exclude the candidates based on Theorem 5.
The follower computation of a candidate edge is conducted by exploring
$\mathcal{L}$ layer-by-layer, as introduced in Section 5.1. The set $N$ is
further filtered by Theorem 6 once a follower computation is completed. After
each iteration, we get the anchor edge with the most followers. The algorithm
terminates after $b$ iterations.

The resulting set of $b$ anchors is same to that in Algorithm 2. The
correctness is guaranteed by the correctness of optimization techniques. The
worst-case time complexity and space complexity are same to Algorithm 1.

##  6 Evaluation

#### Algorithms.

As far as we know, there is no existing algorithm for the edge $k$-core
problem on general graphs. In the experiment, we implement and evaluate 11
algorithms as shown in Table 2, where the bottom 6 algorithms produce the same
result, because they use the same greedy heuristic and candidate visiting
order.

Algorithm |  Description
---|---
Rand |  randomly chooses $b$ anchor edges (each with at least one follower) from $N_{1}=\\{(u,v)~{}|~{}u\in H_{k-1}(G),v\in C_{k-1}(G)\\}\setminus E(C_{k-1}(G))$
Degree |  chooses $b$ anchors from $N_{1}$ with largest degrees in $\mathcal{L}$
Layer |  chooses $b$ anchors from $N_{1}$ at highest onion layers in $\mathcal{L}$
AKC |  the anchored $k$-core algorithm in Zhang et al. (2017a)
Exact |  identifies the optimal solution by exhaustively searching all possible combinations of $b$ anchors, with all the proposed techniques
Naive |  computes a $k$-core on $G$ for each candidate to find best anchor in each iteration (Algorithm 1)
Baseline |  Naive \+ filtering candidates with Theorem 4
BL+O1 |  Baseline \+ filtering candidates with Theorem 5
BL+O2 |  BL+O1 \+ filtering candidates with Theorem 6
BL+OF |  BL+O2 \+ $\mathcal{L}$ based follower computation (Section 5.1)
EKC |  BL+OF \+ intermediate result reuse (Section 5.2)
Table 2: Summary of Algorithms Dataset | Vertices | Edges | $d_{avg}$ | $k_{max}$
---|---|---|---|---
Facebook | 4,039 | 88,234 | 43.69 | 115
Enron | 36,692 | 183,831 | 10.02 | 43
Brightkite | 58,228 | 214,078 | 7.35 | 52
Gowalla | 196,591 | 950,327 | 9.67 | 51
DBLP | 317,080 | 1,049,866 | 6.62 | 113
Twitter | 81,306 | 1,768,149 | 33.02 | 96
Stanford | 281,903 | 2,312,497 | 14.14 | 71
YouTube | 1,134,890 | 2,987,624 | 5.27 | 51
Flickr | 513,969 | 3,190,452 | 12.41 | 309
Table 3: Statistics of Datasets

#### Datasets.

Flickr is from http://networkrepository.com/, and the others are from
https://snap.stanford.edu/data/. Table 3 shows the statistics of the datasets.

#### Settings.

All programs are implemented in C++. All experiments are performed on Intel
Xeon 2.20GHz CPU and Linux System. We vary the parameters $k$ and $b$.

###  6.1 Effectiveness

Figure 4: Effectiveness of the Greedy Heuristic, k=20, b=5

Figure 4 reports the number of followers w.r.t. $b$ anchors. For random
approaches, we report the average number of followers for 100 independent
tests. In Figure 4, Degree and Layer failed to get any followers in more than
half of the datasets, because a large degree vertex or a vertex in higher
onion layer does not necessarily have a follower. A totally random algorithm
cannot get any followers in most settings. Thus, in Rand, we only choose the
edges from the set where each edge has at least one follower. We observe that
the majority of the candidate edges does not have any followers, while our
greedy heuristic always finds effective anchors in the experiments. In Figure
4, we notice that EKC preserves more than 100 followers in Flickr with only
$5$ anchors.

(a) Gowalla, b=3

(b) Brightkite, k=3

Figure 5: Exact vs EKC

We also compare the performance of EKC with the optimal solution from Exact.
Due to extremely high cost, we run Exact on an induced subgraph by 50 random
vertices. The results are from 100 independent settings. The values of $k$ and
$b$ are small due to the small-scale of data. Figure 5 shows that the margins
between EKC and Exact are not unacceptable, considering the non-submodular
property and the significantly better efficiency of EKC.

(a) Gowalla, k=20 (b) Brightkite, b=5

Figure 6: Anchor $k$-core vs Edge $k$-core

Figure 6 shows that EKC produces similar and almost larger numbers of
followers than AKC where the meaning of $b$ is different: the former for edges
and the latter for vertices. Although the anchoring of one vertex means the
preserving of many edges, EKC can produce similar followers with one anchor
edge. It shows that the edge $k$-core can enlarge the $k$-core with less graph
manipulations.

###  6.2 Efficiency

Figure 7: Running time on Different Datasets

Figure 7 reports the performance of three algorithms on all the datasets with
$k=20$ and $b=5$. The datasets are ordered by the number of edges. We find the
runtime among different datasets is largely influenced by the characteristics
of datasets and the onion layer structures in $k$-core computation. Baseline
cannot finish computation on 4 datasets within one week. BL+OF significantly
improve the performance by applying a series of techniques. EKC performs even
better given all the optimizations.

(a) Gowalla, b=20 (b) Brightkite, b=5

Figure 8: Candidate Pruning

In Figure 8, we report the number of candidate edges of different algorithms
by incrementally adding techniques, where the details of each algorithms is
given in Table 3. We can see the improvement brought by each proposed
technique.

(a) Brightkite, b=5 (b) Brightkite, k=20

Figure 9: Runtime with Different $k$ and $b$

Figure 9 studies the impact of $k$ and $r$. Figure 9(a) shows that the runtime
decreases with a larger input of $k$. This is because the number of candidate
edges become fewer with a larger $k$. Figure 9(b) reports the runtime with an
increasing $b$, which is proportional to the value of $b$. The margin between
BL+OF and EKC becomes smaller with a larger $k$, because the number of
connected components becomes less when $k$ increases. It is reported that EKC
largely outperforms the other algorithms under all the settings.

##  7 Conclusion

In this paper, we investigate the problem of edge $k$-core, which aims to add
a set $b$ of edges in a network such that the size of the resulting $k$-core
is maximized. We prove the problem is NP-hard and APX-hard. An efficient
algorithm, named EKC, is proposed with novel optimizations. Extensive
experiments on 9 real-life datasets are conducted to demonstrate our model is
effective and our algorithm is efficient.

## Acknowledgments

Xuemin Lin is supported by 2019DH0ZX01, 2018YFB10035
04, NSFC61232006, ARC DP180103096 and DP170101628. Wenjie Zhang is supported
by ARC DP180103096.

## References

  * Batagelj and Zaversnik [2003] Vladimir Batagelj and Matjaz Zaversnik.  An o(m) algorithm for cores decomposition of networks.  CoRR, cs.DS/0310049, 2003.
  * Bhawalkar et al. [2015] Kshipra Bhawalkar, Jon M. Kleinberg, Kevin Lewi, Tim Roughgarden, and Aneesh Sharma.  Preventing unraveling in social networks: The anchored k-core problem.  SIAM J. Discrete Math., 29(3):1452–1475, 2015.
  * Cheema et al. [2010] Muhammad Aamir Cheema, Ljiljana Brankovic, Xuemin Lin, Wenjie Zhang, and Wei Wang.  Multi-guarded safe zone: An effective technique to monitor moving circular range queries.  In ICDE, pages 189–200, 2010.
  * Chitnis and Talmon [2018] Rajesh Chitnis and Nimrod Talmon.  Can we create large k-cores by adding few edges?  In CSR, pages 78–89, 2018.
  * Chitnis et al. [2013] Rajesh Hemant Chitnis, Fedor V. Fomin, and Petr A. Golovach.  Preventing unraveling in social networks gets harder.  In AAAI, 2013.
  * Feige [1998] Uriel Feige.  A threshold of ln _n_ for approximating set cover.  J. ACM, 45(4):634–652, 1998.
  * Kapron et al. [2011] Bruce M. Kapron, Gautam Srivastava, and S. Venkatesh.  Social network anonymization via edge addition.  In ASONAM, pages 155–162, 2011.
  * Karp [1972] Richard M. Karp.  Reducibility among combinatorial problems.  In Complexity of Computer Computations, pages 85–103, 1972.
  * Kitsak et al. [2010] Maksim Kitsak, Lazaros K Gallos, Shlomo Havlin, Fredrik Liljeros, Lev Muchnik, H Eugene Stanley, and Hernán A Makse.  Identification of influential spreaders in complex networks.  Nature physics, 6(11):888–893, 2010.
  * Lai et al. [2005] Yung-Ling Lai, Chang-Sin Tian, and Ting-Chun Ko.  Edge addition number of cartesian product of paths and cycles.  Electronic Notes in Discrete Mathematics, 22:439–444, 2005.
  * Luce and Perry [1949] R Duncan Luce and Albert D Perry.  A method of matrix analysis of group structure.  Psychometrika, 14(2):95–116, 1949.
  * Luo et al. [2008] Yi Luo, Wei Wang, and Xuemin Lin.  SPARK: A keyword search engine on relational databases.  In ICDE, pages 1552–1555, 2008.
  * Luo et al. [2011] Yi Luo, Wei Wang, Xuemin Lin, Xiaofang Zhou, Jianmin Wang, and Keqiu Li.  SPARK2: top-k keyword query in relational databases.  IEEE Trans. Knowl. Data Eng., 23(12):1763–1780, 2011.
  * Malliaros and Vazirgiannis [2013] Fragkiskos D. Malliaros and Michalis Vazirgiannis.  To stay or not to stay: modeling engagement dynamics in social graphs.  In CIKM, pages 469–478, 2013.
  * Medya et al. [2019] Sourav Medya, Tiyani Ma, Arlei Silva, and Ambuj K. Singh.  K-core minimization: A game theoretic approach.  CoRR, abs/1901.02166, 2019.
  * Morone et al. [2019] Flaviano Morone, Gino Del Ferraro, and Hernán A Makse.  The k-core as a predictor of structural collapse in mutualistic ecosystems.  Nature Physics, 15(1):95, 2019.
  * Natanzon et al. [2001] Assaf Natanzon, Ron Shamir, and Roded Sharan.  Complexity classification of some edge modification problems.  Discrete Applied Mathematics, 113(1):109–128, 2001.
  * Seidman [1983] Stephen B Seidman.  Network structure and minimum degree.  Social networks, 5(3):269–287, 1983.
  * Suady and Najim [2014] Suaad AA Suady and Alaa A Najim.  On edge-addition problem.  Journal of College of Education for Pure Science, 4(1):26–36, 2014\.
  * Ugander et al. [2012] Johan Ugander, Lars Backstrom, Cameron Marlow, and Jon Kleinberg.  Structural diversity in social contagion.  PNAS, 109(16):5962–5966, 2012.
  * Wu et al. [2013] Shaomei Wu, Atish Das Sarma, Alex Fabrikant, Silvio Lattanzi, and Andrew Tomkins.  Arrival and departure dynamics in social networks.  In WSDM, pages 233–242, 2013.
  * Zhang et al. [2017a] Fan Zhang, Wenjie Zhang, Ying Zhang, Lu Qin, and Xuemin Lin.  OLAK: an efficient algorithm to prevent unraveling in social networks.  PVLDB, 10(6):649–660, 2017.
  * Zhang et al. [2017b] Fan Zhang, Ying Zhang, Lu Qin, Wenjie Zhang, and Xuemin Lin.  Finding critical users for social network engagement: The collapsed k-core problem.  In AAAI, pages 245–251, 2017.
  * Zhang et al. [2017c] Yikai Zhang, Jeffrey Xu Yu, Ying Zhang, and Lu Qin.  A fast order-based approach for core maintenance.  In ICDE, pages 337–348, 2017.
  * Zhang et al. [2018a] Fan Zhang, Long Yuan, Ying Zhang, Lu Qin, Xuemin Lin, and Alexander Zhou.  Discovering strong communities with user engagement and tie strength.  In DASFAA, pages 425–441, 2018.
  * Zhang et al. [2018b] Fan Zhang, Ying Zhang, Lu Qin, Wenjie Zhang, and Xuemin Lin.  Efficiently reinforcing social networks over user engagement and tie strength.  In ICDE, pages 557–568, 2018.
  * Zhu et al. [2018] Weijie Zhu, Chen Chen, Xiaoyang Wang, and Xuemin Lin.  K-core minimization: An edge manipulation approach.  In CIKM, pages 1667–1670, 2018.