\documentclass[onecolumn,pra,superscriptaddress,nofootinbib,longbibliography,preprint]{revtex4-2}
\usepackage[letterpaper, left=1in, right=1in, top=1.2in, bottom=1.2in]{geometry}
\usepackage{titlesec}

\usepackage{cmap} % Load before fontenc 
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{bm}
\usepackage{xcolor}
\definecolor{blueviolet}{rgb}{0.2, 0.2, 0.6}
\definecolor{webgreen}{rgb}{0,.5,0}
\definecolor{webbrown}{rgb}{.6,0,0}
\usepackage[pdftex,
	bookmarks=false,
	colorlinks=true, %allcolors=blueviolet,
	urlcolor=webbrown, 
	linkcolor=blueviolet, 
	citecolor=webgreen,
	pdfstartpage=1,
	pdfstartview={FitH},  % FitBH
	bookmarksopen=false
	]{hyperref}
\usepackage{tikz}
\usepackage{natbib}
\usepackage{mathtools}
\usepackage{graphicx}
\usepackage{csquotes}
\usepackage{braket}
\usepackage{dsfont}
\usepackage{listings}
\usepackage{chemformula}
\usepackage{hyperref}
% \usepackage{url}
\allowdisplaybreaks


\usepackage{multirow}
\usepackage{amssymb}

\usepackage{amsthm}

\usepackage{color}
\usepackage{nicefrac}


\DeclareFixedFont{\ttb}{T1}{txtt}{bx}{n}{9} % for bold
\DeclareFixedFont{\ttm}{T1}{txtt}{m}{n}{9}  % for normal

\newcommand{\addsection}[3]{\addtocontents{toc}{\protect\contentsline{section}{\protect\numberline{#1}#2}{#3}}}
\newcommand{\addsubsection}[3]{\addtocontents{toc}{\protect\contentsline{subsection}{\protect\numberline{#1}#2}{#3}}}

% Custom colors
\usepackage{color}
\definecolor{deepblue}{rgb}{0,0,0.5}
\definecolor{deepred}{rgb}{0.6,0,0}
\definecolor{deepgreen}{rgb}{0,0.5,0}

% Python style for highlighting
\newcommand\pythonstyle{\lstset{
language=Python,
basicstyle=\ttm,
morekeywords={self},              % Add keywords here
keywordstyle=\ttb\color{deepblue},
emph={MyClass,__init__},          % Custom highlighting
emphstyle=\ttb\color{deepred},    % Custom highlighting style
stringstyle=\color{deepgreen},
frame=tb,                         % Any extra options here
showstringspaces=false
}}

% Python environment
\lstnewenvironment{python}[1][]
{
\pythonstyle
\lstset{#1}
}
{}

% Python for external files
\newcommand\pythonexternal[2][]{{
\pythonstyle
\lstinputlisting[#1]{#2}}}

% Python for inline
\newcommand\pythoninline[1]{{\pythonstyle\lstinline!#1!}}

%\usepackage{float}

\usepackage{xcolor}
\newcommand{\revise}[1]{\textcolor{black}{#1}}
\newcommand{\modify}[1]{\textcolor{red}{#1}}


\newcommand{\jp}[1]{{\color{blue}\textbf{JP:}~#1}}
\newcommand{\jm}[1]{{\color{blue}\textbf{JM:}~#1}}
\newcommand{\rb}[1]{{\color{red}\textbf{RB:}~#1}}
\definecolor{orange}{RGB}{255,127,0}
\newcommand{\zj}[1]{{\color{orange}\textbf{ZJ:}~#1}}
\newcommand{\change}[1]{{\color{blue}#1}}

\newcommand{\robert}[1]{\textcolor{orange}{[Robert: #1]}}
\newcommand{\masoud}[1]{{\color{blue}\textbf{MM:}~#1}}
%\newcommand{\masoud}[1]{\ignorespaces}
%\newcommand{\masoud}[1]{}

\def\avg#1{\mathinner{\langle{#1}\rangle}}
\def\bra#1{\ensuremath{\mathinner{\langle{#1}|}}}
\def\ket#1{\ensuremath{\mathinner{|{#1}\rangle}}}
\def\dbra#1{\ensuremath{\mathinner{\langle\langle{#1}|}}}
\def\dket#1{\ensuremath{\mathinner{|{#1}\rangle\rangle}}}
\newcommand{\ketbra}[2]{\lvert #1 \rangle \! \langle #2 \rvert}
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
\newcommand{\dbraket}[2]{\langle \langle #1|#2\rangle \rangle}
\newcommand{\expval}[1]{\langle #1\rangle}
\newcommand{\tr}{\text{Tr}}
\newcommand{\ns}{{N_s}}


\newcommand{\cA}{{\mathcal{A}}}
\newcommand{\cI}{{\mathcal{I}}}
\newcommand{\cL}{{\mathcal{L}}}
\newcommand{\cB}{{\mathcal{B}}}
\newcommand{\cT}{{\mathcal{T}}}
\newcommand{\cC}{{\mathcal{C}}}
\newcommand{\cP}{{\mathcal{P}}}
\newcommand{\cW}{{\mathcal{W}}}
\newcommand{\cX}{{\mathcal{X}}}
\newcommand{\cY}{{\mathcal{Y}}}
\newcommand{\cZ}{{\mathcal{Z}}}
\newcommand{\cD}{{\mathcal{D}}}
\newcommand{\cE}{{{\mathcal{E}}}}
\newcommand{\cF}{{\mathcal{F}}}
\newcommand{\cM}{{\mathcal{M}}}
\newcommand{\cN}{{{\mathcal{N}}}}
\newcommand{\CN}{{\mathcal{N}}}
\newcommand{\cS}{{\mathcal{S}}}
\newcommand{\cG}{{{\mathcal{G}}}}
\newcommand{\cU}{{{\mathcal{U}}}}
\newcommand{\cQ}{{\mathcal{Q}}}

\newcommand{\indicator}{\mathds{1}}
\newcommand{\euler}{\mathrm{e}}
\newcommand{\rmi}{\mathrm{i}}

\newcommand{\wt}{\widetilde}

\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\sign}{sign}
\DeclareMathOperator{\SWAP}{SWAP}

\newtheorem{example}{Example}
\newtheorem{proposition}{Proposition}
\DeclareMathOperator{\Tr}{tr}
\DeclareMathOperator*{\E}{{\mathbb{E}}}
\DeclareMathOperator*{\Var}{\mathrm{Var}}
\DeclareMathOperator{\median}{\mathrm{median}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{definition}{Definition}
\newtheorem{lemma}{Lemma}
\newtheorem{fact}{Fact}
\newtheorem{task}{Task}
\newtheorem{remark}{Remark}

\newcommand{\robnote}[1]{\textcolor{blue}{[Robert: #1]}}

\newcommand\blfootnote[1]{%
  \begingroup
  \renewcommand\thefootnote{}\footnote{#1}%
  \addtocounter{footnote}{-1}%
  \endgroup
}

\newcommand{\Id}{I}
\newcommand{\vct}[1]{\mathbf{#1}}

\usepackage[noend]{algpseudocode}
\usepackage{algorithm,algorithmicx}

\algrenewcommand\alglinenumber[1]{\sf\scriptsize\color{blue}{#1}}
\algrenewcommand\algorithmicrequire{\textbf{Input:}}
\algrenewcommand\algorithmicensure{\textbf{Output:}}


\begin{document}

(Chinese abstract here.)

\newpage


% \title{Hardware-Efficient and Accurate Unsupervised Quantum Machine Learning on Electronic Structure Problems}

% \author{Chia-Tung Chu}
%     \email{b07202059@ntu.edu.tw}
% 	\affiliation{Department of Physics and Center for Theoretical Physics, National Taiwan University, Taipei 10617, Taiwan}
% \author{Shih-Kai Chou}
% 	\affiliation{Department of Physics and Center for Theoretical Physics, National Taiwan University, Taipei 10617, Taiwan}
% \author{Hsi-Sheng Goan}
%     \email{goan@phys.ntu.edu.tw}
% 	\affiliation{Department of Physics and Center for Theoretical Physics, National Taiwan University, Taipei 10617, Taiwan}
%     \affiliation{Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan}
%     \affiliation{Physics Division, National Center for Theoretical Sciences, Taipei, 10617, Taiwan}
% \date{\today}
% \date{September 23, 2022}

\begin{abstract}
    As the rapid evolution of quantum hardware industry, various algorithms in quantum computational chemistry are put forward for rigorous implementation. 
    However, due to lack of quantum error correction, a promising hybrid quantum-classical algorithm to find the lowest eigenvalue of a given Hamiltonian for near-term quantum hardware is the variational quantum eigensolver (VQE). 
    To get the best performance out of limited computational resources, using efficient basis sets to represent molecular Hamiltonian in electronic structure problems has become a crucial step. 
    Some particular VQE calculations using Daubechies wavelet basis sets and efficient encoding schemes implemented via hardware-efficient parameterized quantum circuits (PQCs) on real quantum machines or via chemical-inspired PQCs on simulators to date are only accomplishing chemical accuracy on relatively small and simple structure molecules.

    In this paper we present a hybrid quantum-classical \revise{machine learning} computational framework that can be trained unsupervisedly to reliably predict molecular vibrational frequencies.
    A novel \revise{scheme} for encoding the interatomic distances \revise{of the molecules} into \revise{the} quantum circuits \revise{of the quantum-classical framework} is also proposed and can allow for a fewer-shot prediction of the molecular ground state \revise{potential energy curve}.
    We demonstrate with $\ch{H2}$, $\ch{LiH}$ and $\ch{PN}$ molecules that our framework exhibits high energy accuracy \revise{with the calculated results of vibrational frequencies in excellent agreement with experimental data} while remains short in quantum circuit depth as the number of qubits scales up.
    Our approach provides a powerful \revise{hybrid quantum-classical computational} architecture for the calculation of molecular \revise{properties} on near-term quantum \revise{computer} processors.
    The preliminary results have been accepted as poster by the premier conference in Quantum Information Science, 26th Quantum Information Processing.
\end{abstract}


\maketitle

\clearpage

\tableofcontents

\section{Introduction}

Simulating many-body physics is one of the promising applications on near-term quantum computers, among which is solving the electronic structure problems efficiently, as compared to classical computers.
\revise{For example, the computational complexity of exact quantum
chemical computations scales exponentially with system
size on classical computers, while it is in a polynomial
manner for variational hybrid quantum-classical algorithms \cite{georgescu2014quantum, mcclean2016theory, mcardle2020quantum, abrams1999quantum}.}
Most recent studies on solving electronic problem with quantum computers focus primarily on approximating the ground state energy curve of given molecules.
\revise{Here we present a hybrid quantum-classical machine learning computational framework that can be trained unsupervisedly to reliably predict molecular vibrational frequencies to compare with experimental data.}

%For instance, quantum algorithms providing exponential speedup for finding eigenvalues as well as the theory of variational hybrid quantum-classical algorithms \cite{georgescu2014quantum, mcclean2016theory, mcardle2020quantum, abrams1999quantum} have nourished the era of quantum computational chemistry.~\cite{armaos2019computational, mcardle2020quantum}
%We are, however, elevating the emphasis of this work to the calculation of molecular vibrational frequencies, which is a more relevant physical quantity regarding real-world purpose.
%Specifically speaking, it allows us to gain better understanding of infrared absorption as well as the mechanisms and kinetics of chemical reactions. 

%In general, quantum computers are naturally adept at calculating molecular vibrational frequencies since it can store and manipulate many-body wave function with an exponential speedup over its classical counterparts.

\revise{In the noisy intermediate-scale quantum (NISQ) era~\cite{preskill2018quantum}, }
the available quantum devices to date are relatively small and prone to errors.
Without the presence of fault-tolerant quantum computing, we need algorithms that require less resource, while not sacrificing the accuracy when solving electronic structure problems.
Common hardware-efficient algorithms generally embed variational quantum circuits (VQC) as part of their workflows, and often feature a two-local~\cite{brandao2016local} design \revise{in the VQC ansatz} in order to approximate the eigenvalue of a given qubit Hamiltonian.
The overall circuit with purely unitary operations, under which quantum states evolve, are necessarily linear.
We are, however, engaging in a quantum machine learning (QML) task, which not only ground state energies of the given Hamiltonian are obtained, but also the energies of more unseen Hamiltonians, needed to be predicted.
Therefore, non-linearity must be introduced into our quantum circuit to increase its expressibility.


Various ways exist for realizing non-linear operations in quantum circuit. For example, in the quantum convolutional neural network (QCNN)~\cite{cong2019quantum}, they utilize partial measurements to inject non-linearity into the circuit.
Circuit of this kind, extending toward the last measurement, has a much longer depth compared to that of ours with global measurement, \revise{and} hence in our work the circuit is more likely to run on a NISQ device.
Besides utilizing the global measurement as means for achieving nonlinear operations as stated in~\cite{holmes2022connecting}, several other techniques have been introduced in our work to increase the trainability.
%As stated in~\cite{holmes2022connecting}, 
\revise{For example, we propose a new encoding scheme and post an additional constraint on our circuit parameters during the encoding stage as a kind of a} regularizer, which eventually leads to a better convergence, and will be addressed later in the paper.


Generally speaking, our workflow can unsupervisedly train a hybrid \revise{variational} quantum circuit to calculate the ground state energy curve of a molecule, under a given set of bond lengths.
After training, the circuit with fixed parameters can predict the ground state \revise{energies} at unseen bond lengths, within a promising accuracy.
Combined with our regularizer during the encoding stage, accurate molecular vibrational frequency can be obtained, while only few quantum resource are used.
We will demonstrate our frequency calculation subroutine with $\ch{H2}$, $\ch{LiH}$ and $\ch{PN}$ molecules and discuss their performance.


The outline of the paper is as follows. 
In Sec.~\ref{sec:formulation}, we formally define the electronic structure problem and the cost function associated with the optimization process.
In Sec.~\ref{sec:workflows}, we give an illustration of our general framework regarding an advanced problem-solving procedure.
In Sec.~\ref{sec:encoding}, we introduce the notion of adaptive encoding, and in Sec.~\ref{sec:hamiltonian} we mention the better quality of the Hamiltonian used in our work as compared to the Gaussian minimal basis sets.
In Sec.~\ref{sec:restore}, we show the capability of our hybrid network to approximate the molecular dissociation curve within chemical accuracy, with relatively few quantum resource required.
In Sec.~\ref{sec:freq_predict} and \ref{sec:compare_encoding}, our main results show great \revise{agreement} with experimental data, and different encoding schemes are compared to demonstrate the advantage of our proposed method.
We conclude and aim on several future directions regarding our work in Sec.~\ref{sec:conclusion}.



\section{General framework}

\subsection{Problem formulation}
\label{sec:formulation}

The framework we proposed in this paper mainly resolves a predefined electronic structure problem.
The problem is formulated as follows: given a Hamiltonian of a many-body quantum system, the goal is to approximate its associated minimum eigenvalue, which is essentially the ground state energy of the system.
Whenever there is a \revise{change in} configuration of the quantum system, a new Hamiltonian must be obtained, hence a new ground state energy. 
In our problem, only the distance between the constituent atoms is varied. 
Therefore the calculation of ground state energies at all possible bond lengths is equivalent to recovering the entire ground state \revise{potential} energy curve, or the dissociation curve of that specific molecule. 
Once we have calculated adequate number of points near the global minimum, we can obtain the molecular vibrational frequency by taking the second order derivative of the ground state energy curve with respect to bond lengths.


In previous works~\cite{mcardle2019digital}, one would normally obtain the ground state energy of the concerned molecule at several bond lengths individually, and then derive the molecular vibrational frequency from five to seven consecutive energy points near \revise{the} global minimum. 
Considering only the hardware-efficient methods other than the Unitary Coupled Cluster (UCC) method~\cite{lee2018generalized}, normal workflow requires the invoking of new parametrized quantum circuit (PQC) for each individual Hamiltonian, 
that is, the determination of a unitary $U(\pmb\theta)$ which represents the overall quantum operation that our circuit has done on the input state, with $\pmb\theta$ being a vector of parameters that completely specifies the action of every parametrized gate in our PQC.


A common way of constructing $U(\pmb\theta)$ is to utilize the set of \revise{parameterized}single-qubit rotational Pauli gates $\{R_x(\theta), R_y(\theta), R_z(\theta)\}$, combined with some multi-qubit \revise{entangling} gates (such as \texttt{CNOT}) as the building blocks for our circuit.
The circuit can then be trained on its parameters according to some optimization subroutines.
Notice that the trained vector of rotation angles $\pmb\theta = \{\theta_j\}$ has to be discarded every time when one is calculating the ground state energy $C_{\rho, H}(\pmb\theta)$ of a new Hamiltonian $H$,
\begin{equation}\label{eq:hamiltonian}
    C_{\rho, H}(\pmb\theta) = \Tr[HU(\pmb\theta)\rho U(\pmb\theta)^\dagger],
\end{equation}
where $\rho = (\ket0\bra0)^{\otimes n}$ is an $n$-qubit input state.
In order to obtain enough points for vibrational frequency calculation, the minimization of cost function $C_{\rho, H}(\pmb\theta)$ must be repeated several times, and the same number of parameter vectors $\pmb\theta$ all serve no purpose after a series of point-wise energy calculations.



\subsection{Workflows}
\label{sec:workflows}


In our work, we decide to take the physical value of the trained parameter vector $\pmb\theta$ into account. 
We realize it by considering the loss function comprised of multiple Hamiltonian at once,
\begin{equation}\label{eq:sumover_h}
    C_{\rho, \{H_i\}}(\pmb\theta) = \sum_i\Tr[H_i\tilde U(\pmb\theta)\rho \tilde U(\pmb\theta)^\dagger],
\end{equation}
which is evaluated on a single quantum circuit. 
Notice that we have also introduced non-linear operations to the circuit for more expressibility, and hence the operator $\tilde U(\pmb\theta)$ need not be unitary. 
In addition, we anticipate our hybrid quantum-classical circuit capable of predicting ground state energies of unseen Hamiltonian $H(l)$ at different bond lengths after training, therefore the information of bond length $l$ must somehow be encoded into the circuit. 
We denote the bond length encoding scheme by
\begin{equation}\label{eq:encode}
    \rho = \rho(l) = \tilde G(l)(\ket0\bra0)^{\otimes n}\tilde G(l)^\dagger,
\end{equation}
where the encoding circuit $\tilde G(l)$, which depends on specific bond lengths, again need not be unitary. 
Now we are ready to train our hybrid network based on the above discussed new cost function formalism.



Fig.~\ref{fig:fig1} shows the unrolled general framework of our computational model. 
The model is a multi-layer structure, each layer is essentially a parametrized quantum circuit, followed by a global Pauli-$Z$ measurement at the end. 
The only information transmitting between any two adjacent layers are the set of Pauli-$Z$ expectation values $\{\expval{\sigma_z}_p\}, \forall 1\leq p\leq n$, for a $n$-qubit quantum circuit. 
In Fig.~\ref{fig:fig1}(a) we illustrate a $K$-layer structure, where all parameters $\{\theta\}$ in each layer can be updated independently and simultaneously, according to the multi-dementional gradient $\frac{\partial L}{\partial\pmb\theta}$ of the loss function
\begin{equation}\label{eq:loss}
    L_{\{l_i\}, \{H_i\}}(\pmb\theta) = \sum_i\Tr[H_i\tilde U(\pmb\theta) \tilde G(l_i) \rho \tilde G(l_i)^\dagger \tilde U(\pmb\theta)^\dagger],
\end{equation}
which is calculated by the classical optimizer. 
During each iteration, a new set of parameters $\{\theta\}$ are obtained from some pre-defined optimization algorithm, and are updated to the quantum circuits for the evaluation of loss function during next iteration.


\begin{figure}[!h]
    \centering
    \includegraphics[width=.9\textwidth]{Fig1.pdf}
        \caption{
            Schematic illustration of the hybrid quantum circuit framework with $n$ qubits, solid lines indicate quantum wires, dotted lines indicate classical interactions with the quantum circuit. 
            \textbf{(a)} The high-level general framework with loss-evaluation and parameter-updating routine. Note that the three stacks of $\expval{H}$ calculation share the same set of parameters. 
            \textbf{(b)} Each single layer is comprised of an encoding part (quantum circuit), a variational quantum circuit, and a Pauli-$Z$ measurement for each qubit. The measurement outcome is then passed to the next layer. 
            \textbf{(c)} The general circuit for encoding the input scalars from previous layer (bond lengths or $\expval{\sigma_z}$. The red arrows represent functions that maps from input to rotation angles. 
            \textbf{(d)} The variational quantum circuit with hardware-efficient ansatz.  Linear entanglement pattern is adopted, and its repetition can be adjusted according to the desired expressibility.
        \label{fig:fig1}
        }
\end{figure}


Note that we have encoded the information carried by bond lengths into the circuit at the beginning of the framework, where technical details will be discussed in Sec.~\ref{sec:encoding}. 
At the end of the $K^\text{th}$ layer, a series of Pauli string representing the qubit Hamiltonian is attached to the circuit such that we can measure the expectation value of $\{H_i\}$, thus able to calculate the loss function defined in Eq.~\ref{eq:loss}. 
Notice that all trainable parameters are embedded in VQC part, serving as the angle for $R_y$ gates. The hardware-efficient ansatz illustrated in Fig.~\ref{fig:fig1}(d) can be fine-tuned so that the repetition number $m$ of the 2-local circuit gives our network sufficient expressibility, while maintaining its trainability.


\subsection{Encoding}
\label{sec:encoding}

There are various ways for encoding classical bond lengths to a quantum state. 
In our work, the encoding part is universal in all layers, that is, it handles not only molecular bond length, but also the vector of Pauli-$Z$ measurement outcomes from the previous layer. 
Therefore our encoding scheme is equivalent to the task of finding a suitable set of functions $\{f_i\}$, such that the encoded $n$-qubit state
\begin{align}\label{eq:psi_encode}
    \ket\psi_{\text{encode}} &= \Pi_{i=0}^{n-1}\Big(\otimes R_y(\theta_i) H \Big) \ket{\psi_0} \nonumber \\
    &= \Pi_{i=0}^{n-1}\Big(\otimes R_y\big(f_i(\vec a)\big) H \Big) \ket0^{\otimes n},
\end{align}
has peak performance in vibrational frequency prediction, given a vector of scalar $\vec a$. 
The vector can be comprised of either specific bond length $\vec a = (l,\cdots,l)$, or measurement outcomes from the previous layer $\vec a = (\expval{\sigma_z}_0, \cdots ,\expval{\sigma_z}_{n-1})$.


Consider the encoding formalism described in Eq.~\ref{eq:encode}, one would naturally come up with a function that represents the direct encoding scheme~\cite{xia2020hybrid}
\begin{equation}\label{eq:direct}
    \theta_i = f_i(\vec a) = a_i,
\end{equation}
which is relatively simple, yet effective. 
Nevertheless, as we observe the continuity of the exact molecular ground energy curve, a concern has risen. 
Generally we need to ensure the smoothness of the ground energy curve, in order to have a reliable frequency calculation. 
This implies that the energy difference between adjacent bond lengths must be suppressed under a threshold, and this tiny energy gap can be accounted to the degree of overlapping of $\ket\psi_{\text{encode}}$ from different bond length, $l_i$ and $l_{i+1}$. 
Hence our goal is to find a new set of function that takes 
\begin{equation}\label{eq:new_function}
    \arg\max_{\{f_i\}}\sum_j \big|\braket{\psi_j | \psi_{j+1}}_\text{encode}\big|^2,
\end{equation}
where $j$ sums over all dinstinct bond lengths. Hence, we propose the \emph{adaptive encoding} scheme.

In the adaptive encoding, one makes use of the bond length resolution $\Delta l$ before training. 
Suppose we have picked five bond lengths for frequency prediction, $\vec l = \{0.70, 0.72, 0.74, 0.76, 0.78\}$ Å. 
The resolution in this case would be $\Delta l = |l_0-l_1| = 0.02$ Å, and the set of encoding functions are
\begin{equation}\label{eq:adaptive}
    \theta_i = f_i(\vec a) = a_i + (i-\frac{n-1}{2})(\frac{\Delta l}{2})
\end{equation}
for a $n$-qubit circuit. 
For instance, the parameter vector $\vec\theta$ which encode $\vec a = (0.74,\cdots,0.74)$ in a direct scheme will remain $(0.74,0.74,0.74,0.74)$ for a 4-qubit circuit, 
but in the adaptive encoding scheme, the two bond lengths $l_2 = 0.74$ and $l_3=0.76$ will be encoded as $\vec\theta(l_2) = (0.71, 0.73, 0.75, 0.77)$ and $\vec\theta(l_3) = (0.73, 0.75, 0.77, 0.79)$, respectively. 
We will see in Sec.~\ref{sec:compare_encoding} that the latter encoding scheme has significant improvement on the stability of vibrational frequency prediction. 



\subsection{Hamiltonian Preparation and Methods}
\label{sec:hamiltonian}

In this work we adopted the Daubechies wavelet basis function \cite{hong2022accurate} for the construction of molecular Hamiltonian, 
due to its capability of yielding accurate results with excellent computational efficiency, as compared with contracted Gaussian basis sets on large molecules.~\cite{mohr2015accurate, ratcliff2020flexibilities}
We first calculate the matrix elements, which are the orbital integrals in second quantized Hamiltonians, by using the BigDFT~\cite{genovese2008daubechies, mohr2014daubechies, genovese2011daubechies} package on a classical computer. 
Next we transform them into qubit Hamiltonians, hence compatible with quantum processors, during which we also adopt the two-qubit reduction~\cite{steudtner2018fermion}.
Therefore, a $n$ qubit Hamiltonian is transformed into a $n-2$ qubit Hamiltonian after encoding.
Overall, we obtain three sets of 2, 4, and 6 qubit Hamiltonians for $\ch{H2}$, $\ch{LiH}$, and $\ch{PN}$ molecules, respectively.

The optimization is performed by the Broyden–Fletcher–Goldfarb–Shanno algorithm~\cite{broyden1970convergence} with maximum 500 iterations and gradient norm tolerance to stop as $10^{-5}$. 
%In the simulation, the expectation of the operator is simulated by the \texttt{qiskit.opflow} module and the Hamiltonian can be treated as a single operator, the circuit architecture is essentially an extension of Qiskit's 2-local ansatz.
We run our hybrid networks on IBM's 27-qubit quantum machine, \texttt{ibmq-montreal}, for all 3 types of molecules ($\ch{H2}$, $\ch{LiH}$, and $\ch{PN}$).


We calculate the vibrational frequency based on the harmonic oscillator model of diatomic molecule, which requires the coefficient from the second-order fitting of the energy curve. 
To obtain a reliable coefficient, we need at least 5 consecutive points of energy on the energy curve. 
We also initialize all the variational parameters with a Gaussian distribution of mean 0 and standard deviation $\pi/2$ before training. 
All other molecular configuration datas are from the \href{https://cccbdb.nist.gov}{Computational Chemistry Comparison and Benchmark DataBase}.

\section{Results}
\label{sec:results}

\subsection{Restoring ground energy curve}
\label{sec:restore}

\begin{figure}[!h]
    \centering
    \includegraphics[width=0.8\textwidth]{Fig2.pdf}
        \caption{
            The capability of our proposed framework to recover the entire ground state energy curve, showcased for $\ch{H2}$ (2-qubit Hamiltonian) and $\ch{LiH}$ (4-qubit Hamiltonian). 
            \textbf{(a)} We used $m=1$ and $K=2$ architecture for the $\ch{H2}$ molecule, and 
            \textbf{(b)} We used the same $m=1$ and $K=2$ architecture for the $\ch{LiH}$ molecule. 
            The trained energies, as well as the predicted ones, fitted well with the exact diagonalization of Hamiltonian. 
            All predicted energies (testing points) converge within the chemical accuracy (indicated as shaded regions in the two bottom plots).
        \label{fig:fig2}
        }
\end{figure}

\begin{table}[!h]
    \centering
    \caption{Framework configurations for the \revise{potential} energy curve calculation of $\ch{H2}$ and $\ch{LiH}$. 
    For the VQC part of the circuit we adopted linear entanglement throughout the entire experiment.
    We trained both molecules with 10 energy points.
    ($K$ denotes the number of layers used in our framework, and $m$ denotes the repetition of ansatz in the VQC part of each layer, as illustrated in Fig.~\ref{fig:fig1}.)}
    \begin{ruledtabular}
    \begin{tabular}{cccc}
        Molecule & Hamiltonian & Range (Å) & $K, m$ \\
        \hline
        $\ch{H2}$ & 2 qubit & 0.4 $\to$ 2.48 & 2, 1 \\
        $\ch{LiH}$ & 4 qubit & 1.4 $\to$ 1.8 & 2, 1 \\
    \end{tabular}
    \end{ruledtabular}
    \label{tab:tab1}
\end{table}

Fig.~\ref{fig:fig2} shows the capability of our circuit to recover the entire ground state energy curve with only few training points. 
Table~\ref{tab:tab1} displays the hyperparameters we used for the calculation. 
In both cases ($\ch{H2}$ and $\ch{LiH}$), we only utilize 2 layers and one repetition of hardware-efficient ansatz in our circuit setting. 
We can see from the plots that all of the energy points, whether corresponds to the bond lengths for training or testing, converge within the chemical accuracy ($0.0016$ Ha). 
As for the purpose of predicting vibrational frequency, we require even fewer training points, since the relevant prediction range shrinks toward the proximity of global minimum. 
That is, a qualified calculation of frequency need not the information of the entire curve.


\subsection{Predicting Frequency}
\label{sec:freq_predict}

\begin{figure*}[!h]
    \centering
    \includegraphics[width=1.0\textwidth]{Fig3.pdf}
        \caption{
            Ground state energy prediction near global minimum for frequency calculations. 
            We train our hybrid network on only two distinct bond lengths for each molecule (indicated as hexagonal markers), and let it predict the energy of neighbor bond lengths (testing points, indicated as star markers) after training. 
            We fit a set of consecutive energy points to second-order polynomial functions (shown as wide solid lines), then extract the coefficient and derive the harmonic vibrational frequency for each molecule. 
            The results display great match with frequencies that are derived from the exact diagonalization of Hamiltonians.
            \textbf{(a)} For $\ch{H2}$, we computed energy for 5 points, after fitting we obtain $f_{\text{H}_2} = 4460.78 \pm 2.16 \text{cm}^{-1}$.
            \textbf{(b)} For $\ch{LiH}$, we computed energy for 7 points, after fitting we obtain $f_\text{LiH} = 1424.37 \pm 4.09 \text{cm}^{-1}$.
            \textbf{(c)} For $\ch{PN}$, we computed energy for 5 points, after fitting we obtain $f_\text{PN} = 1368.55 \pm 9.70 \text{cm}^{-1}$.
            We conducted the experiment 100 times for $\ch{H2}$ and $\ch{LiH}$ molecules to mitigate the error from randomness (hence 100 slightly distinct frequencies generated for both molecules), 
            and 50 times for $\ch{PN}$ molecule, due to the growing queueing time as the qubit number increases. 
            Next, we calculate the average vibrational frequency of each molecule, and serve them as our final predictions.
            Overall, the above statistics show good align with our anticipated network performance.
            (Note that we have rescaled the $y$-axis of $\ch{H2}$ and $\ch{PN}$ for better visualization, where the scaling factor for recovering the true energy ($Y$) from the printed energy ($y$) is anchored at the top-left corner of each plot.
            The relations are $Y = 10^{-5}\cdot y - 1.133 (\text{Hartree})$ for $\ch{H2}$, and $Y = y - 16.11 (\text{Hartree})$ for $\ch{PN}$.)
        \label{fig:fig3}
        }
\end{figure*}


\begin{table*}[!hbt]
    \centering
    \caption{Framework configurations as well as results for the vibrational frequency calculation of $\ch{H2}$ , $\ch{LiH}$ and $\ch{PN}$. 
    We adopted linear entanglement \revise{for the VQC ansatz throughout the entire simulation process to reduce the circuit depth.}
    $f_\text{predict}$, $f_\text{diag}$, and $f_\text{exp}$ represent the vibrational \revise{frequencies from our network prediction, the exact diagonalization of Hamiltonian, and the experimental data}, respectively, and $\text{error}_\text{diag}$ and $\text{error}_\text{exp}$ denote the \revise{deviations} of our predicted frequency from \revise{the values of the exact diagonalization and experimental data,} respectively.
    }
    \begin{ruledtabular}
    \begin{tabular}{cccc}
        Molecule & $K, m$ & $l_\text{training}$ (Å) & $l_\text{testing}$ (Å) \\
        \hline
        $\ch{H2}$, 2 qubit & 2, 1 & 0.73, 0.74 & 0.725, 0.735, 0.745 \\
        $\ch{LiH}$, 4 qubit & 2, 1 & 1.51, 1.53 & 1.49, 1.50, 1.52, 1.54, 1.55 \\
        $\ch{PN}$, 6 qubit & 3, 2 & 1.47, 1.49 & 1.46, 1.48, 1.50 \\
    \end{tabular} 
    \begin{tabular}{cccccc}
        Molecule & $f_\text{predict}\,(\text{cm}^{-1})$ & $f_\text{diag}\,(\text{cm}^{-1})$ & $f_\text{exp}\,(\text{cm}^{-1})$ & $\text{error}_\text{diag}$ & $\text{error}_\text{exp}$ \\
        \hline
        $\ch{H2}$, 2 qubit & 4460.78 $\pm$ 2.16 & 4462.40 & 4401.21 & -0.03\% & 1.35\% \\
        $\ch{LiH}$, 4 qubit & 1424.37 $\pm$ 4.09 & 1426.23 & 1405.50 & -0.13\% & 1.34\% \\
        $\ch{PN}$, 6 qubit & 1368.55 $\pm$ 9.70 & 1365.63 & 1336.95 & 0.21\% & 2.36\% \\
    \end{tabular}
    \end{ruledtabular}
    \label{tab:tab2}
\end{table*}

We present our results on vibrational frequency prediction in Fig.~\ref{fig:fig3}.
The framework settings we used during the work are listed in Table~\ref{tab:tab2}.


For $\ch{H2}$ and $\ch{LiH}$ molecule, we repeated the predicting procedure for 100 times on quantum machine. 
For $\ch{PN}$ we repeated for 50 times due to the long queueing time. 
Table~\ref{tab:tab2} shows that the results are consistent with the exact values obtained from Hamiltonian diagonalization. 
We train our hybrid network on only 2 bond lengths for all three molecules (indicated as blue hexagonal markers in Fig.~\ref{fig:fig3}), and left the other for prediction.
Overall, we generated 5, 7, and 5 energy points for $\ch{H2}$, $\ch{LiH}$ and $\ch{PN}$ respectively, which are essential for second-order function fitting, and hence the calculation of frequency.
We can see from Table~\ref{tab:tab2} that the predicted frequencies are very accurate.
Note that even as we increased the ansatz repetition and the layer number for 6-qubit $\ch{PN}$ to increase the network's expressibility,
the overall circuit depth still remains shallow.


\subsection{Comparison of Encoding}
\label{sec:compare_encoding}

\begin{figure}[!h]
    \centering
    \includegraphics[width=.9\textwidth]{Fig4.pdf}
        \caption{
            \textbf{(a)} Scatter plot for the vibrational frequency calculation of LiH molecule on a 4-qubit hybrid network. 
            We conduct the experiment for 100 times on both direct encoding and adaptive encoding scheme, the advantage of the latter over the former has been statistically visualized on the plot. 
            In terms of frequency error rate, adaptive encoding has a value of $-0.014\pm0.018$, which is significantly more stable and more accurate than that of the direct encoding ($-0.034\pm0.169$).
            \textbf{(b)} The loss function trajectory over iterations for both encoding schemes. 
            We can see that the two methods take essentially the same order of magnitude of steps to converge, making the adoption of adaptive encoding from direct encoding with nearly no cost.
        \label{fig:fig4}
        }
\end{figure}

We are curious about the degree of performance improvement that the adaptive encoding scheme brings about eventually. 
In Fig.~\ref{fig:fig4} we show the difference between our proposed encoding method and the direct encoding scheme. 
From the scatter plot we notice that the adaptive encoding scheme guarantees a more stable prediction of vibrational frequency than the original direct encoding scheme.
Such stableness is visualized from a much higher concentration of blue star markers near the zero-error line, as compared to the sparse distribution of orange dots in Fig.~\ref{fig:fig4} (a).
Whereas the iteration counts remain same for both schemes, making the adoption of adaptive encoding from direct encoding with nearly no cost, as illustrated in Fig.~\ref{fig:fig4} (b).


\section{Conclusion and Outlook}
\label{sec:conclusion}

In this work, we proposed a hybrid framework that shows great potential in frequency prediction for various molecules. 
We have shown that our hybrid quantum-classical neural network, combined with adaptive encoding scheme can be trained to obtain the electronic energies at certain bond lengths and then generate the whole potential energy curve. 
The results of $\ch{H2}$, $\ch{LiH}$ and $\ch{PN}$ are very accurate and the corresponding errors are over an order of magnitude smaller compared to previous work~\cite{xia2020hybrid}.
The circuit's width and depth are also shortened due to the advantage of Daubechies wavelet basis set over the Gaussian based basis set, and the adoption of global measurement over partial measurement as non-linear operations.
The fact that our network succeeds at predicting ground state energy on previously unseen bond lengths may implies that it has been 
trained to memorize the relation between individual bond length and their corresponding ground state wave function (or at least the state with energy approximating the true ground state energy).


Furthermore, we show the importance of problem-specific encoding technique in comparison with direct encoding scheme. 
Often we are not really concerned about the absolute value of the ground state energy,
\begin{figure}[h]
    \centering
    \includegraphics[width=.8\textwidth]{Fig5.pdf}
        \caption{
            The adaptive encoding method proposed in our work serves as one of the means to establish correlation between energies of neighboring bond lengths.
            As one can see, high correlation makes energies (blue marks) distribute less randomly along the fitted curve (blue curve).
            Thus, the figure shows the advantage of adaptive encoding over direct encoding, under the same environmental settings ($\ch{LiH}$, 4-qubit Hamiltonian).
        \label{fig:fig5}
        }
\end{figure}
but are rather more interested in the energy difference, when the interatomic distance is varied, 
for example, the molecular vibrational frequency near equilibrium bond length. 
Through the design of new encoding methods, we aim to maintain the correlation between energies under various bond lengths, which is illustrated in Fig.~\ref{fig:fig5}. 
The idea is similar to the regularization technique in the context of classical machine learning.
This approach may also be generalized to accomplish the encoding of two-dimensional molecular configuration data in future work. 
On the other hand, the intermediate nonlinear measurements can reduce the circuit depth , which are more suitable for NISQ devices, and the layer-structure proposed in this work may serve as a fundamental building block for hybrid VQCs in the near-term future. 


During our work, Barren plateau has emerged as the qubit number scales up. 
From Ref.~\cite{holmes2022connecting} we know that \emph{strongly correlating parameters} and/or initializing close to the solution (and then restricting the ansatz to explore the region close to the initialization) are the most effective approaches to avoid exponentially vanishing cost gradients. 
Thus, in the future work, we can post a stronger constraint on the circuit parameters to maintain their correlation, to avoid low trainability. 
We will also keep on to investigate the performance of our proposed framework on quantum hardware with the presence of noise in the future.


\vspace{-1em}
\subsection*{Acknowledgments:}
\vspace{-1em}

The author thanks Prof. Hsi-Sheng Goan for invaluable guidance and support throughout this project.
The author also thanks Jyh-Pin Chou, Yuan-Chung Cheng for great suggestions and discussions.
The author also thanks IBM Quantum Hub at NTU for providing computational resources and accesses for conducting the experiments on real quantum device.
Chia-Tung Chu is supported by the National Science and Technology Council under Grants No. MOST 111-2813-C-002-009-M.


\bibliography{references}

\clearpage

\section{Appendix}

\subsection{Hardware specifications}

All quantum circuits in our work are run on the IBM's 27-qubit quantum computer, \texttt{ibmq-montreal}.
The calibration datas of relevant qubits are listed in Table~\ref{tab:tab3}, while Fig.~\ref{fig:fig6} shows the hardware's qubit map view.
The experiments are chosen to conduct on qubits \#\{1,2,3,4,5,8\}, so that the actual CNOT gates number can be minimized hence mitigate the errors.

\begin{table}[!h]
    \centering
    \caption{
        Calibration data for 27-qubit quantum machine \texttt{ibmq-montreal} with regard to all qubits used.
        For each qubit, six types of relevant properties are listed, including relaxation time (T$_1$), dephasing time (T$_2$), qubit operation frequency, anharmonicity, Pauli-X error and CNOT error.
        We run our quantum circuits on qubits \#\{1,4\} for $\ch{H2}$, qubits \#\{1,2,3,4\} for $\ch{LiH}$, and qubits \#\{1,2,3,4,5,8\} for $\ch{PN}$, respectively. 
        All listed values have been averaged within the time span of our experiments, from January 14th, 2023 to February 26th, 2023. 
        The notation for CNOT error $A_B$ represents the error of CNOT gate connecting qubits \#A and \#B.
        }
    \begin{ruledtabular}
    \begin{tabular}{ccccccc}
        Qubit & T$_1$       & T$_2$       & Frequency & Anharmonicity  & Pauli-X error & CNOT error \\
        \#    & ($\mu$s) & ($\mu$s) & (GHz)     & (GHz)          & $(10^{-4})$   & $(10^{-3})$ \\
        \hline
        1 & 161.8354 & 30.7795 & 4.8327 & -0.3247 & 2.1961 & $1_4$:9.336354071; $1_2$:8.006437467 \\
        2 & 119.6453 & 140.7872 & 4.9821 & -0.3386 & 2.3030 & $2_3$:7.8294181379; $2_1$:8.006437467 \\
        3 & 64.8695 & 41.4315 & 5.1040 & -0.3352 & 3.8469 & $3_5$:11.082183326; $3_2$:7.8294181379 \\
        4 & 146.0053 & 131.9639 & 5.0024 & -0.3387 & 2.0058 & $4_1$:9.336354071 \\
        5 & 124.2324 & 105.1834 & 5.0339 & -0.3372 & 2.6063 & $5_3$:11.082183326; $5_8$:6.1232952963 \\
        8 & 92.3221 & 99.0577 & 4.9072 & -0.3243 & 2.4756 & $8_5$:6.1232952963
    \end{tabular} 
    \end{ruledtabular}
    \label{tab:tab3}
\end{table}

\begin{figure}[!h]
    \centering
    \includegraphics[width=.8\textwidth]{Fig6.pdf}
        \caption{
            Map view of the 27-qubit \texttt{ibmq-montreal} quantum machine, 
            we exploited six connected qubits (\#\{1,2,3,4,5,8\}) during our experiments.
        \label{fig:fig6}
        }
\end{figure}


\subsection{Error mitigation}

Quantum computers are subject to various types of errors that can occur due to environmental factors such as thermal noise, imperfect control pulses, and decoherence. 
Therefore, quantum error mitigation is a crucial part in our work, as it is essential for achieving reliable circuit operations, 
since untreated errors can significantly affect the accuracy of quantum computations, leading to incorrect results.
To address these challenges, we exploited various quantum error mitigation techniques in \texttt{Qiskit Runtime} primitive (provided by IBM) to minimize the impact of errors on our proposed framework. 
In the following subsections, we briefly elaborate all three techniques that are used in our experiments.

\subsubsection{Twirled Readout Error eXtinction (TREX)}

The Twirled Readout Error eXtinction method reduces measurement error by diagonalizing the noise channel associated with measurement by randomly flipping qubits via $X$ gates immediately prior to measurement, 
and flipping the corresponding measured bit if an $X$ gate was applied. 
A rescaling term from the diagonal noise channel is learned by benchmarking random circuits initialized in the $\ket{0}$ state. 
Consequently, the bias from expectation values which is resulted from readout noise can be removed. 


\subsubsection{Zero Noise Extrapolation (ZNE)}

The Zero Noise Extrapolation method computes an expectation value of the observable for different noise factors (amplification stage), 
and then uses the measured expectation values to infer the ideal expectation value at the zero-noise limit (extrapolation stage). 
This procedure is visualized in Fig.~\ref{fig:fig_resiliance}.
This approach tends to reduce errors in expectation values, but is not guaranteed to produce an unbiased result.
The overhead of this method scales with the number of noise factors (which is 3 times overhead in our experiment's setting).

\begin{figure}[h]
    \centering
    \includegraphics[width=.55\textwidth]{resiliance-2.pdf}
        \caption{
            Illustration of the Zero Noise Extrapolation (ZNE) method.
        \label{fig:fig_resiliance}
        }
\end{figure}


\subsubsection{Probabilistic Error Cancellation (PEC)}

The Probabilistic Error Cancellation method mitigates error by learning and inverting a sparse noise model that is able to capture correlated noise. 
PEC returns an unbiased estimate of an expectation value so long as learned noise model faithfully represents the actual noise model at the time of mitigation. 
In practice, the experimental procedure for learning the noise model has ambiguities due to certain error terms that cannot be independently distinguished. 
These are resolved by a symmetry assumption, which depending on the true underlying noise may lead a biased estimate of the mitigated expectation values due to using an imperfect noise model.
In our experiments' setting, this method specifically addresses noise in self-inverse two-qubit gates, so it first stratifies each input circuit into an alternating sequence of simultaneous 1-qubit gates followed by a layer of simultaneous 2-qubit gates. 
Then it learns the noise model associated with each unique ywo-qubit gate layer, as illustrated in Fig.~\ref{fig:fig_resiliance}.

\begin{figure}[h]
    \centering
    \includegraphics[width=.75\textwidth]{stratified.png}
        \caption{
            An example of a stratified circuit to achieve Probabilistic Error Cancellation, where the layers of two-qubit gates are labeled layer $1$ through $n$. 
            Note that each $\mathcal{U}_i$ is composed of two-qubit gates on the native connectivity graph of the quantum processor. 
            The open boxes represent arbitrary single-qubit gates.
        \label{fig:fig_stratified}
        }
\end{figure}

The overhead of this method scales with the number of noise factors. 
In our experiments the expectation value is sampled at three noise factors, leading to a roughly 3-times overhead.

Since PEC uses a quasi-probability method to mimic the effect of inverting the learned noise,
it requires sampling from a randomized circuit family associated with the user’s original circuit. 
Therefore, applying PEC will increase the variability of the returned expectation value estimates unless the number of samples per circuit is also increased for both input and characterization circuits. 
The amount of samples required to counter this variability scales exponentially with the noise strength of the mitigated circuit, thus it is a relatively expensive error-mitigation method.



\subsection{Unrolled circuit diagram}


To give a concrete example of the circuit being run on the hardware, in Fig.~\ref{fig:fig6}
we explicitly print out the quantum part of our hybrid network for $\ch{PN}$ molecule with the same parameters ($K=3$ and $m=2$) listed in Table~\ref{tab:tab2}. 

\begin{figure}[!h]
    \centering
    \includegraphics[width=1.0\textwidth]{unroll_1.pdf}
        \caption{
            The quantum part (unrolled) of our hybrid network for $\ch{PN}$ molecule, which is essentially a 6-qubit circuit
            with framework parameters $K=3$ and $m=2$ (referring to Fig.~\ref{fig:fig1}).
        \label{fig:fig7}
        }
\end{figure}

\clearpage

\section{Publication List}

\noindent[1] \underline{Chia-Tung Chu}, Shi-Kai Chou, and Hsi-Sheng Goan,
"Hardware-Efficient and Accurate Unsupervised Quantum Machine Learning on Electronic Structure Problems", 
in \textit{proceedings of 26th Quantum Information Processing}, 
Ghent, Belgium, February 2023. 
\href{https://indico.cern.ch/event/1175020/page/28188-accepted-posters}{[link]}

\end{document}  


