\section{Introduction}
The accurate and efficient simulation of quantum many-body systems is a key challenge in computational physics, at the heart of fields as diverse as quantum chemistry, materials science, and quantum information. Open problems in the simulation of correlated materials include frustrated magnets \cite{Coey2010} and high-temperature superconductors \cite{Hasan2010}, for which the simulation of the interacting problem remains intractable to date despite significant recent progress. Dynamic correlation functions are central quantities in the understanding of quantum systems. These objects characterize the dynamics of quasiparticles through the (correlated) system, and are intimately tied to its resulting optical, magnetic, and transport properties. They also allow access to many static observables, notably including the total energy of the system. The single-particle Green's function, describing the propagation of a single electron or hole, is a particularly important example. This quantity, in its real-frequency representation, is directly accessible via (inverse) photoemission experiments, and is therefore a highly sought-after quantity in numerical simulations.

It is however a challenging quantity to compute, formally requiring access to the exponentially large set of eigenstates of the electron-attached and electron-removed systems. The Green's function has also taken on a more formal role within the Dynamical Mean-Field Theory~(DMFT), which at its heart requires the single-particle Green's function of a correlated impurity model. DMFT has emerged in the last couple of decades as one of the most promising approaches to describe strong correlation effects within a post density-functional theory framework \cite{Imada1998, Bednorz1986, Droghetti2017,Kovaleva2008, Weber2014}, thereby increasing the demand for methods that can accurately simulate the Green's function of correlated models.

A new frontier in the simulation of quantum systems has arisen from the advent of usable digital quantum computers. It has been argued since the days of Feynman that quantum systems should be used for the simulation of quantum problems, and while there are obvious advantages in this approach, it also requires the reformulation of algorithms for implementation on quantum devices, and for extraction of these quantities with the restricted (unitary) operation set intrinsic to quantum computation \cite{Lloyd1073, doi:10.1080/00268976.2011.552441, Cao2019, Bauer2020, mcardleQuantumComputationalChemistry2018}. While several methods such as quantum phase estimation (QPE) \cite{Cleve1998, Kitaev1996, Abrams1999} could provide exponential speed-up for the direct computation of the spectra of quantum systems, in practice, they require fault-tolerant quantum computing and is out of reach of current quantum devices. The current era of devices have been termed the noisy intermediate-scale quantum~(NISQ) computation\cite{Preskill2018}, where the emphasis for practical algorithms is on shallow gate depth and avoidance of significant numbers of auxiliary qubits to reduce the impact of uncorrected quantum gate/readout errors, decoherence, and paucity of coherent qubit arrays.

With these limitations in mind, the most applicable methods to date for current and near-term devices have relied on a hybrid quantum-classical optimization of the problem, termed the variational quantum eigensolver (VQE) \cite{VQE}. The VQE has also been extended into several methods to compute low lying excited states of systems in a state-specific fashion \cite{McClean2017, Colless2018, Santagati2018, Higgott2019, Jones2019, McArdle2019, Ollitrault2020, Tilly2020, zhangAdaptiveVariationalQuantum2021}. In these, a variational functional is found, which is sampled according to a parameterized quantum circuit, defining a wave function ansatz. The parameters of this circuit are then optimized classically in order to minimize the variational functional and find an approximation to the desired state. In general, this variational functional is the Ritz functional for the energy, and therefore minimizing the energy converges the circuit to an approximation to the ground state for a given Hamiltonian. With this method, the circuit can be parameterized to produce highly entangled quantum states (albeit with a compact number of parameters) which could not be efficiently sampled via classical approaches, such as the unitary coupled-cluster ansatz. It therefore enables accurate and efficient representations of many-body quantum states on quantum devices.

In this paper, we outline a modification of the VQE, which relies on a different variational functional that can be efficiently sampled, such that a parameterized state can be used on a quantum device to find arbitrary dynamical correlation functions of a given Hamiltonian. This is formulated directly in the frequency domain (either real or Matsubara frequency) and requires the optimization of a circuit defining the system correction vector at each desired frequency (although several optimization techniques help us avoiding to do a full computation at each frequency). This correction vector has previously been used for the optimization of frequency-domain equilibrium correlation functions for other parameterized states in electronic structure theory for `classical' computers, from neural quantum states \cite{ML_Corrvec} to tensor networks \cite{DMRG_Corrvec,DDMRG}, as well as being at the core of other electronic structure algorithms \cite{doi:10.1021/acs.jctc.8b00454, PhysRevB.101.045126,PhysRevB.85.205119}. It defines the response of the ground state to a perturbation at a given frequency, is prepared by an ansatz quantum circuit and optimized classically based on a sampled, variational cost function. The approach limits the requirement to a single ancillary qubit, and allows for relatively short circuit depths for use on near-term quantum devices.

The remainder of the paper is outlined as follows. In Sec.~\ref{sec:lit_rev} we present a review of relevant approaches for the extraction of equilibrium dynamical correlation functions on quantum devices. This is followed by a mathematical description of the proposed VQE for dynamic correlation functions (specifically the single-particle Green's function), with details of the required circuits for application on quantum devices. Furthermore, in Sec.~\ref{sec:qchem-cas} we discuss an approach to restrict the simulation of a correlated Green's function to an complete active orbital subspace, embedded in the mean-field response of low-energy unentangled orbitals. We present an approach to combining the subspace Green's functions rather than self-energies, in order to avoid the inversion of the noisy sampled active space Greens function matrix. In Sec.~\ref{sec:numerical_results} we numerically demonstrate the scheme on a quantum circuit simulator, applied to dihydrogen and lithium hydride with fully {\em ab~initio} correlated dynamics with interactions and quartic terms in the Hamiltonian. We introduce a noise model to mimic the two-qubit gate errors on near-term devices, which constitute the current bottleneck of the approach. Finally, Sec.~\ref{sec:discussion} provides a discussion on the implementation, advantages, and drawbacks of the method.

\section{Dynamical correlation functions on quantum devices} \label{sec:lit_rev}

There are a number of proposed approaches to access equilibrium dynamical correlation functions (and more specifically the equilibrium Green's function) on quantum devices, suitable both for longer-term and non-fault-tolerant near-term devices. These can be broadly categorized into whether they simulate the dynamics in the time or frequency domain.
For computing Green's function in the time domain, a number of approaches have been proposed based on Hamiltonian simulation following a Trotter-Suzuki \cite{Wecker2015, Bauer2016, Kreula2016, Kreula2016_,Keen2020}, or qubitization decomposition\cite{Low2019hamiltonian}. Another method based on a quantum generative model allows the computation of infinite temperature correlation functions \cite{PhysRevB.103.014301}.
While some of these methods have been successfully applied to impurity problems~\cite{Keen2020}, they generally require the gate depth (simulated time) to increase for a sufficiently fine energy resolution, and are therefore often out of reach of current devices.

A more noise-resilient method to compute Green's function in the time domain has been proposed in \citet{Endo2019}, which uses variational quantum simulation algorithm for the state evolution. Alternatively, imaginary time simulation of response properties has recently been demonstrated on current quantum devices \cite{PRXQuantum.2.010317}.

In the frequency domain (in which the current work is formulated), the computation of the Green's function is commonly recast as the solution to a quantum linear system problem \cite{PhysRevLett.103.150502}. These approaches generally rely on block encoding of the operator in order to represent it as a unitary, as well as amplitude amplification and phase estimation steps, and require a number of ancilla qubits and significant gate depths. They are therefore out of near-term capabilities of quantum computers. It is nonetheless worth noting that the convergence properties of these approaches depend on the condition number of the resulting operator and approaches to minimize the gate depth and apply the approach to dynamical quantities in quantum system have seen progress in recent years \cite{10.1145/3313276.3316366,doi:10.1137/16M1087072,PhysRevLett.122.060504,LinLin,Roggero2019, Kosugi2020, caiQuantumComputationMolecular2020}. Alternatively, QPE steps can be exchanged for a more NISQ-friendly VQE algorithm to calculate individual states or construct response functions \cite{DMFTIvan, Endo2019, zhuCalculatingGreenFunction2021, caiQuantumComputationMolecular2020}.

Here, we took a different approach which computes Green's function by solving the quantum linear system problem, but eliminates the requirement for the recasting of the operator into a unitary form. The approach extends the use of the variational quantum eigensolver (VQE) which leverages the quantum device in order to repeatedly sample an objective function of a parameterized correlated quantum state, and uses classical optimization strategies to iteratively refine that state. It employs a solution to the quantum linear system which requires only a single auxiliary qubit, has a gate depth which is only twice that of ground-state VQE algorithms, and is potentially more robust to noise on NISQ-era quantum devices compared with non-variational algorithms\cite{vqela}. A similar scheme has also recently been proposed in Ref.~\onlinecite{caiQuantumComputationMolecular2020}. Furthermore, the approach does not access individual states, and so exact spectral properties at a given frequency are obtained with increasing expressibility of the quantum ansatz.

\subsection{Green's functions}
\label{sec:method}

We first review the formulation of Green's functions in zero-temperature, time-independent fermionic systems, in order to formalize notation and terminology, and before describing the approach on quantum devices. We write the Hamiltonian in the grand-canonical ensemble as
\begin{equation} \label{eq:second_quant}
\mathcal{H}=\sum_{\alpha\beta} h_{\alpha\beta} c^\dagger_\alpha c_\beta
+ \frac{1}{2} \sum_{\alpha\beta\gamma\delta} g_{\alpha\beta\gamma\delta} c^\dagger_\alpha c^\dagger_\beta c_\delta c_\gamma - \mu \sum_\alpha c^\dagger_\alpha c_\alpha,
\end{equation}
where $c_\alpha^{(\dagger)}$ are fermionic single-particle annihilation (creation) operators, $h_{\alpha\beta}$ and $g_{\alpha\beta\gamma\delta}$ are one- and two-body matrix elements between the degrees of freedom, $\mu$ is the chemical potential, and where $\alpha, \beta, \gamma, \delta$ are indices for the fermionic modes. The retarded Green's function~$G^\mathrm{R}(t)$ represents the linear response of the system to the addition and removal of an electron, defined as
\begin{equation} \label{eq:green_function}
G^\mathrm{R}_{\alpha\beta}(t) = -i \Theta(t)\langle \psi_0 | c_{\alpha}(t) c_{\beta}^{\dagger}(0)+c_{\beta}^{\dagger}(0) c_{\alpha}(t)| \psi_0 \rangle,
\end{equation}
where $c_{\alpha}(t)=e^{i\mathcal{H}t} c_{\alpha} e^{-i\mathcal{H}t}$ are fermionic operators in the Heisenberg representation, $\Theta(t)$ is the Heaviside step function, and $|\psi_0\rangle$ is the ground state for the Hamiltonian~$\mathcal{H}$ at a given chemical potential.

In order to formally evaluate the Fourier transform of Eq.~\ref{eq:green_function}, a regularization factor of $\lim_{\eta\to 0^+} e^{-\eta |t|}$ is included. Since $\eta$ is infinitesimal, this doesn't affect the Green's function at any physical, finite value of $t$. The regularized Fourier transform then takes the form
\begin{align} \label{eq:fourrier_trans}
G^\mathrm{R}_{\alpha\beta}(\omega) = \int_{-\infty}^{\infty} G^\mathrm{R}_{\alpha\beta}(t) e^{i\omega t - \eta |t|} \mathrm{d} t,
\end{align}
where the regularization can be understood as a small imaginary shift of the target frequency $\omega$. The real-frequency retarded Green's function is then obtained as
\begin{equation}
\begin{split}
G^\mathrm{R}_{\alpha\beta}(\omega) =&
\langle \psi _{0} | c_{\alpha} \frac{1}{
\omega +i\eta + E_{0} -\mathcal{H}} c^{\dagger}_{\beta}|\psi _{0} \rangle \\
&+ \langle \psi _{0} |c^{\dagger}_{\beta} \frac{1}{
\omega +i\eta + \mathcal{H} -E_{0}} c_{\alpha}|\psi _{0} \rangle,
\end{split}
\end{equation}
where $E_0$ is the ground state energy of $\mathcal{H}$. It is often convenient to also consider the Green's function in imaginary time $\tau=it$,
\begin{align}
G_{\alpha\beta}( \tau ) =
&-\Theta ( \tau ) \langle \psi _{0} |c_{\alpha}( \tau ) c^{\dagger }_{\beta}| \psi _{0} \rangle
\nonumber\\
&+\Theta ( -\tau ) \langle \psi _{0} |c^{\dagger }_{\beta} c_{\alpha}( \tau )| \psi _{0} \rangle ,
\end{align}
where $c_{\alpha}( \tau ) =e^{\tau \mathcal{H}} c_{\alpha} e^{-\tau \mathcal{H}}$. The Fourier transform to the imaginary-frequency domain defines the Matsubara Green's function, as
\begin{equation}
\begin{split}
G_{\alpha\beta}(i\omega_n) =&
\langle \psi _{0} |  c_{\alpha} \frac{1}{i\omega_n + E_{0} -\mathcal{H} } c^{\dagger}_{\beta} |\psi _{0} \rangle \\
&+\langle \psi _{0} |c^{\dagger}_{\beta} \frac{1}{i\omega_n + \mathcal{H} -E_{0} }c_{\alpha} |\psi _{0} \rangle ,
\end{split}
\end{equation}
where $\omega_n = (2n+1)\pi /\beta$ are the Matsubara frequencies,  $\beta$ is the inverse temperature, and $n\in\mathbb{Z}$. Note that since $G_{\alpha\beta} (i\omega_n) = G^{*}_{\alpha\beta} (-i\omega_n)$, $n$ can be restricted to $n\geqslant 0$.
It should be noted that the Matsubara frequencies form a continuous domain in the zero-temperature limit.

Both the real- and imaginary-frequency Green's function can be written in one general definition as
\begin{equation}
\label{eq:master_GF}
\begin{split}
G_{\alpha\beta} (z) =&
\langle \psi _{0} | c_{\alpha} \frac{1}{z + E_{0} -\mathcal{H}} c^{\dagger}_{\beta} |\psi _{0} \rangle\\
&+\langle \psi _{0} | c^{\dagger}_{\beta} \frac{1}{z + \mathcal{H} - E_{0}} c_{\alpha} |\psi _{0} \rangle,
\end{split}
\end{equation}
such that $G^R(\omega)=G(\omega + i\eta)$ and the Matsubara Green's function is obtained by setting $z=i\omega_n$. From this form, we can define an operator $Q_{\pm}(z) = z \pm (\mathcal{H} -E_{0})$, with the sign depending on which term of Eq.~\ref{eq:master_GF} is being considered.
By introducing $V_i$ as an operator defining the perturbation of the system for a general dynamic correlation function (and taking the form of $c_{\alpha/\beta}^{(\dagger)}$ for the single-particle Green's function), we can write the separate terms which contribute to the Green's function in the form
\begin{equation}
\begin{aligned}
\label{eq:G_ij}
G_{i j}(z) =  \langle \psi_0 | V_i^{\dagger} \frac{1}{Q(z)} V_j | \psi_0 \rangle
=  \langle \psi_0 | V_i^{\dagger} | \chi_j(z) \rangle.
\end{aligned}
\end{equation}
The state $|\chi_j(z) \rangle$ defines the `correction vector' of the system, as
\begin{align}
\label{eq:def_chi_linear_eq}
|\chi_j(z) \rangle = \frac{1}{Q(z)} V_j| \psi_0 \rangle ,
\end{align}
and describes the linear response of the ground state to the dynamic perturbation at the chosen frequency. This correction vector is also found by many other algorithms in quantum simulation, including approaches in tensor networks and neural quantum states \cite{DMRG_Corrvec,DDMRG,ML_Corrvec}. We can further simplify the notation, by denoting the state $V_i|\psi_0\rangle $ as $|V_i\rangle$, such that once the correlation vector is obtained, the computation of $G_{ij}$ simplifies to $G_{ij}(z)=\langle V_i | \chi_j (z)\rangle$.

\subsection{Variational eigensolvers for correction vectors} \label{sec:corrvec}

In order to calculate $G_{ij}$ of Eq.~\ref{eq:G_ij}, we first map all fermionic operators to spin operators that can be directly measured on a quantum computer. Any operator can be mapped to a weighted sum of Pauli operators acting over the qubits, $\sum_{a} w_{a} P_{a}$, with $P_{a} \in \{X, Y, Z, I\}^{\otimes m}$, where $m$ is the number of qubits. There exist a number of possible mappings \cite{bravyi_fermionic_2002, Seeley2012, bravyi_tapering_2017, Setia2019, mcardleQuantumComputationalChemistry2018}; in this work we use the Jordan--Wigner mapping \cite{Jordan1928}. After this mapping, we use the VQE to find the ground state energy $E_0$ of a state represented by the parameterized quantum circuit, $U_{\mathrm{GS}}$\cite{VQE}. This circuit will prepare the ground state wave function $|\psi_0 \rangle$ from an initial state $\ket{0}$ of the quantum device as $\ket{\psi_0} = U_{\mathrm{GS}}\ket{0}$, though alternative ground state preparation algorithms could also be used \cite{QITE,Wang2019,yeteraydeniz2021benchmarking}.

We then aim to prepare a quantum circuit in order to solve for the correction vector via the linear equations $Q(z) |\chi_j(z)\rangle = |V_j\rangle$. In general, $|\chi_j(z)\rangle$ is not a normalized vector and therefore we rewrite it as $|\chi_j(z)\rangle =\gamma_j(z) |\tilde{\chi}_j(z)\rangle$, where $|\tilde{\chi}_j(z)\rangle$ is a normalized state and $\gamma_j(z) \in \mathbb{C}$ controls the true normalization of the correction vector for each frequency point ($z$) and perturbation (indexed by $j$). Our task is to find a second parameterized quantum circuit at each desired frequency point, $U_j(\theta)$, where $\theta$ represents the set of adjustable rotation angles in the chosen gates of the ansatz, such that the resulting circuit transforms the initial state $|0\rangle$ of the quantum device to $|\tilde{\chi}_j(z)\rangle$, i.e. $|\tilde{\chi}_j(z)\rangle=U_j(\theta(z))|0\rangle$. Once this is obtained, $\gamma_j(z)$ can be calculated as
\begin{align}
\label{eq:gamma}
\gamma_j(z) = \frac{\left\langle V_j|V_j \right\rangle}{\langle V_j | Q(z) U_j(\theta(z)) |0\rangle },
\end{align}
since at convergence
\begin{equation}
    |V_j\rangle = Q(z) |\chi_j \rangle = Q(z)\gamma_j(z) U_j(\theta)|0\rangle.
\end{equation}

To find $U_j(\theta)$, we follow a method similar to the variational linear equation resolution method presented in \citet{vqela}. We define a cost function, $g(\theta;z, j)$, as
\begin{align}
g(\theta;z, j) =&
\langle 0 | U^\dagger_j(\theta) Q^\dagger(z) Q(z) U_j(\theta) |0\rangle
\nonumber\\
& - \langle 0 | U^\dagger_j(\theta) Q^\dagger(z) \frac{ |V_j\rangle \langle V_j |}{\left\langle V_j|V_j \right\rangle} Q(z) U_j(\theta) |0\rangle \nonumber \\
=&
\langle 0 | U^\dagger_j(\theta) Q^\dagger(z) Q(z) U_j(\theta) |0\rangle
\nonumber\\
& - \frac{1}{\left\langle V_j|V_j \right\rangle} |\langle V_j | Q(z) U_j(\theta) |0\rangle|^2 \nonumber \\
=& \langle \psi | H'(z,j) |\psi \rangle, \label{eq:g}
\end{align}

where $|\psi \rangle = U_j(\theta) |0\rangle $ and $H'(z,j)$ is a new hermitian operator, defined as
\begin{equation}
\begin{aligned}
    H'(z,j) = Q^\dagger(z) \left(1 - \frac{|V_j\rangle \langle V_j |}{\left\langle V_j|V_j \right\rangle} \right) Q(z).
\end{aligned}
\end{equation}

This formalism is convenient, as $|\tilde{\chi}_j(z)\rangle$ is the normalized ground state of $H'(z,j)$, and thus $|\psi\rangle$ can be used as a variational approximation to it. When the cost function is minimized, the positive semi-definite projector $\left(1-\frac{|V_j\rangle \langle V_j |}{\left\langle V_j|V_j \right\rangle}\right)$ projects to the zero vector, as
\begin{align}
&\left(1 - \frac{|V_j\rangle \langle V_j |}{\left\langle V_j|V_j\right\rangle}\right) Q(z) |\tilde{\chi}_j(z)\rangle \nonumber \\
=&\,
\frac{1}{\gamma_j(z)}\left(1 - \frac{|V_j\rangle \langle V_j |}{\left\langle V_j|V_j\right\rangle}\right)  |V_j\rangle \nonumber \\
=&\, \vec{0},
\end{align}
ensuring that the minimal value of $g(\theta;z,j)=0$, corresponding to the lowest eigenvalue of $H'(z,j)$. Therefore, given a sufficiently flexible circuit parametrization, $U_j(\theta)$, minimizing $g(\theta;z,j)$ will yield the solution $\theta^*(z) =\mathrm{argmin}_{\theta} g(\theta;z,j)$, such that  $Q(z) U_j(\theta^*(z))|0\rangle \propto |V_j\rangle$, up to the constant factor $\gamma_j(z)$.
Hence, the final correction vector is found as
\begin{equation}
    |\chi_j\rangle=\gamma_j(z) Q(z) U_j(\theta^*(z))|0\rangle,
\end{equation}
with the value of $\gamma_j(z)$ at each frequency found via Eq.~\ref{eq:gamma}.

\subsection{Implementation on a quantum device} \label{sec:implementation}

\begin{figure}
    \begin{center}
    \begin{adjustbox}{width=0.8\linewidth}
        \begin{quantikz}[column sep=1mm]
            & \lstick{$|0\rangle $} & \gate{H} & \gate{R_z(\phi)} & \gate{X}   & \ctrl{1}   & \gate{X} & \ctrl{1}   & \ctrl{1}& \gate{H} & \meter{} \\
            & \lstick{$|0\rangle^{\otimes m} $} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}      &\ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}              & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}        & \gate{U_1} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}      & \gate{U_2} & \gate{P} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}    & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}
        \end{quantikz}
    \end{adjustbox}
    \end{center}
\caption{The quantum circuit used to compute overlaps of the form given in Eq.~\ref{eq:def_overlapterm}.
The upper horizontal line represents the ancillary qubit, with the lower line representing the $m$ qubits of the physical degrees of freedom.
The expectation value of the $Z$-measurement on the ancillary qubit yields
$\mathrm{Re}\left( e^{i\phi}\langle 0 |U_1^\dagger P U_2 |0\rangle  \right)$ \cite{Endo2019}.
Hence, by choosing $\phi=0, \pi/2$, both the real and imaginary part of Eq.~\ref{eq:def_overlapterm} can be sampled.}
\label{fig:overlap circuit}
\end{figure}

For a practical quantum implementation, we require an efficient way to sample from $g(\theta;z,j)$ at each desired $z$ and $V_j$, compute the norm of the correction vector $\gamma_j(z)$ from Eq.~\ref{eq:gamma}, and ultimately the overlap $\langle V_i|\chi_j(z)\rangle$ in order to compute $G_{ij}(z)$.

Recall that both $V_{i}$ and $Q(z)$ have been mapped to a weighted sum of Pauli operators.
Symmetric expectation values can therefore be computed in a straightforward fashion. These include e.g. the first term of Eq.~\ref{eq:g}, as
\begin{equation}\label{eq:expect_pauli_ops}
\begin{split}
&\langle 0|U^\dagger(\theta) Q^\dagger(z) Q(z) U(\theta)|0\rangle \\
=& \sum_{\alpha,\beta} w^*_{Q,a}w_{Q,b} \langle 0|U^\dagger(\theta) P_{a} P_{b} U(\theta)|0\rangle,
\end{split}
\end{equation}
since we can decompose $Q(z)$ into Pauli operator strings as $Q(z)=\sum_{a} w_{Q,a} P_{a}$ for some $w_{Q,a}\in \mathbb{C}$. Note that $P_{a} P_{b}$ is equal to another Pauli operator times a constant. Quantities in Eq.~\ref{eq:expect_pauli_ops} can be sampled in a straightforward manner in keeping with traditional VQE implementations \citep{VQE}. However, calculations of overlaps between different states, including the calculation of $\gamma_j(z)$, $G_{ij}(z)$, and the second contributing term to $g(\theta;z,j)$ in Eq.~\ref{eq:g}, require a more involved computation, similar to a Hadamard test circuit. All of these contributions can be decomposed into the form
\begin{equation}
\label{eq:def_overlapterm}
\left\langle 0 |U_1(\theta_1) P U_2(\theta_2) |0\right\rangle,
\end{equation}
where $U_1$ and $U_2$ are two (potentially parametrized) quantum circuits. We use the method described in \citet{Endo2019} to calculate the real and imaginary parts of Eq.~\ref{eq:def_overlapterm}, requiring the use of an additional single ancillary qubit, with the circuit shown schematically in Fig.~\ref{fig:overlap circuit}.
This circuit evolves the state as follows.
\begin{enumerate}
    \item After $H$ and $R_z(\phi)$ gates:
    \begin{equation*}
    2^{-1/2}(|0\rangle +e^{i\phi}|1\rangle ) |0\rangle
    \end{equation*}
    \item After $X$, the controlled $U_1$, and $X$ gates:
    \begin{equation*}
        2^{-1/2}(|0\rangle U_1 |0\rangle +e^{i\phi}|1\rangle |0\rangle )
    \end{equation*}
    \item After the controlled $U_2$ and $P$ gates:
    \begin{equation*}
        2^{-1/2}(|0\rangle U_1 |0\rangle +e^{i\phi}|1\rangle P U_2 |0\rangle )
    \end{equation*}
    \item Before measurement:
    \begin{equation*}\begin{split}
        2^{-1} |0\rangle \left( U_1|0\rangle  + e^{i\phi} P U_2 |0\rangle \right) &\\ + 2^{-1} |1\rangle  \left( U_1|0\rangle  - e^{i\phi} P U_2 |0\rangle  \right)
    \end{split}\end{equation*}
\end{enumerate}

Measurement of the ancillary qubit in the $Z$ basis can return $0$ or $1$, with probabilities $p(0)$ and $p(1)$ respectively. If $\phi=0$, then it can be verified that $p(0)-p(1)=\mathrm{Re} \langle 0| U_1^\dagger P U_2|0 \rangle$, while if $\phi=\pi$ it follows that $p(0)-p(1)=-\mathrm{Im}\langle 0|U_1^\dagger P U_2|0\rangle$. This approach allows us to compute all remaining quantities of interest and avoids the requirement of inversion of the $U_1$ circuit of the traditional Hadamard test.

\subsection{Active space embedding of Green's functions}
\label{sec:qchem-cas}

In order to overcome a restriction to small system sizes due to the limitations of qubit number, we can also seek a multiscale approach, whereby the quantum solution to the Green's function is solved in a strongly correlated subset of the total degrees of freedom, in the presence of a mean-field description of the response outside this. These multi-resolution or embedded approaches to the correlated dynamics are common in many approaches where the full solution is intractable, and where a subspace of the most important correlated physics can be identified, often involving subspaces chosen via locality or energetic arguments \cite{tilly2021reduced,Bauer2016,GalliQuantumEmbedding,yamazaki2018practical,dhawan2021dynamical,Takeshita2020,Yao2020,rossmannek2020quantum, Yalouz_2021}. In this work, we will demonstrate a simple approach to embed a noisy Green's function, e.g. sampled from a quantum device, within a mean-field environment described efficiently on a `classical' computer.
This will employ a `complete active space' partitioning of the molecular orbitals of the {\em ab initio} system.

In this approach the ground state wave function and response vectors are expressed as the tensor product of a highly entangled quantum state within a complete active space~(CAS) of molecular orbitals, and a simple product state in the remaining `core' degrees of freedom, as
\begin{equation}
    \ket{\psi_0} = \ket{\psi_\mathrm{CAS}} \otimes \ket{\phi_\mathrm{core}} .
\end{equation}
The core state is therefore restricted to a single Slater determinant. In this work, the CAS space is simply chosen to span the highest-energy occupied and lowest-energy unoccupied single-particle orbitals which arise from a prior mean-field calculation. This space is therefore expected to span the most important low-energy quantum fluctuations required to describe the correlated physics of the system.
As a result of this approximation, many-body correlated states only have to be calculated in the CAS (in the presence of the static mean-field Coulomb and exchange potentials from the core electrons), reducing the number of qubits which are needed to describe a system.
Within the active space, an effective Hamiltonian is constructed according to
\begin{multline}\label{eq:h_cas}
    \mathcal{H}_\mathrm{CAS}
    = \sum_{ij \subset \mathrm{CAS}} (h_{ij} + V^\mathrm{core}_{ij}) c^\dagger_i c_j \\
    + \frac{1}{2} \sum_{ijkl \subset \mathrm{CAS}} g_{ijkl} c^\dagger_i c^\dagger_j c_l c_k
    ,
\end{multline}
where $V^\mathrm{core}$ describes the Coulomb repulsion and exchange interaction with the electrons of the core space, in a partitioning which is common in the field of quantum chemistry \cite{Roos1980,olsen11}.
The CAS Green's~function can then be calculated using the correction vector method described in Sec.~\ref{sec:corrvec}, but
with the reduced dimensionality subspace Hamiltonian $\mathcal{H}_\mathrm{CAS}$ of Eq.~\eqref{eq:h_cas} used in place of~$\mathcal{H}$.

To go from the active space Green's~function to a Green's~function of the entire system, the core part has to be accounted for.
The `formally correct' approach to this involves an inversion of the resulting $G_\mathrm{CAS}$, to calculate an active space self-energy according to the Dyson equation, as
\begin{equation}\label{eq:s_cas}
    \Sigma_{\mathrm{CAS}}(z) = G_{0,\mathrm{CAS}}^{-1}(z) - G_{\mathrm{CAS}}^{-1}(z) ,
\end{equation}
where~$G_{0,\mathrm{CAS}}$ is the non-interacting active space Green's~function \cite{PhysRevB.101.045126}.
The self-energy of Eq.~\eqref{eq:s_cas} can then be added to the (generalized) Fock matrix, $F$, of the CAS $\oplus$ core system, thus defining the full system Green's~function
\begin{equation}\label{eq:g_casci_dyson}
G(z) = \frac{1}{z - F - \Sigma_\mathrm{CAS}(z)}
.
\end{equation}

However, while this is formally correct (and similar to a self-consistent step in dynamical mean-field theory \cite{Bauer2016}), there are some considerations to take into account in this context. The quantum-derived Green's function in the CAS space will contain random noise, from the sampling of the required expectation values such as those in Eqs.~\ref{eq:expect_pauli_ops} and \ref{eq:def_overlapterm}. The result is that this random error will manifest as systematic error in the self-energy due to the non-linear inversion operations required in Eqs.~\ref{eq:s_cas} and \ref{eq:g_casci_dyson} \cite{PhysRevB.98.085118}.
For this reason, we also consider a simpler, but approximate scheme. In this, the non-interacting Green's function~$G_0(z) = \left[ z - F \right]^{-1}$, and the CAS Green's function are combined directly,
i.e.
\begin{equation}\label{eq:g_casci_nondyson}
G(z) = G_0(z) + P_\mathrm{CAS} \left[ G_\mathrm{CAS}(z) - G_0(z) \right] P_\mathrm{CAS}
,
\end{equation}
where $P_\mathrm{CAS}$ is the projector onto the CAS space.
The two approaches become identical, if all of the off-diagonal elements of the generalized Fock matrix~$F$ which couple the CAS and core spaces are zero. In this case, the operations of matrix inversion and projection commute and the approaches should be identical.
While this is not generally the case, $F$ is often still block-diagonally dominant, and therefore the Green's~function of Eq.~\eqref{eq:g_casci_nondyson}
becomes a good approximation to Eq.~\eqref{eq:g_casci_dyson}. Importantly however, this avoids the necessity of matrix inversions, which would themselves introduce a systematic error when performed on stochastically sampled active space Green's functions due to the matrix inversions.
In Sec.~\ref{sec:numerical_results}, we will investigate the use of these active space approaches to Greens functions, and compare the Green's function calculated from Eq.~\ref{eq:g_casci_dyson}, which we will term `Dyson-CAS', with that of the approximation of Eq.~\ref{eq:g_casci_nondyson} (`Non-Dyson-CAS'), in order to quantify this approximation which avoids matrix inversion.

\section{Numerical results\label{sec:numerical_results}}

In this section, we present the results of quantum simulation experiments to test the quantum correction vector VQE on two different molecules, hydrogen (\ce{H_2}) and lithium hydride (\ce{LiH}). Their Hamiltonians are obtained using the STO-3G basis. For the \ce{LiH} molecule, we have further tested the embedding of an active space sampling of the correlated Greens function in a low-energy subspace, detailed in Sec.~\ref{sec:qchem-cas}, to reduce the required number of active spin-orbitals to four. The Hamiltonian is then mapped to a weighted sum of Pauli operators through the Jordan--Wigner mapping \citep{Jordan1928}.

We focus our attention on the accuracy of the single particle spectrum of the systems, for which we define $G(z)$ as the trace of the imaginary part of the full Green's function matrix on either the real or Matsubara axis, as
\begin{equation}
    G(z) = \sum_{\alpha} \mathrm{Im} \left[G_{\alpha \alpha}(z)\right].
\end{equation}
Furthermore, we restrict our simulations to conserve spin symmetry, by working in a restricted formalism where the two spin channels are constrained to be the same, precluding single-particle spin symmetry breaking.
The simulation is realized using the computing platform Hyrax developed by Rahko, with the matrix elements defining the molecular Hamiltonian of each system computed using the {\tt PySCF} simulation package \cite{pyscf1,pyscf2}.

\subsection{Noise model}

In order to model the performance of realistic quantum devices, we include the effects of noises coming from both the variance in the measurements, as well as a two-qubit depolarization noise arising from simulated decoherence of the qubits due to environmental interactions. In particular, this noise model in the two-qubit depolarizing channel is defined as
\begin{align}
    \label{eq:2qb_dep}
    \mathcal{N}(\rho) =
    (1-p_2) \rho+\frac{p_2}{15}\left(\sum_{a,b} \sigma_{a}^{(1)} \sigma_{b}^{(2)} \rho \sigma_{a}^{(1)} \sigma_{b}^{(2)}\right)-\frac{p_2}{15} \rho,
\end{align}
where $\sigma^{(1)}$ ($\sigma^{(2)}$) denotes a Pauli matrix operating on the first (second) qubit that each two-qubit gate acts on, and $a,b$ sums over all possible Pauli matrices (including the identity) \cite{nielsen_chuang_2010}. This two-qubit depolarizing noise model is applied after every two-qubit gate in the circuit. The $p_2$ parameter is set to $10^{-3}$ for our numerical experiments in Sec~\ref{sec:numerical_results}. On current generation quantum devices, the one-qubit error rate is at least one order of magnitude smaller than the two-qubit error rate and therefore we assume that the one-qubit errors are negligible in the experiments. Furthermore, due to the large number of two-qubit gates in the circuit of Fig.~\ref{fig:overlap circuit}, we anticipate two-qubit errors to be the dominant errors in the computation. We note that when three-qubit gates CCX and CCZ are used for the controlled unitary gates in Fig.~\ref{fig:overlap circuit}, they are decomposed into one and two-qubit gates as done in the work of \citet{on_cnot_cost}.

In order to mitigate the effect of this two-qubit noise, we use a simple exponential extrapolation scheme in our sampling of expectation values, as introduced in Ref.~\onlinecite{EndoExp}, and realized experimentally in Ref.~\onlinecite{DigitalZNE}. To do this, we first execute the circuit with the original noise strength of $p_2=10^{-3}$, and subsequently boost the noise in the simulation by setting $p_2$ to twice its original value. On a physical quantum device this boosting of noise in particular channels can be achieved in several different approaches \cite{Cincio2018,DigitalZNE}. The expectation values of the measurements at these different noise levels are fit to an exponential form as $E(p_2)=E(0) e^{-b p_2}$ to obtain the extrapolated noise-free estimate of $E(0)$. Finally, to converge the effect of measurement noise due to the intrinsic variance in the expectation values, we repeat measurements $10^6$ times to control the statistical accuracy.

\subsection{Simulation details}

\begin{figure}
\begin{center}
\begin{quantikz}
    & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \gategroup[wires=3,steps=4,style={inner sep=6pt}]{Repeated $d_{\mathrm{HEA}}$ times} \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \gate{S} & \ctrl{1} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \gate{S} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}  \\
    & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \gate{S} & \targ{} & \ctrl{1}  & \gate{S} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} \\
    & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \gate{S} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \targ & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{} & \gate{S} & \ifthenelse{\the\pgfmatrixcurrentcolumn>1}{\arrow[arrows]{l}}{}
\end{quantikz}
\caption{An illustration of the hardware efficient ansatz for three qubits.
Each $S$ represents a series of single qubit rotation
gates. Each $S$ gate maintains its own angle, and can be adjusted
independently. The depth of the ansatz is denoted by $d_\mathrm{HEA}$.}
\label{fig:HEA_ansatz}
\end{center}
\end{figure}

For both the ground state and correction vector circuits, we used the Hardware Efficient Ansatz (HEA) \cite{HEAAnsatz}, illustrated in Fig.~\ref{fig:HEA_ansatz}. The single qubit gates ($S$ in Fig.~\ref{fig:HEA_ansatz}) and number of repeated blocks (depth $d_\mathrm{HEA}$ in Fig.~\ref{fig:HEA_ansatz}) were dynamically adjusted to converge the results without overparameterization of the circuit. The angles, $\theta$, in these single-qubit gates were adjusted to minimize the cost function for the ground state and correction vector circuits (Eq.~\ref{eq:g}) using the Rotosolve algorithm \citep{Rotosolve}.
We note that for the optimization of the correction vector, the cost function $g(\theta;z,j)$ has a theoretical minimum of $0$, and therefore the optimization is continued until a small
threshold value of $\varepsilon > 0$ is achieved. Typically, we use $\varepsilon=0.05$, and decrease this value if the
obtained Green's function at specific frequency points deviates significantly from its values at neighbouring frequencies. Furthermore, the circuit depths for the correction vectors, $d_{\mathrm{HEA}}$ can vary for different values of $z$. If significant residual values of $g(\theta;z,j)$ remain, then the circuit depth can be increased to provide a more expressive ansatz for the correction vector.

These variations in the precision required for the optimization and flexibility of the circuit as $z$ is changed are closely tied to the condition number of the operator~$H'(z,j)$. Large condition numbers are found near the poles of the Green's function on the real-frequency axis, which is regularized by the choice of $\eta$, and therefore convergence is faster and requires shallower circuits away from these poles, or on the Matsubara frequency axis. To accelerate the optimization, we initialize the circuit angles for the correction vector optimization by reusing optimized angles from nearby frequencies for subsequent optimization points.

\subsection{Hydrogen molecule}

\begin{figure*}
    \centering
    \subfloat[Retarded Green's function at $d=2\text{\;\r{A}}$]{%
        \label{fig:h2-d2-real}%
        \includegraphics[width=0.45\linewidth]{H2d2_Real_Greens_function_ExpExtra_final.pdf}
    }%
    \hspace{8pt}%
    \subfloat[Matsubara Green's function at $d=2\text{\;\r{A}}$]{%
        \label{fig:h2-d2-imag}%
        \includegraphics[width=0.45\linewidth]{H2d2_Matsubara_Greens_function_ExpExtra_final.pdf}
    }%
    \\
    \subfloat[Retarded Green's function at $d=3\text{\;\r{A}}$]{%
        \label{fig:h2-d3-real}%
        \includegraphics[width=0.45\linewidth]{H2d3_Real_Greens_function_ExpExtra_final.pdf}
    }%
    \hspace{8pt}%
    \subfloat[Matsubara Green's function at $d=3\text{\;\r{A}}$]{%
        \label{fig:h2-d3-imag}%
        \includegraphics[width=0.45\linewidth]{H2d3_Matsubara_Greens_function_ExpExtra_final.pdf}
    }%
    \caption[Matsubara Green's function and retarded Green's function calculated for \ce{H_2} at different distance $d$]{
    Single-particle spectrum of \ce{H_2} computed at two different bond lengths, on both the real-frequency and Matsubara axis. This spectrum is computed from the trace of the imaginary part of the retarded Green's function, calculated via a quantum simulation of the VQE correction vector algorithm  (VQE Correction Vector) and compared to
        the exact solution and uncorrelated mean-field limit (Hartree--Fock).
    }
    \label{fig:H2_GF_plots}%
\end{figure*}

Figure~\ref{fig:H2_GF_plots} shows the simulated single-particle spectrum for the \ce{H_2} molecule at two different bond lengths. These spectra are computed on both the real-frequency (where a broadening of $\eta=0.05$\;Ha is used) and Matsubara axis. In both of these domains the effect of the correlated physics of the system on the single-particle spectrum is clear, with the interacting physics opening the fundamental gap of the system. We expect the correlation effects and entanglement of the orbitals to be larger at the stretched geometry of 3\;\AA, where this opening of the charge gap is most significant and where the situation is similar to the opening of a Mott gap due to interactions from a condensed matter perspective. At less stretched geometries, correlation-induced substructures are observed at this energy resolution, splitting the peaks compared to the mean-field solution at this basis size.

The quantum correction vector approach reproduces the exact spectrum to high fidelity, where the exact results are computed via an exact inversion of the full Hamiltonian on classical computers for these small system sizes. The discrepancy on the Matsubara axis is almost unobservable on the scale of the plot, and only has a maximum of $0.22\;\text{Ha}^{-1}$ over the entire domain of the Matsubara axis. Note that the analytic properties of the Green's function ensure that the same information is contained in both the Matsubara and real-frequency domains. On the real-frequency axis, the discrepancies primarily manifest as noise in the transition amplitudes at the excitation energies. All of these observations are consistent with the condition number of $H'(z,j)$ being the determining factor in the accuracy of the method, and hence why an adaptive approach to ansatz expressibility and convergence thresholds is useful in these regimes. Furthermore, we should note that formal properties such as sum rules and causality which should be present in the final Green's function are obeyed to good accuracy, and at full convergence should be obeyed exactly.

\subsection{Lithium Hydride}

\begin{figure}
    \centering
    \subfloat[$d=2$, Retarded Green's function]{%
        \label{fig:lih-d2-real}%
        \includegraphics[width=0.45\textwidth]{LiH_Real_Greens_function_ExpExtra_final_inc_Dyson.pdf}
    }%
    \hspace{8pt}%
    \subfloat[$d=2$, Matsubara Green's function]{%
        \label{fig:lih-d2-imag}%
        \includegraphics[width=0.45\textwidth]{LiH_Matsubara_Greens_function_ExpExtra_final_inc_Dyson.pdf}
    }%
    \caption[Matsubara Green's function and retarded Green's function calculated for the lithium hydride molecule at distance $d=3$ angstrom]{
        Single-particle spectral function for \ce{LiH} on both the real-frequency and Matsubara axis. Results are presented from the embedding of a two-orbital, low-energy active space, in the deep-lying core $\sigma$ orbital at the mean-field level. Results are shown from quantum simulation with the VQE correction vector `VQE Correction Vector (Non-Dyson)', exact classical approaches (both with `Dyson' embedding and `Non-Dyson' embedding), and the uncorrelated mean-field result (`Hartree--Fock'). The inset in Fig.4.a emphasizes the region of the retarded Green's function where there is a noticeable difference between the data obtained from VQE correction vector method in the two different approaches.
    }
    \label{fig:LiH}%
\end{figure}

In Figure~\ref{fig:LiH} we apply the method to the \ce{LiH} molecule, at a slightly stretched bond length of 2\;\AA.  For this system we take an active space of the highest energy occupied and lowest energy unoccupied degrees of freedom from a prior Hartree--Fock calculation as a four-qubit space in which to solve for the correlated VQE correction vector. This leaves two electrons in a low-energy $\sigma$-bonding orbital which are considered inert for the purposes of the entanglement and correlation effects with the rest of the electrons, and therefore just contribute a Coulomb and exchange static potential to the active space Hamiltonian, as described in Sec.~\ref{sec:qchem-cas}. The correlation is found to induce a satellite peak at energy ranges significantly outside the energy window of the orbitals chosen.

There are two approaches to combine the response of this inactive core orbital with the rest of the system, which were termed `Dyson' and `Non-Dyson' in Sec.~\ref{sec:qchem-cas}. The former approach combines the self-energies in the different spaces and is formally correct. However, it requires inversion of Green's function estimates from the quantum solver, which will introduce systematic error in the final result for any non-zero random error. In contrast, the `Non-Dyson' approach combines the Greens functions directly, neglecting the effect of the off-diagonal block of the generalized Fock matrix coupling the spaces in the full system mean-field spectrum. However, this does not require any inversion of the sampled Green's function, and since the initial orbitals are designed to have no coupling in the mean-field picture, it is only the presence of correlations which can subsequently induce a small coupling between the active space and core orbital. Furthermore, from a computational perspective, the `Dyson' embedding requires the calculation of the off-diagonal Green's function elements in the active space, however this does not require the optimization of any further correction vectors.

The assumption that the difference between the `Dyson' and approximate `non-Dyson` embedding is slight is corroborated by the results of Fig.~\ref{fig:LiH}, where we can see that for this system the two approaches for an exact active space Green's function embedded in the wider system are virtually indistinguishable. This is expected for this choice of active space embedding, where the correlated subspace is chosen from canonical Hartree--Fock orbitals and where the hybridization in the final state between the active space and core orbital is expected to be small. From another perspective, the differences are expected to be minor in the case where there is not significant relaxation of the full system electronic density due to the correlation effects.

We can also observe the differences between these approaches when the active space Green's function is solved with the quantum simulator and significant simulated quantum noise. Both of the quantum-simulated VQE solutions to the embedded active space Green's function again reproduce the correlated spectrum to a high degree of accuracy, with only relatively minor differences arising from the magnitude of the transition amplitudes at the excitation energies. However, close to the poles, the accuracy of the `Non-Dyson' approach is marginally superior, indicating that the bias due to the inversion of the noisy Green's function matrix is more significant than the neglect of the hybridization between the active space and core orbital for this particular choice of active space orbitals. While this is an interesting observation, this conclusion is likely to still be sensitive to the specific choice of active space orbitals, molecular system, and noise in the sampled active space Green's function. The low-energy satellite peak in Fig.~\ref{fig:LiH} at $\sim-0.45$~Ha is however underestimated slightly in both quantum approaches.

\subsection{Perspective}\label{sec:discussion}

It is important to consider the scalability and applicability of the VQE correction vector algorithm to current and near-term
quantum devices \citep{Preskill2018}.
The algorithm prepares both the ground state and the correction vector for different frequency points
using circuits representing a variational ansatz. The approach of parameterized quantum circuits has been widely used for the practical simulation of both ground and excited states in recent years \cite{mcardleQuantumComputationalChemistry2018, McClean2017, Colless2018, Santagati2018, Higgott2019, Jones2019, McArdle2019, Ollitrault2020, Tilly2020}, although the scalability of the optimization problem for the gate parameters is a source of debate in the literature \cite{PhysRevB.102.075104,wang2021noiseinduced,Grant2019initialization}. This is expected to be compounded by the condition number of the~$H'(z,j)$ operator generally being larger than that of its parent Hamiltonian used to simulate the ground state.
The correction vector~$|\chi_j(z)\rangle$ at each frequency contains contributions from all the excited states of $\mathcal{H}$ in the probed symmetry sector. However, on the real-frequency axis, this combination of states will be dominated by the excited state corresponding to the nearest excitation energy to the sampled frequency.
In this case, we expect the required ansatz complexity (depth, for example) to express the correction vector to be similar to the required ansatz complexity to describe these excited states, which will again depend on the level of entanglement between the degrees of freedom.
Therefore, we expect the scalability of our correction vector ansatz to be comparable
to other state-specific VQE-based algorithms \cite{DMFTIvan,tilly2021reduced}.

We can also consider the overall circuit complexity in this approach compared to ground-state VQE implementations. This complexity is increased due to the requirement to sample overlaps of the type described in Sec.~\ref{sec:implementation}, which require a significant number of controlled-unitary gates. This has inhibited experiments on the current generation of quantum computers, due to their relatively high two-qubit gate error rates. However, the largest circuit depths required in the correction vector VQE algorithm (in particular, required for the calculation of $g(\theta)$ in Eq.~\ref{eq:g}) is only twice that of the circuit depths of a ground state VQE algorithm. Furthermore, only a single ancillary qubit is required for the measurements of overlaps. This certainly puts it in the range of near-term devices when the bottleneck of the two-qubit gate fidelity is improved with respect to noise. Furthermore, all optimization algorithms designed to perform the classical parameter updates can be straightforwardly ported from ground state VQE codes \cite{Rotosolve}, as well as potentially more sophisticated approaches such as the AdaptVQE for ansatz selection \cite{Adaptvqe}.

\section{Conclusions}

In this paper, we propose an approach for the simulation of arbitrary dynamic correlation functions on quantum devices, solving directly in the (real or Matsubara) frequency domain via a coupled quantum-classical variational approach.
The method combines the conventional VQE and a variational linear equation solver to prepare the correction vector via a parameterized circuit ansatz in order to calculate the single-particle Green's function of correlated molecular systems.
The approach can converge to the exact Green's function, assuming a sufficient expressibility in the circuit representing the ground and linear response states, with a circuit depth which is then independent of the resolution on the frequency axis. This contrasts with state-specific approaches to constructing Green's functions which start from the spectral representation, and which suffer from a cutoff in the energy resolution \cite{DMFTIvan,Endo2019}, or time-evolution methods where the length of time the state is propagated for (which generally changes the circuit depth) controls the frequency resolution of the resulting dynamical function \cite{Bauer2016, Kreula2016, Kreula2016_, Keen2020}. Although the overall gate depths are only a factor of two longer than a ground-state VQE approach, the increase in the number of two-qubit entangling gates results in the algorithm being particularly sensitive to the levels of two-qubit gate noise.

In our numerical simulations, the method reproduced the Green's function over both real and imaginary frequencies for the \ce{H_2} and \ce{LiH} molecules to high accuracy, under an simulated noise model for errors in the dominating two-qubit gates. We focus both on simulation of the whole Green's function, and also on an approach to isolate a correlated low-energy `active' subspace in which to simulate the correlated dynamics. These dynamics can be embedded into the mean-field spectrum of the whole system either via a combination of Green's functions, or self-energies. While the latter approach is formally correct, it requires inversion of the noisy sampled Green's function matrix. We demonstrate that the Green's function approach to combining the response of the active and external spaces is sufficiently accurate for the \ce{LiH} system, where the correlation effects are not too strong, and the hybridization between these canonical degrees of freedom is small. Future work will look to expand this approach to other dynamical correlation functions of interest, as well as investigate the use of the approach within a self-consistent dynamical mean-field embedding theory \cite{Bauer2016,tilly2021reduced,Keen2020,DMFTIvan,GalliQuantumEmbedding} and assess the viability of the approach on physical quantum devices.

\begin{acknowledgments}

H.C. is supported through a Teaching Fellowship from UCL.
J.T. is supported by and industrial CASE (iCASE) studentship, funded by and UK EPSRC [EP/R513143/1], in collaboration with University College London and Rahko Ltd. G.H.B. gratefully acknowledges support from the Royal Society via a University Research Fellowship, as well as funding from the European Research Council (ERC) under the European Union\textquotesingle s Horizon 2020 research and innovation programme (Grant Agreement No. 759063). This research project is in great part funded by Innovate UK, project Quantifi from UK Research and Innovation, as part of the UK National Quantum Technologies Programme and Industrial Strategy Challenge Fund. We thank the members of the Amazon Quantum Solutions Lab at Amazon Web Services (AWS) for their valued support of project Quantifi, and provision of computational resources used in this work. The views expressed are those of the authors and do not reflect the official policy or position of AWS or the members of the Amazon Quantum Solutions Lab.

\end{acknowledgments}

\bibliographystyle{unsrt}


\section{Introduction}

With the popularity of high-speed internet and 5G mobile communication, low-cost video capture devices, and social media, a large number of videos are generated and uploaded online every second.  This brings many issues including copyright, privacy, and child protection.  In this work, we focus on partial video copy detection (PVCD), which identifies videos and localizes copy segments in a large database of long untrimmed reference videos.

The PVCD is a very challenging problem because of a large variety of editing and transformations such as rotation, scaling, cropping, captioning, watermarking, FPS change, picture in picture.  For the background, the readers are referred to \cite{VCDB2},\cite{TEMPORAL-NETWORK} and references therein.

There are three main parts to the PVCD algorithm. The first one is to identify the videos which have copy segments in them. The second one is to represent a video with a feature. The last one is to localize the copy segment in a candidate video. We explore PVCD in all these three directions and our contributions are summarized as follows.

Traditional PVCD algorithms using hand-crafted features can use features saved in an inverted file and use a fast search such as k-nearest-neighbor (KNN) approaches \cite{VCDB1, MYLiu, r2,r13}. However, since the deep CNN features are explored for PVCD, there are no works studying to use KNN approaches likely because the deep CNN features are hard to encode or compress further, and require a lot of resources in computation and storage. Instead, reference videos are scanned one by one \cite{VCDB2}.

Our first contribution is to propose a PVCD framework to use fast KNN search on a database of deep CNN features. To our knowledge, this is the first work that uses a KNN search method on deep CNN features for PVCD.  We first construct a global feature database consisting of all frame features of all reference videos. Given a query segment, we use the frame features of the query segment to search in the global database and return the top-K answers for every query frame.  These answers' scores are accumulated by their reference video ID.  The reference videos whose scores are the highest are returned as the candidate videos.

Our second contribution is to propose to use a transformer \cite{Transformer} encoder to improve the frame CNN feature. In \cite{VCDB2} and \cite{LAMV}, CNN is used to extract the feature for a video frame. This feature is for one individual frame image without considering the temporal context. Inspired by the transformer \cite{Transformer} for its power to model the context in the language model, image classification \cite{VIT}, object detection \cite{DETR}, video classification \cite{ViViT}, video action recognition \cite{Transformer-action-recognition}, video re-localization \cite{Transformer-VideoReloc}, video point cloud \cite{Fan_2021_CVPR} etc., particularly by \cite{VISIL-Ext}, we explore its application in an location-sensitive video representation for the PVCD.

Our last contribution is on the copy segment localization algorithm. We still use the temporal network algorithm \cite{TEMPORAL-NETWORK}, but we make some major modifications to improve the performance as well as the speed. Since the feature and the similarity matrix are optimized, we can tighten the parameters such that the algorithm runs a lot faster while achieving the same or better performance.

Please note that, our contributions are not about exploring brand-new ideas. Instead, we optimize existing techniques into a well-engineered framework to solve a practical problem. We evaluate the new PVCD framework on the VCDB dataset \cite{VCDB1}. The best F1 score is improved from the state-of-the-art 0.8025 to 0.8764 by a large margin. We also evaluate the speed in a few critical parts of the PVCD and show the expedition of PVCD.

\section{Related Work}

Video copy detection includes near duplicate video detection or retrieval (NDVD or NDVR), PVCD, and content-based copy detection (CCD). For an example of NDVR, please refer to \cite{kordopatis2017dml} and references therein. For CCD, please refer to the TRECVID 2011 report \cite{TRECVID2011}.

Two important parts of PVCD are the video representation and the copy localization. In early times, a global feature is generated to represent a video for its simplicity, fast speed, and low data storage \cite{r11,r3,r4}. However, since its performance is poor, we only consider PVCD using frame features. Traditional local features for PVCD include LBP, SIFT, and variants, which are then aggregated into frame features using BOW, Fisher vectors, or similar methods \cite{r2,r12,r13,TEMPORAL-NETWORK}.

Before CNN feature is used, Hamming or binary embedding and some simple local sensitivity hashing methods are used \cite{VCDB1} on these handcrafted features. Other works using compact feature embedding include \cite{MYLiu,r2,r13}. In \cite{MYLiu}, a video fingerprint-based inverted search is implemented. In \cite{r2}, a multiple level search in inverted files is  implemented on a database of a single reference video.  Due to the significant inferior performance of these hard crafted features, we only consider CNN feature in this work.

The work \cite{VCDB2} is the first one that explores deep CNN features on PVCD.  Since then, the LAMV \cite{LAMV} implements the PVCD in the frequency domain using the temporal matching kernel \cite{TMK}. Since CNN is used in PVCD and similar tasks \cite{VCDB2}, \cite{VISIL}, \cite{VISIL-Ext}, the reference videos are scanned one by one, and no fast KNN search approaches are available except for \cite{MMTA}. In \cite{MMTA}, a low-cost global descriptor in combination with a decision strategy adapted from reinforcement learning is used.

For the copy localization, early PVCD works include \cite{HOUGH-VOTE}, and \cite{TEMPORAL-NETWORK}, and the references therein. In these two works, the authors propose the Hough voting algorithm and the temporal network (or network flow) algorithm, which find the location of the copy in the similarity matrix between the query video and the reference video. These two algorithms are evaluated in \cite{VCDB2} on CNN features.

All these algorithms use the similarity matrix directly. For PVCD, this similarity matrix is ideally a diagonal sub-matrix at the location of the copy segment. In reality, however, it is not. The ViSiL \cite{VISIL} first explores the information hidden in the similarity matrix. The authors propose a simple CNN to train the noisy sub-matrix to a better diagonal sub-matrix for similar tasks such as NDVR, event video retrieval. In these tasks, they only care if there is a same event in the reference video, but do not consider the event's location. Another disadvantage is it has to scan every reference video to find the candidate video, which cannot be used for a real video retrieval problem.

In 2019, there was a PVCD challenge on Data Fountain organized by China Computer Federation \cite{BDCI-CCF}. In this challenge, the query was an entire video, instead of a video segment as in the VCDB standard protocol. In one of the forum discussions, one team used a KNN search in a global feature database, but the localization part was unclear \cite{HYData}. Part of our work is inspired by this scheme. Based on the KNN search in the global feature database, we develop a whole framework for PVCD.

\section{Proposed Method}

\begin{figure}[t]
    \centering
    \includegraphics[scale=0.7]{diagram.png}
    \caption{Diagram of our PVCD algorithm}
    \label{fig1}
\end{figure}

\subsection{System Diagram}

The diagram of our PVCD algorithm is shown in Figure \ref{fig1}. There are three main functions. The first one is the video representation. In this work, the frame is extracted in the video uniformly Then; the CNN feature is generated for the frame. This feature is for a single frame without considering the temporal context. To overcome this deficiency, we propose to use a transformer encoder, highlighted in blue in Figure \ref{fig1}, to improve the representation of the video, which we call representation learning in this paper.

The second function is the global feature database - the foundation for a fast KNN search and video scoring. The output of this function is a shortlist of candidate videos that have the closest matches to the query video. More details will be given shortly in Section 3.2.
The third function is the copy segment localization, which finds the location of the copy segment in the candidate videos. We propose to use a diagonalization CNN, highlighted in green in Figure 2, to refine the similarity matrix. After that, the temporal network algorithm is used to localize the copy segment.  The outputs are the starting and ending location of the copy segments and their similarity scores.

Please note that there are offline and online processes in this diagram, denoted with solid and dashed lines in Figure \ref{fig1}. In the offline process,  the global feature database is created from all the reference videos. This offline process does not affect the online query speed. In the online process, the query video is processed to generate frame features, then the KNN search and copy detection. Please note that this query approach does not scan videos. It first finds candidate videos from the KNN search results. Then copy detection is done on this shortlist of candidate videos.

\subsection{Global Feature Database}

In this algorithm, frames are uniformly extracted out of every reference video at one frame per second (FPS=1).  A CNN feature is generated for every frame. These features are all inserted into the global feature database. The video ID and frame index are combined to index the global feature database.

Assume there are $N$ reference videos, denoted by $r_i,i=1,2,...,N$. Let the extract frames of the video $r_i$ be $r_{i,j}$, where $j=1,2,...,M_i$, and $M_i$ is the number of extracted frames in video $r_i$. For simplicity, we also use $r_{i,j}$ to represent the CNN feature of this frame, without ambiguity. So the features of the reference video $r_i$ can be described with a matrix $R_i = [r_{i,1}; r_{i,2}; ... ; r_{i,M_i}]$, which has a dimension $M_i \times d$, where $d$ is the length of the feature.  The global feature database for all the reference videos consists of features ${r_{i,j}}$, where $i=1,2,...,N$, $j=1,2,...,M_i$ in ascending order on $j$ then on $i$. The index $i$ is mapped to a video ID if a video ID instead of number index is preferred. The index for the database is a composite index of the video ID $i$ and frame index $j$. This function is pre-processed on the reference video database and is indicated as offline in Figure 2.
Since all CNN features are $L_2$ normalized, the cosine distance, a measure of the similarity of a pair of frames, is equivalent to the inner product or Euclidean distance. These two distances are the most commonly used metric in KNN algorithms. So this global feature database is KNN-searchable with any KNN algorithm.

Similarly, assume $q_t, t=1,2,...,T$ be the frame feature of the query video, where $T$ is the number of extracted frames. The feature matrix is $Q = [q_1; q_2; ... ; q_M]$. When this query video detects a copy in the reference video database, the frames and features are generated in the same way on the query video, and the features are used to search in the global feature database.  This is indicated as online process in Figure 2.  For every frame $q_t, t=1,2,...,M$, the first top-K (denoted $top\_K\_all$) matches  $r^t_{vid(l),j(l)}$ are returned, where $l=1,2,..,K$, $vid(l)$ is the video ID of the $l-th$ match, and the $j(l)$ is the frame index of this match. In other words,
\begin{multline}
    sim(r^t_{vid(1),j(1)})>=sim(r^t_{vid(2),j(2)})>=\\
    sim(r^t_{vid(3),j(3)})>=...>=sim(r^t_{vid(K),j(K)})
\label{EQ1}
\end{multline}
where $sim()$ is the frame similarity score.

Next in the video scoring function, every returned answer is assigned to a video according to answer's video ID as follows,

\begin{equation}
sim(vid) = \sum_{t} \sum_{l} sim(r^t_{vid(l),j(l)}) \; s.t. \;  vid(l)=vid \end{equation}

The videos whose total similarity scores $sim(vid)$ are the highest are the candidate videos.  We return the top-K (denoted $top\_K\_video$) videos for every query, so the returned candidate videos $vid_i$ have,

\begin{equation}
sim(vid_1)>= sim(vid_2)>=...>=sim(vid_K)
\end{equation}
The value of $top\_K\_video$ can be chosen appropriately for a particular dataset.

\subsection{Modified Temporal Network Localization}

After the candidate videos are available, the temporal network algorithm \cite{TEMPORAL-NETWORK} is used to localize the position of the copy in the candidate videos.  Since the number of candidate videos is usually very small, the running speed is acceptable, close to real-time.

There are a few key parameters that affect the temporal network localization performance. The first one is the top-K matched (denoted $top\_K\_one$) frames for every frame in the query. Only the top K frames with the highest matching scores are used, and the rest are ignored. The second parameter is the $max\_step$, which is the allowed maximum step in both row and column directions in the similarity matrix when a network flow path is formed, which is the same as the constraint $wnd$ in \cite{TEMPORAL-NETWORK}. The third one is a similarity threshold $sim_{th}$ where similarities smaller than this threshold are not considered.

Consider the similarity matrix between a pair of videos ($q$,$r$), where $q$ is the query video segment, and $r$ is one reference video found in the previous step. Let us call every element in the matrix a 2D point $p$, where $p[0]$ and $p[1]$ are its row and column index. The similarity matrix is expressed as  $sim(p)=sim(q_{p[0]},r_{p[1]})$. Please note that for simplicity, we drop the video index of $q$ and $r$ since there is only one $q$ and one $r$ in this pair. The index $i$ is now the frame index.  The temporal network localization becomes a network flow optimization problem,

\begin{multline}
path(p),score,L = \max_p {{\sum_{i}^L{sim(p_i)}}} \; s.t. \; \\
(a): 0<p_{i+1}[0]-p_{i}[0]<=max\_step    \\
(b): 0<p_{i+1}[1]-p_{i}[1]<=max\_step    \\
(c): p_i[1] \in top-K(p_i[0])   \\
(d): sim(p_i)>= sim_{th}   \\
(e): |p_{i+1}[0]-p_{i}[0] - (p_{i+1}[1]-p_{i}[1])| < max\_diff \\
\end{multline}

The conditions (a) and (b) are in the original temporal network algorithm \cite{TEMPORAL-NETWORK}.  The conditions (c)-(e) are our modifications.  These modifications show a big impact on both speed and performance. In the original temporal network \cite{TEMPORAL-NETWORK}, the network flow runs in the whole similarity matrix. In conditions (c) and (d), we propose to only keep the reliable top-K matched reference video frames for every query frame. As a result, the similarity matrix becomes a sparse matrix, and the network flow runs more accurately and efficiently.

The condition (e) is added to force the network flow to run along with a narrow band along the diagonal direction.  Since we optimize the similarity matrix to near a diagonal matrix in the copy segment location, we can tighten the conditions by using smaller $K$ and $max\_diff$ and expedite the algorithm.

The returned result includes the longest $path$ of length $L$ whose total $score$ is maximized. From the path ${p_1,p_2,...,p_L}$, the starting and ending location of copy clip in $q$ and $r$ are defined, as well as the video similarity score,

\begin{equation}
\begin{aligned}
& q_{start} = p_1[0],\;  q_{end} = p_L[0] \\
& r_{start} = p_1[1],\;  r_{end} = p_L[L] \\
& sim(q,r) = score/L \\
\end{aligned}
\end{equation}

\subsection{CNN Feature with Transformer}

The transformer is introduced in the paper \cite{Transformer}. It is one of the breakthroughs in machine learning in recent years. It has been studied extensively on different applications, including \cite{VIT, DETR, ViViT,Transformer-action-recognition, Transformer-VideoReloc}, etc. Particularly in \cite{VISIL-Ext}, transformer is proposed to improve representation CNN feature in a video retrieval task. Therefore, we explore its application in PVCD - a video retrieval and re-localization task, where we add the supervision of location.

The self-attention mechanism of the transformer is effective at modeling long-term dependencies within a sequence input. The encoder can be used to aggregate temporal context in representing videos. Even though the encoder output keeps the same dimension as the input, the contextual information within a longer range of each frame feature is accounted for.

Let the video descriptor $X$ of dimension ${M \times d}$.  Denote the parameter matrices of the transformer encoder by $W_Q$,$W_K$, $W_V$. The video descriptor X is encoded into Query Q (please do not confuse with the query video feature matrix), Key K and Value V by three different linear transformations,

\begin{equation}
Q = X^TW_Q, K = X^TW_K, V = X^TW_V
\end{equation}
The self-attention is calculated as,

\begin{equation}
Att(Q, K, V ) = softmax ( \, \frac{QK^T}{\sqrt{d_k}} ) \,V
\end{equation}
where $\sqrt{d_k}$ is used to scale the dot product. This self-attention value is then passed to a LayerNorm layer and a Feed-Forward Layer to get the output of the Transformer encoder. The multi-head attention is usually used, where the multi-head attention outputs are concatenated before the Feed-Forward layer.
This transformer encoder can be used on different video frame features. In this work, it is applied to the Resnet-RMAC \cite{LAMV} feature. Our proposed transformer encoder for the CNN feature is shown in Figure \ref{fig2}.

\begin{figure}[t]
    \centering
    \includegraphics[scale=0.5]{transformer.png}
    \caption{CNN feature using a transformer encoder.}
    \label{fig2}
\end{figure}

An example of a similarity matrix before and after applying the transformer is shown in Figure \ref{fig3}. In this example, it can be clearly seen that after the transformer encoder, the similarity matrix becomes more like a diagonal sub-matrix at the copy segment location.

\begin{figure}[t]
    \centering
    \includegraphics[scale=1.1]{transformer-effect.png}
    \caption{An example of an original similarity matrix and its transformer encoded one.}
    \label{fig3}
\end{figure}

\subsection{Loss Function}

We use the notation defined in Section 3.1. For query video feature matrix $Q$ of dimension $T \times d$, and the reference video feature matrix $R$ (we omit index $i$ for simplicity) of dimension $M \times d$, the similarity matrix $P$ is,
\begin{equation}
P =[P_{i,j}] = QR^T
\end{equation}

Assume there is a copy segment of the query video Q in the reference video R. Let the starting location and length of the copy segment in the query video be $(s_q,L)$, in the reference video be $(s_r,L)$. Then the ground-truth (GT) similarity matrix $P_{GT}$ is,

\begin{equation}
\begin{aligned}
&P_{GT} = [P'_{i,j}]\\
&P'_{s_q+i,s_r+i} = 1,i = 0,1,...,L-1\\
&P'_{i,j} =0 \; elsewhere
\end{aligned}
\end{equation}

The loss function is defined as the mean square error (MSE) of all the elements between matrix $P$ and $P_{GT}$,

\begin{equation}
loss = \frac{1}{MT} \sum_{i,j}{[(P_{i,j}-P'{i,j})w_{i,j}]^2}
\end{equation}
where the $w_{i,j}$ is a weight factor to balance the number of ones and number of zeros in the similarity matrix, since the number of zeros is a lot more than the number of ones, we use a weight of 0.1 for zeros and 1.1 for ones, as an example.

Please note that there are some arguments that PVCD may not have an identity sub-matrix at the copy segment location due to video rate (FPS) change or reshuffling. This is partially true. However, in PVCD, we use the segment F1 as a metric and treat detection of any frames as a successful detection of the copy segment. Furthermore, since adjacent frames are usually very similar, an identity submatrix is likely showing up at some locations of the copy segment. The evaluation results on the VCDB dataset will demonstrate our algorithm's effectiveness.

\section{Experiments}

\subsection{Implementation Details}

All the videos are first uniformly extracted with 1 FPS into frames, as widely used in many works including \cite{VCDB1,VCDB2,HengYang}. We test other FPS and find 1 FPS gives a good trade-off between performance and speed. When the scale of the reference videos gets very large, larger FPS can be used to speed up the generation of the CNN features. The impact is that the resolution of the PVCD is reduced.

Then the Resnet29-RMAC CNN feature \cite{LAMV} is extracted for every frame. In one ablation study, we also extract the VGG16 max-pooling feature. The dimension of feature $d=512$.

We use the global frame feature database and the temporal network for copy segment localization. The parameters we use are ($top\_K\_all$, $top\_K\_one$, $top\_K\_video$, $sim_{th}$, $max\_step$, $max\_diff$), whose meanings are defined in Sections 3.1 and 3.2. For this study, we do not use the KNN search method. Instead, matrix multiplication and sorting are used to return the top-K matches.

In the training of the transformer encoder, we use a multi-head transformer with inner dimension 128. We use the Adam optimizer, learning rate 1E-6 and 30 epochs. The weight factor 0.1 and 1.1 are used on the ones and zeros elements in the similarity matrix.

In the training of the diagonalization CNN, we use the Adam optimizer, learning rate 1E-4, and 30 epochs. No weight factors are used on the ones and zeros elements in the similarity matrix.

For the speed test, we use an Intel i7-6850K CPU at 3.60GHz and a Ubuntu 16.04 LTS OS. We use an Nvidia GTX 1080TI GPU in the GPU mode of the KNN speed test.

\subsection{Dataset}

 VCDB \cite{VCDB1} is the only public PVCD dataset.  VCDB includes 528 videos in the core dataset and a very large background dataset. There are a total of 9236 annotated copy segments. In this section, we evaluate our algorithm on the VCDB dataset. We only use the videos in the core dataset and use all the annotated copy segments as the query. As in the standard VCDB protocol, the detected segment with any frame overlap with the ground truth is called a true positive. We only use segment F1 score in this study. One frame per second is extracted uniformly in all videos.

 However, for training the transformer, we cannot use the same VCDB dataset for both training and test since the whole VCDB dataset is used in the benchmark performance test \cite{VCDB1,VCDB2,LAMV}. Instead, we use a private PVCD dataset as a training dataset and use the whole VCDB dataset as a test dataset. This makes the benchmark on the VCDB dataset more reasonable and also demonstrates this algorithm's cross-domain performance. This dataset has more than 300 reference videos, most of which are longer than the VCDB reference videos. There are a few thousand query videos, each of which has a copy segment in multiple reference videos.

\subsection{Training Data Collection for Transformer}

In ViSiL \cite{VISIL}, \cite{VISIL-Ext}, and \cite{Transformer-VideoReloc}, triplet loss or contrastive learning is used, therefore triplet videos including anchor video, positive video, and negative video are prepared. In our work, there is no triplet video needed. We use the standard supervised learning using training samples (positive or negative) and their ground truth. So the collection of training data is easier.

The positive samples are prepared from the annotations of the training dataset. We borrow the idea of using hard negative samples in ViSiL \cite{VISIL}. We first do the PVCD on the features before the transformer is applied, collect all the false-positive results as the hard negative samples. If the numbers are not enough, we can collect more negative samples by randomly choosing video segments in the reference videos in the training dataset. This has been proved much better than only randomly choosing video segments in the reference videos.

Furthermore, since the number of negative samples is a lot more than the positive samples, we randomly pick negative samples from its pool in every epoch of the training process and make the numbers of positive and negative video segments equal in any epoch.

\subsection{KNN Search in Global Feature Database}

In the following subsection, we study the global feature database and the CNN feature without using the transformer or the similarity CNN.

After the candidate video is found, there are two methods to form the similarity matrix.  We can use use the original similarity matrix calculated from every pair-wise frame between the query segment and the candidate video, or we can use the reconstructed similarity matrix using the returned KNN results. Only the frames found by the KNN belonging to this candidate video have a non-zero value in the similarity matrix; others are all filled with zero. As a result, this similarity matrix is very sparse, and the temporal network localization runs a lot faster than the original similarity matrix.  If the performance is comparable, then the reconstructed similarity matrix is preferred for its fast speed.

Other parameters we try to optimize include $top\_K\_all$, and ($top\_K\_one$, $top\_K\_video$, $sim_{th}$, $max\_step$). To expedite the tuning, we only use the first 500 query segments.  Please note that, for a fair comparison, we must use the same set of query segments in this comparison. We should not compare these results with results for a different set of query segments.  The F1 score on the full query dataset may be lower than that on the first 500 segments. This is understandable because the distribution of the video similarity score of the first 500 segments is not stable yet.

We first use the reconstructed similarity matrix to optimize parameters, and the results are listed in the first panel in Table \ref{T0}.  From the results, we see that $top\_K\_video$=20 have better results than $top\_K\_video$=10, because in VCDB dataset many segments have copy in more than 10 reference videos. The similarity threshold $sim_{th}$ = 0.5,0.6 give good results. The best result of F1=0.9010 is achieved with ($top\_K\_one$, $top\_K\_video$, $sim_{th}$, $max\_step$)=(20.20.0.5,5).

Next, we check the performance of the original similarity matrix with a few good parameter configurations from the previous step, and the results are listed in the second panel in Table \ref{T0}. The best F1=0.8988 is achieved with $top\_K\_one$=20, $top\_K\_video$=20, $sim_{th}$=0.6, $max\_step$=5. From these results, we see that the reconstructed similarity matrix can achieve the same or even better performance than the original similarity matrix.

Thirdly, we test scanning the reference videos. In this mode, the original similarity matrix is used, and $top\_K\_video$ is not used. The results denoted by S are shown in the third panel in Table \ref{T0}. We notice the big improvement from $sim_{th}$=0.0 to $sim_{th}$=0.5, which demonstrates the effectiveness of the modified temporal network localization algorithm. Furthermore, the speed of $sim_{th}$=0.5 is a lot ($>50$ times) faster than $sim_{th}=0.0$ or using no threshold in the original temporal network.

\begin{table}
	\begin{center}
		\caption{Tuning parameters using the first 500 query segments. In the matrix column, R indicates reconstructed similarity matrix, O indicates original similarity matrix, and S indicates scanning videos. The parameters are  ($top\_K\_one$,$top\_K\_video$,$sim_{th}$,$max\_step$).}
		\label{T0}
		\begin{tabular}{ccc}
		    \hline
			Matrix & Parameters & F1\\
			\hline
			R & (5, 10, 0, 10) & 0.7917 \\
			R & (10, 10, 0, 10) & 0.8015 \\
			R & (10, 20, 0, 10) & 0.8584 \\
			R & (20, 20, 0, 10) & 0.8609 \\
			R & (20, 20, 0, 5) & 0.8609 \\
			R & (20, 20, 0.25, 10) & 0.8637 \\
			R & (20, 20, 0.5, 10) & 0.8980 \\
			R & (20, 20, 0.5, 5) &  \textbf{0.9010} \\
			R & (20, 20, 0.6, 10) & 0.8959 \\
			\hline
			O & (20, 20, 0, 10) & 0.8774 \\
			O & (20, 20, 0.25, 10) & 0.8654 \\
			O & (20, 20, 0.5, 5) & 0.9010 \\
			O & (20, 20, 0.6, 10) &  0.8959\\
			O & (20, 20, 0.6, 5) &  \textbf{0.8988}\\
			\hline
			S& (20, - , 0.5, 10) &  0.7352\\
			S& (20, - , 0.0, 10) &  0.6694\\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\subsection{Other CNN Features}

Other than the Resnet-29 RMAC feature from the LAMV \cite{LAMV}, we also use the VGG-16 pool5 feature map \cite{VGG16}. We do not directly use the FC output since it is widely believed that the Conv layer feature works better than the FC feature in the image retrieval task. We test the following pooling method: Max pooling, average pooling, Crop, SPoC, and find that the max-pooling gives the best F1 score. The results are summarized in Table \ref{T1}. For all these different pooling methods, please refer to \cite{Hu2020PyRetri} for details. The Resnet-29 RMAC F1 score outperforms the VGG-16 max-pooling by a big margin.

We analyze the performance improvements in the F1 score. Using the better feature is a big factor in it. Changing the feature from VGG-16 max pooling to the Resnet-29 RMAC, the F1 score improves from 0.7657 to 0.8613.  However, with the same Resnet-29 RMAC feature, in Table 2, the F1 score is only 0.7352 by scanning all reference videos with a $sim_{th}=0.5$. With $sim_{th}=0$, the F1 score is even worse 0.689. This shows the improvement due to the modified temporal network localization algorithm.  While with the global feature database, the F1 score improves to 0.8613. With all these factors added together, our PVCD framework achieves significantly improved performance.

\begin{table}
	\begin{center}
		\caption{Different CNN features and pooling methods.}
		\label{T1}
		\begin{tabular}{cc}
		    \hline
			CNN Feature & F1 Score\\
			\hline
			AlexNet FC, scanning videos \cite{VCDB2} & 0.6503\\
			\hline
			VGG-16 Max pooling & 0.7657 \\
            VGG-16 Avg pooling & 0.7656 \\
            VGG-16 SPoC & 0.7538 \\
            VGG-16 Crow & 0.7538 \\
            Resnet-29 RMAC & \textbf{0.8613} \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\subsection{CNN Features with Transformer}

In this section, we add the transformer encoder on the Resnet-29 RMAC CNN feature. We study to determine if this extra transformer encoder can improve the F1 score further over the previous two subsections. We use the previous F1 = 0.8613 as a baseline.

We test a 2-head and an 8-head transformer on the full VCDB query dataset. The only difference is that the original CNN feature in Figure \ref{fig1} is replaced with the transformer encoded feature. The F1 score results are listed in Table \ref{T2}, as well as the parameters used to achieve this F1 score. Please note that the speed does not include the extraction of frames from videos and the extraction of the original CNN features from frames, which are pre-processed and saved.

\begin{table}
	\begin{center}
		\caption{F1 scores and speed test results of transformer (TF) encoded feature on the VCDB query dataset. The parameters are $top\_K\_all$=200, ($top\_K\_one$, $top\_K\_video$, $sim_{th}$, $max\_step$, $max\_diff$). Transformers with 2-head (2H) and 8-head (8H) are tested. }
		\label{T2}
		\begin{tabular}{cccc}
		    \hline
			Method & Parameters & F1 score & speedup\\
			\hline
            Baseline & (20,20,0.5,5,5) & 0.8613 & 1\\
            2H TF & (3,20,0.55,2,2) & 0.8709 & 1.64 \\
            2H TF & (3,20,0.55,2,0) & 0.8714 & 1.68 \\
            8H TF & (3,20,0.55,2,2) & 0.8760 & 1.42 \\
            8H TF & (3,20,0.55,2,0) & \textbf{0.8764} & 1.59 \\
            \hline
            VGG16-MP & (20,20,0.5,5,5) & 0.7657 & - \\
            VGG16-MP,2H & (20,20,0.5,5,5) & 0.8077 & - \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

From the results, we notice that the transformer improves the F1 score from 0.8613 to 0.8760.  The improvement is not significant, perhaps because the previous baseline is already very good. In fact, the temporal network can model the temporal information in the copy localization algorithm indirectly, but it has to work on a pair of videos. The transformer, on the other hand, models the temporal information directly. The results show that the combination of the transformer and the temporal network localization algorithm brings some extra benefit, demonstrating that some temporal information is missed in the temporal network localization algorithm alone.

However, if we pay attention to the parameters, we notice that a lot smaller parameters ($top\_K\_all$, $top\_K\_one$, $max\_step$, $max\_diff$) are used, which can help increase the speed. Please note that $max\_diff$ = 0 means the network flow is only allowed along the diagonal in the similarity matrix. It indicates that a simple dynamic programming algorithm can be used for copy segment localization.

\textbf{Ablation study}: We use the 2-head and an 8-head transformer as an ablation study. The results show that a 2-head transformer is good enough. More heads get only a small performance gain while the model size gets larger quickly. We also test different $d_k$, inner dimensions and find that they have a very small impact on the performance.

We also test the transformer encoder on the VGG16 max-pooling (denoted VGG16-MP in Table \ref{T2}) feature, and the results are listed in Table \ref{T2}. The original F1 score is 0.7657. After the transformer encoder, the F1 score increases to 0.8077. This more significant improvement shows the effectiveness of transformer encoder on CNN feature not as good as the Resnet29-RMAC feature \cite{LAMV}.

\begin{figure*}[t]
    \centering
    \includegraphics[scale=0.4]{pvcd-example.png}
    \caption{Examples, the first row is a successful detected case, the second row is a failed case.}
    \label{fig4}
\end{figure*}

\textbf{Qualitative example}: Shown in Figure \ref{fig4} are examples of a successful and a failed case. In the successful case in the first row, the video on the left is significantly shaded and distorted. The target video on the right is flipped but is still partially detected. In the failed case in the second row, the target video has a relatively small actual content window. This failed case indicates that there is still work to do in PVCD.

\subsection{Comparison with the State of the Art}

We use our best F1 score on the full query dataset in Table \ref{T1} and Table \ref{T2} as our benchmark and compare it with the previous state-of-the-art results.  The comparison is listed in Table \ref{T3}. We notice that our CNN feature with transformer outperforms all existing methods. The F1 score improves from the previous 0.8025 to 0.8764. In addition, our algorithms improve the copy detection speed, as already shown in Table \ref{T2}, and more details of the speed test of KNN with GPU will be given in the next subsection.

\begin{table}
	\begin{center}
		\caption{Comparison with existing methods on VCDB dataset.}
		\label{T3}
		\begin{tabular}{cc}
		    \hline
			Method &  F1 Score\\
			\hline
			CNN \cite{VCDB2} & 0.6503 \\
			LAMV \cite{LAMV} & 0.6740 \\
			Ref \cite{MYLiu} & 0.6440 \\
			Ref \cite{MMTA} & 0.8025  \\
			Ours (No-transformer) & 0.8613  \\
			Ours (Transformer) & \textbf{0.8764}  \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\subsection{KNN Speed Test}

In this section, we test the query speed of different search methods.  The time reported here is only the query search time, not including the network flow localization time.

The baseline search method we test is a brute force search and sorting in matrix multiplication. Then we test the FAISS index types: Flat-CPU, Flat-GPU, IVF-Flat-GPU \cite{FAISS}. Since the total number of frames in the VCDB dataset is less than 100,000, we do not test the scalar or product quantization. The average search time per frame, the speedup over the matrix multiplication are summarized in Table \ref{T4}.  We see that the speedup of the Flat-CPU is about eight times, and that of Flat-GPU and IVF256-Flat-GPU 15 ad 11 times, respectively.

\begin{table}
	\begin{center}
		\caption{Average search time per frame}
		\label{T4}
		\begin{tabular}{ccc}
		    \hline
			Method &  time(ms) & Speed up \\
			\hline
			Matrix multiplication & 20.0 & 1  \\
			Flat-CPU & 2.55 & 7.86 \\
			Flat-GPU & 1.48 & 14.79 \\
			IVF256-Flat-GPU & 1.82 & 10.81 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\section{Conclusion}

In this work, we study a global feature database and a fast KNN search method for PVCD to find candidate videos that likely have copy segments in them. Using Resnet29-RMAC \cite{LAMV} feature, we improve the previous F1 score around 0.80 to 0.8613.

Then we study an algorithm to improve the similarity of video pairs in PVCD. The transformer encoder enhances the representation of the frame CNN feature. The transformer improves the F1 score on the full VCDB query dataset further to 0.8764.

Lastly, we study a modified temporal network to localize copy segments in the candidate videos. The modifications corporate with the improved representation or video similarity to improve both the F1 score and the speed.

Overall, we improve the F1 score of PVCD on the VCDB dataset from the previous state-of-the-art 0.8025 to the new 0.8764.

\bibliographystyle{ieee_fullname}