\section{Introduction}

Automatic differentiation is a technique for constructing a set of
instructions for computing the derivatives of a function defined by
a computer program. It is implemented in software packages
such as Autograd \cite{autograd}, PyTorch \cite{pytorch}, TensorFlow \cite{tensorflow} and Jax \cite{jax2018github}, and is now a crucial component in
machine learning libraries.
It is natural that this functionality
in classical computation is being replicated in quantum computing
frameworks. An increasing number of libraries, e.g., Qiskit \cite{Qiskit}, TensorFlow Quantum \cite{tfq}, and PennyLane \cite{bergholm2020pennylane}, provide automatic computation of gradients and higher derivatives
mainly for optimization and machine learning applications. Quantum differentiable
programming has also been implemented specifically for quantum chemistry
in PennyLane \cite{arrazola2021differentiable}.

Automatic differentiation has a broader range of uses than the applications
mentioned above. In this work we consider the use of automatic differentiation
to solve a problem in condensed matter physics, namely, the calculation
of fidelity susceptibility \cite{PhysRevE.76.022101,PhysRevB.76.104420}.
This is a quantity
related to derivatives of a wave function overlap and has been used in the
study of quantum phase transitions in a variety of condensed matter models.
Our focus here is not on the physics of this problem but primarily to provide
an end-to-end demonstration of quantum automatic differentiation to solve it.
The PennyLane library \cite{bergholm2020pennylane} is used and the
cross-platform capability of this framework is illustrated by running
PennyLane code on Qiskit simulators. Finally, this study exposes some of the
issues that have to be dealt with in the NISQ era of quantum computation.
The computation had to be designed to reduce sensitivity to noise as much as
possible since gate error mitigation with automatic differentiation is still
a significant challenge. We implement such mitigation and perform an analysis
of the complications that arise when doing so in a practical setting. To that end this work presents, to the best of our knowledge, the first study of how applying error mitigation affects the computation of energy derivatives with respect to Hamiltonian parameters, which we find to be a very sensitive process even for small systems and a limited amount of mitigation.

In Sect. \ref{subsec:quantum_gradients} we briefly review how quantum gradients can be computed automatically using parameter-shift rules. In Sect. \ref{sec:fs} the concepts of fidelity \cite{PhysRevE.74.031123} and fidelity susceptibility \cite{PhysRevE.76.022101,PhysRevB.76.104420}
as used in condensed matter physics are introduced. They are a measure
of how a wave function changes when Hamiltonian parameters change.
Therefore, they can be useful in the study of quantum systems which
exist in different phases depending on interaction parameters. The
calculation of fidelity susceptibility is done in a variational framework
and its derivation in terms of derivatives with respect to variational
parameters is outlined. In order to complete the calculation, one
needs to know how variational parameters change with respect changes in
Hamiltonian parameters. For this, derivatives of energy expectation values are
required. This is also discussed in Sect. \ref{sec:tfim}.

The model used in this study is the one-dimensional transverse field
Ising model. It is the simplest model which displays a quantum phase
transition and is discussed, for example, in \cite{Suzuki2013}. The
model is exactly solvable by transforming it to a model of free fermions
\cite{PFEUTY197079}
but that would defeat our purpose, so we work with the spin degrees of
freedom directly. In simulations including a hardware noise model
only small systems (4 or 6 spins) are feasible. The techniques
used to simplify the calculation and reduce noise sensitivity are
discussed in Sect.~\ref{sec:tfim} and results of simulations are
presented in Sect.~\ref{sec:simres}. A demonstration of automatically-differentiable error
mitigation is presented in Sect.~\ref{sec:error_mitigation}. Alternate algorithms for calculating wave function overlaps are discussed in Appendix A and supplementary error mitigation results are presented in Appendix B.

\section{Background}

\subsection{Quantum gradients and automatic differentiation}
\label{subsec:quantum_gradients}

Let $U(\theta)$ be a parametrized (variational) quantum circuit with $\theta = \{\theta_i\}$ a set of real-valued parameters. Consider a quantum computation that begins with a set of qubits in the all-zeros state $\ket{0\cdots0}$, executes $U(\theta)$, and measures the expectation value of a Hamiltonian $H$. The analytical expression of the expectation value is a function of $\theta$,
\begin{equation}
E(\theta) = \langle H \rangle = \langle 0 \cdots 0 | U^\dagger (\theta)  H U(\theta) |0\cdots 0 \rangle.
\end{equation}

\noindent This function can be differentiated with respect to the variational parameters. A common means of achieving this is using \emph{parameter-shift rules}.

Suppose we would like to compute the gradient of $E(\theta)$ with respect to a particular $\theta_i$. If $\theta_i$ appears in the circuit as the parameter of a single-qubit rotation, it is by now well-known that the gradient can be obtained by running the circuit at two different values of the parameter and combining the results according to the shift rule \cite{mitarai2018quantum, li2017hybrid, wierichs2022general, schuld2019evaluating, mari2021estimating}:
\begin{equation}
\frac{\partial E(\theta)}{\partial \theta_i} = \frac{1}{2}\left[ E(\ldots, \theta_i + \pi / 2, \ldots) - E(\ldots, \theta_i - \pi/2, \ldots) \right],\label{eq:param_shift}
\end{equation}
\noindent where it is implicit that all parameters other than $\theta_i$ are held constant during this process.

Such rules are advantageous as they allow for computation of gradients in both software and on hardware without modifying the actual structure of the circuit, as only the input parameters change. For more sophisticated unitary operations, generalized parameter-shift rules can be used \cite{wierichs2022general}, which may require evaluating the circuit at four, or even more shifted values, depending on the operation.

Parameter-shift rules provide a straightforward recipe for computing quantum gradients; they can be further differentiated in order to compute similar rules for higher-order derivatives, such as Hessians \cite{mari2021estimating}. The availability of these rules thus allows gradient computation to be incorporated into automatic differentiation frameworks, wherein parameters can be modified through classical processing (e.g., multiplication by a constant or by another parameter), and the chain rule can be applied. Results of fully-differentiable quantum computations can also be incorporated into more extensive hybrid classical-quantum computations.

\subsection{Fidelity susceptibility}\label{sec:fs}

Phase transitions are a ubiquitous feature of many-body systems. The
long-standing paradigm of describing phase transitions in terms of a
local order parameter associated with symmetry breaking \cite{itzykson_drouffe_1989}
is not adequate for describing many to the systems studied in contemporary
condensed matter physics \cite{sachdev_2011,Bernevig+2013}. Many alternative
ways of describing phase transitions have been proposed including
fidelity \cite{PhysRevE.74.031123} and fidelity susceptibility \cite{PhysRevE.76.022101,PhysRevB.76.104420}.

Consider a Hamiltonian with a single variable parameter
\begin{equation}
H(r)=H_{0}+rH_{1}. \label{eq:basic_ham}
\end{equation}
 When $r$ is varied the system may be in different phases. The fidelity,
as defined by Chen \emph{et al.} \cite{PhysRevA.77.032111} is
\begin{equation}
F(r,\delta)=\left|\left\langle \psi(r)|\psi(r+\delta)\right\rangle \right|,
\end{equation}
 where $\psi(r)$ and $\psi(r+\delta)$ are eigenstates (only ground
states are considered here) at two different values of the Hamiltonian
parameter. The susceptibility is the second derivative of this quantity
w.r.t. $\delta$ evaluated at $\delta$ = 0. By expanding $\psi(r+\delta)$
in eigenstates of $H(r)$ one can derive a kind of spectral representation
of the fidelity susceptibility which involves the sum over a complete
set of states. For the ground state, the expression for the fidelity
susceptibility is \cite{PhysRevA.77.032111}
\begin{equation}
\mathcal{S}(r)=\sum_{n\neq0}\frac{\left|\left\langle \psi_{n}(r)\left|H_{1}\right|\psi_{0}(r)\right\rangle \right|^{2}}{\left|E_{0}(r)-E_{n}(r)\right|^{2}}.
\end{equation}
This, in principle, requires knowledge of the complete spectrum. However,
the use of automatic differentiation allows for the calculation of
the fidelity susceptibility using only the ground state wave function.

We suppose that the ground state wave function has been obtained by
a variational calculation. Let $\theta = \{\theta_{i}\}$ denote the variational
parameters which at the variational minimum are implicit functions
of $r.$ The quantity to be computed is
\begin{equation}\label{eq:d2over}
\mathcal{S}(r)=\left|\partial_{r}^{2}\left\langle \psi|\psi(\theta(r),r)\right\rangle \right|
\end{equation}
where it is understood that differentiation is applied only to the
right side of the overlap. By the chain rule
\begin{equation}
\partial_{r}\left\langle \psi|\psi(\theta(r),r)\right\rangle =\sum_{i}\frac{\partial\left\langle \psi|\psi(\theta(r),r)\right\rangle }{\partial\theta_{i}}\frac{\partial\theta_{i}}{\partial r},
\end{equation}
 and
\begin{equation}\label{eq:d2chain}
\partial_{r}^{2}\left\langle \psi|\psi(\theta(r),r)\right\rangle =\sum_{i,j}\frac{\partial^{2}\left\langle \psi|\psi(\theta(r),r)\right\rangle }{\partial\theta_{j}\partial\theta_{i}}\frac{\partial\theta_{i}}{\partial r}\frac{\partial\theta_{j}}{\partial r}+\sum_{i}\frac{\partial\left\langle \psi|\psi(\theta(r),r)\right\rangle }{\partial\theta_{i}}\frac{\partial^{2}\theta_{i}}{\partial r^{2}}.
\end{equation}
With a parametrization of the variational wave function and a circuit
to compute the overlap, the derivatives of
the overlap w.r.t. the variational parameters can be computed in a straightforward manner through parameter-shift rules. Obtaining derivatives
of the variational parameters w.r.t. the Hamiltonian parameter $r$
is more involved.

For the transverse field Ising model (considered in this work) or other Hamiltonians for which
the expansion of the ground state wave function involves coefficients
which are real, the gradient of the overlap w.r.t. variational parameters
will vanish so the second term in (\ref{eq:d2chain}) will not contribute.
If the expansion coefficients are complex this may not be the case.

The Hadamard test is a standard method to compute the overlap of two wave
functions; automatic differentiation of the Hadamard test circuit to
get the Hessian needed in the first term in (\ref{eq:d2chain}) is a way to
get to
the fidelity susceptibility. Alternatively, consider $\left|\left\langle \psi|\psi(\theta(r),r)\right\rangle \right|^{2}$
and, if the gradient $\frac{\partial\left\langle \psi|\psi(\theta(r),r)\right\rangle }{\partial\theta_{i}}$
vanishes, then
\begin{equation}\label{eq:d2ovsqr}
\frac{\partial^{2}\left|\left\langle \psi|\psi(\theta(r),r)\right\rangle \right|^{2}}{\partial\theta_{j}\partial\theta_{i}}=2\frac{\partial^{2}\left\langle \psi|\psi(\theta(r),r)\right\rangle }{\partial\theta_{j}\partial\theta_{i}}
\end{equation}
 is another means by which to compute overlap Hessian. The choice of algorithm
for computing the right-hand side of (\ref{eq:d2ovsqr}) is discussed
in Sect. \ref{sec:simres}  and in Appendix A.

In order to complete the calculation of the fidelity susceptibility
the quantities $\frac{\partial\theta_{i}}{\partial r}$ are needed.
To get these, derivatives of the Hamiltonian expectation value are
used. See the article by Pulay \cite{pulay1987} for a systematic discussion.
The energy is
\begin{equation}
E=\left\langle \psi(\theta(r),r)\left|H\right|\psi(\theta(r),r)\right\rangle
\end{equation}
 and the condition for the variational minimum is
\begin{equation}\label{eq:var}
\frac{\partial E}{\partial\theta_{i}}=\frac{\partial\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}}=0
\end{equation}
 for all $i$. The first energy derivative is
\begin{eqnarray}
\frac{\partial E}{\partial r} & = & \left\langle \psi\left|\frac{\partial H}{\partial r}\right|\psi\right\rangle +\sum_{i}\frac{\partial\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}}\frac{\partial\theta_{i}}{\partial r}\\
 & = & \left\langle \psi\left|H_{1}\right|\psi\right\rangle
\end{eqnarray}
at the variational minimum. The second energy derivative is
\begin{equation}\label{eq:d2E}
\frac{\partial^{2}E}{\partial r^{2}}=\sum_{i}\frac{\partial\left\langle \psi\left|H_{1}\right|\psi\right\rangle }{\partial\theta_{i}}\frac{\partial\theta_{i}}{\partial r},
\end{equation}
for a Hamiltonian which is linear in the variable parameter.

To get the derivatives of the variational parameters w.r.t. $r$, take
the derivative of the variational condition (\ref{eq:var}). The first
order response is
\begin{equation}
\frac{\partial}{\partial r}\left(\frac{\partial\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}}\right) = \frac{\partial\left\langle \psi\left|H_{1}\right|\psi\right\rangle }{\partial\theta_{i}}+\sum_{j}\frac{\partial^{2}\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}\partial\theta_{j}}\frac{\partial\theta_{j}}{\partial r} = 0 \label{eq:first_order_response}.
\end{equation}
 This gives a set of equations that can be solved to get the $\frac{\partial\theta_i}{\partial r}.$
This is sufficient for the calculation of the first and second energy
derivatives as discussed in \cite{arrazola2021differentiable, mitarai2020theory} and for
the fidelity susceptibility calculated here. If second derivatives of the
variational parameters are required they can be obtained by taking another
derivative of the variational condition \cite{pulay1987}
\begin{eqnarray}
\frac{\partial^{2}}{\partial r^{2}}\left(\frac{\partial\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}}\right) & = & \sum_{j}\frac{\partial^{2}\left\langle \psi\left|H_{1}\right|\psi\right\rangle }{\partial\theta_{i}\partial\theta_{j}}\frac{\partial\theta_{j}}{\partial r}+\sum_{j,k}\frac{\partial^{3}\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}\partial\theta_{j}\partial\theta_{k}}\frac{\partial\theta_{j}}{\partial r}\frac{\partial\vartheta_{k}}{\partial r}\nonumber \\
 &  & +\sum_{j}\frac{\partial^{2}\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}\partial\theta_{j}}\frac{\partial^{2}\theta_{j}}{\partial r^{2}},\nonumber \\
 & = & 0.
\end{eqnarray}

\noindent In summary, to compute the fidelity susceptibility (second derivative of the
overlap) and the second derivative of energy, one needs to compute:
\begin{enumerate}
 \item The Hessian of the overlap w.r.t. the variational parameters
 \item The Hessian of the expectation value of the full Hamiltonian w.r.t. the
   variational parameters
 \item The gradient of the expectation value of the $r$-dependent part of the
   Hamiltonian ($H_1$) w.r.t. the variational parameters
 \end{enumerate}

\noindent In this work we will demonstrate how all of these gradients can be computed using quantum automatic differentiation software.

\section{Fidelity susceptibility and the Transverse Field Ising Model}\label{sec:tfim}

This section outlines the analytical description of the quantum computation that will be performed. For the one-dimensional transverse field Ising model the Hamiltonian
\begin{equation}
H(r)=-\sum_{i=0}^{L-1}\left(\sigma_{i}^{x}\sigma_{i+1}^{x}+r\sigma_{i}^{z}\right)
\end{equation}
will be used with periodic boundary conditions, so there is a translation
invariance. Furthermore, the Hamiltonian is invariant under a transformation
by $\sigma^{z}$ at every site, that is,
\begin{equation}\label{eq:HP}
\left[H(r),\mathcal{P}\right]=0,
\end{equation}
 where
\begin{equation}
\mathcal{P}=\bigotimes_{i} \sigma_{i}^{z}.
\end{equation}

Consider a one-dimensional spin chain with an even number $L$ of sites.
The total number of possible spin configurations is $2^{L}.$ However,
due to (\ref{eq:HP}), only half of the possible spin configurations, i.e.,
$2^{(L-1)}$ contribute to the eigenstate of $H.$ For the ground state,
the contributing spin configurations are those where both the number
of up and down spins is even. Due to translation invariance many spin
configurations will have the same coefficient. These features are
used to simplify the quantum computation.

As an illustration of how the computation can be simplified consider
the case of 4 spins. The number of possible spin configurations is
16 but only 8 appear in the ground state wave function. These spin
configurations can be grouped into 4 composite basis states:
\begin{eqnarray*}
 & \left|0000\right\rangle, \\
 & \frac{1}{2}\left(\left|0011\right\rangle +\left|1001\right\rangle +\left|1100\right\rangle +\left|0110\right\rangle \right),\\
 & \frac{1}{\sqrt{2}}\left(\left|0101\right\rangle +\left|1010\right\rangle \right),\\
 & \left|1111\right\rangle .
\end{eqnarray*}
 Calculating the matrix elements of the full Hamiltonian using the
above states gives a reduced $(4\times4)$ Hamiltonian matrix
\begin{equation}
\left(\begin{array}{cccc}
-4r & 0 & -2 & 0\\
0 & 0 & -2\sqrt{2} & 0\\
-2 & -2\sqrt{2} & 0 & -2\\
0 & 0 & -2 & 4r
\end{array}\right).
\end{equation}
 It can be verified by exact diagonalization that the ground state energy
and fidelity susceptibility using the reduced Hamiltonian are identical
to the results obtained with the full Hamiltonian.

\begin{figure}[tb]
\centering
\includegraphics[scale=0.35,trim={6cm 4cm 6cm 4cm},clip=true,angle=-90]{./figs/fourspin.pdf}
\caption{Circuit for the 4-spin variational wave function encoded in two qubits. }
\label{fourspin}
\end{figure}

The quantum computation using the reduced Hamiltonian can be done
using two qubits with a variational wave function depending on three
parameters. The circuit for the unitary operator which instantiates the 4-spin
variational wave function is shown in Fig. \ref{fourspin}. In terms of Pauli operators
the Hamiltonian (19) acting on two qubits is
\[
-X(0)-X(1)-X(0)Z(1)+Z(0)X(1)-\sqrt{2}(X(0)X(1)+Y(0)Y(1))-2r(Z(0)+Z(1)).
\]

The above procedure can be extended to systems with more spins. For
6 spins, the 64 possible spin configurations can be reduced to a
composite basis of eight states. The quantum computation can be done
using three qubits. The circuit for the 6-spin variational wave function
depending on seven parameters is shown in Fig. \ref{sixspin}. This circuit has been optimized for a quantum architecture with only nearest-neighbour qubit connectivity.

\begin{figure}[tb]
\centering
\includegraphics[scale=0.50,trim={6cm 1cm 6cm 1cm},clip=true,angle= -90]{./figs/sixspin.pdf}
\caption{Optimized circuit for the 6-spin variational wave function encoded in three qubits. }
\label{sixspin}
\end{figure}

\section{Simulation results}\label{sec:simres}

Computations were carried out for 4- and 6-spin systems using quantum
simulators. Code was written using the PennyLane differentiable programming
framework. For production work, the PennyLane \texttt{default.qubit} and Qiskit
Aer simulators were used (through the PennyLane-Qiskit plugin)\footnote{To check cross-platform capability, code was a tested on Cirq and
Rigetti Forest simulators also.}. All code and data files used for this work is available on GitHub \cite{ourgithub}.

The first step in the computation was the determination of optimal
variational wave function parameters for a number of different values
of the Hamiltonian parameter $r$ spanning a range from 0.5 to 1.4 (this region was chosen as it is centered around the point of a phase transition at $r=1$).
This was done using a standard variational quantum eigensolver (VQE).
VQE was performed using fully analytical simulation. All variational parameters were
initialized to 0. A gradient descent optimizer
with step size 0.1 was used, with the maximum number of iterations set to 1000, however optimization would terminate early if the energy computed at a
particular iteration was within $10^{-8}$ of the true value obtained by exact diagonalization.

Both the fidelity susceptibility, Eq. (\ref{eq:d2over}), and
the second derivative of the energy (\ref{eq:d2E}) are of interest
as indicators of a phase transition
(see Ref.\cite{PhysRevA.77.032111}). As per Eqs. (\ref{eq:first_order_response}) and (\ref{eq:d2ovsqr}), to determine these quantities the
gradient
\begin{equation}\label{eq:grad}
\frac{\partial\left\langle \psi\left|H_{1}\right|\psi\right\rangle }{\partial\theta_{i}}
\end{equation}
 and the Hessians
\begin{equation}\label{eq:hess}
\frac{\partial^{2}\left\langle \psi\left|H\right|\psi\right\rangle }{\partial\theta_{i}\partial\theta_{j}},\enskip \thinspace\frac{\partial^{2}\left|\left\langle \psi|\psi(\theta(r),r)\right\rangle \right|^{2}}{\partial\theta_{j}\partial\theta_{i}}
\end{equation}
 are required. PennyLane computes these automatically by applying parameter-shift rules to circuits
for the Hamiltonian expectation value and squared wave function overlap. The circuit
employed for the square of the wave function overlap is shown in generic
form in Fig. \ref{uudag}. Starting in a state of all 0's, the probability
of measuring all 0's gives the squared overlap. If one could not make
use of Eq. (\ref{eq:d2ovsqr}) and the overlap, instead of its square,
was really required, a Hadamard test could be used. This is discussed in
Appendix A.

\begin{figure}[tb]
\centering
\includegraphics[scale=0.45,trim={4cm 10cm 1cm 10cm},clip=true,angle=0]{./figs/uudagger.pdf}
\caption{Circuit for calculating the squared overlap of the wave functions instantiated by the
unitary operators $U_i$ and $U_f$. }
\label{uudag}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/exactFS.pdf}
\caption{Fidelity susceptibility per site of the transverse field Ising model as a function
of the Hamiltonian parameter $r$ computed using a simulator of an ideal quantum
computer. The circles (squares) are for a lattice of four (six) sites. The solid
lines are exact results.}
\label{exactFS}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/exactSED.pdf}
\caption{Second derivative of the energy per site of the transverse field Ising model as a function
of the Hamiltonian parameter $r$ computed using a simulator of an ideal quantum
computer. The circles (squares) are for a lattice of four (six) sites. The solid
lines are exact results.}
\label{exactSED}
\end{figure}

To verify that automatic differentiation was working properly the fidelity
susceptibility (\ref{eq:d2over}) and the second derivative of the energy
(\ref{eq:d2E}) were computed at each value of $r$ using the PennyLane \texttt{default.qubit} device which
simulates an ideal (i.e., noise-free) quantum computer. The results, averaged
over twenty trials using 8192 shots for each measurement per trial, are shown
in Figs. \ref{exactFS} and \ref{exactSED} for spin chains of length 4 and 6.
The solid lines are the results obtained by exact solution of the model.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/L4FS.pdf}
\caption{Fidelity susceptibility per site of the transverse field Ising model with four
sites as a function of the Hamiltonian parameter $r$ computed using a simulator of a noisy
quantum computer. The circles (squares) are computed using qubits with the minimum
(maximum) CNOT error. The solid line is the exact result.}
\label{L4FS}
\end{figure}

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/L4SED.pdf}
\caption{Second derivative of energy per site of the transverse field Ising model with four
sites as a function of the Hamiltonian parameter $r$ computed using a simulator of a noisy
quantum computer. The circles (squares) are computed using qubits with the minimum
(maximum) CNOT error. The solid line is the exact result.}
\label{L4SED}
\end{figure}

On a real quantum device, gate errors and readout errors will degrade the
accuracy of the computation. Readout errors can be mitigated in a fairly
straightforward way so in this study these effects were not
included. The focus here is in quantifying gate error effects and the
mitigation of gate errors is discussed in the next section. The
Qiskit Aer simulator is used since it provides access to
noise models for actual (IBM Q) hardware devices.

The four spin computation is done on two qubits using publicly available
calibration data, so we can choose a device and a qubit pair which has the minimum
CNOT gate error. At the time the computation was done, this was the \texttt{ibmq\_manila}
device. The results of noisy simulation averaged over twenty trials are shown in Figs. \ref{L4FS} and \ref{L4SED}.
For comparison the results of using the qubit pair with the largest CNOT error are
also plotted. For this calculation, where the wave function ansatz contains only a
single two-qubit gate, there is very little sensitivity to gate errors.

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/L6FS.pdf}
\caption{Fidelity susceptibility per site of the transverse field Ising model with six
sites as a function of the Hamiltonian parameter $r$ computed using a simulator of a noisy
quantum computer (circles).  The solid line is the exact result.}
\label{L6FS}
\end{figure}

\begin{figure}[tbh]
\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/L6SED.pdf}
\caption{Second derivative of energy per site of the transverse field Ising model with six
sites as a function of the Hamiltonian parameter $r$ computed using a simulator of a noisy
quantum computer (circles).  The solid line is the exact result.}
\label{L6SED}
\end{figure}

The six spin model, even simplified to run on three qubits, is more challenging. In order to entangle all qubits more CNOT gates are required.
We choose to use CNOT gates only between neighbouring qubits
(Fig. \ref{sixspin})
rather than relying on transpiler-generated swap operations. This is to facilitate the error mitigation study discussed in the next section.
The simulation results using
the \texttt{ibmq\_manila} noise model and coupling map, which specifies connectivity, are shown in Figs. \ref{L6FS} and \ref{L6SED}. Note that no transpiler
optimization was applied in doing these calculations.

\section{Error mitigation and automatic differentiation}
\label{sec:error_mitigation}

\subsection{Richardson extrapolation}

Zero-noise extrapolation is a common means of performing error mitigation
on the results of quantum computations run on hardware. It involves
scaling up the amount of noise in a computation through the systematic addition of noisy gates, followed by extrapolation back to a noise-less value.
Richardson extrapolation (RE) is often used to perform such extrapolation \cite{temme2017error, krebsbach2022optimization}; it is based on
polynomial interpolation using \emph{Lagrange basis polynomials}. Suppose we
measure a set of expectation values at $n + 1$ different points (e.g.,
different noise scale factors), labelled $(x_j, E_j)$, for $j=0,
\ldots, n$. The Lagrange polynomial for this set of points is the lowest-degree
polynomial that passes through all of them exactly. It can be constructed as a linear combination of the expectation values and Lagrange basis polynomials,
\begin{equation}
P_n (x) = \sum_{j=0}^{n} E_j \ell_j(x),
\end{equation}
\noindent where
\begin{equation}
    \ell_j(x) = \prod_{\substack{0 \leq m \leq n, \\ m \neq j}} \frac{x - x_m}{x_j - x_m}.
\end{equation}

\noindent To perform error mitigation with RE, first construct the Lagrange polynomial, and then evaluate it at $x = 0$ to obtain an estimate for the mitigated value. A somewhat simplified expression of the polynomial can be used at this point:
\begin{equation}
    P_n (0) = \sum_{j=0}^{n} E_j \ell_j(0)
    =  \sum_{j=0}^{n} E_j  \prod_{\substack{0 \leq m \leq n, \\ m \neq j}} \frac{ - x_m}{x_j - x_m}
    = \sum_{j=0}^{n} E_j \gamma_j, \label{eq:richardson_extrap}
\end{equation}
as per the notation in Eqs. (3) and (4) of
\cite{krebsbach2022optimization}. From hereon we will denote this value by
$\mathcal{E}$ (mitigated values will in general be represented using a calligraphic font).

In this work, we are interested in performing error mitigation not only on expectation values, but also on gradients of them. Suppose that $\theta$ is the parameter with respect to which we would like to compute a derivative.  It is straightforward to show that for $\mathcal{E} = P_n(0)$,
\begin{equation}
    \frac{\partial \mathcal{E}}{\partial \theta} = \frac{\partial}{\partial \theta} \left( \sum_{j=0}^{n} E_j \gamma_j \right) = \sum_{j=0}^{n} \frac{\partial E_j}{\partial \theta} \gamma_j,
\end{equation}
\noindent and so the gradient of an error-mitigated value is simply the error-mitigated computation of a gradient. Similarly, second derivatives are required, and work the same way:
\begin{equation}
    \frac{\partial^2 \mathcal{E}}{\partial \theta^2} = \sum_{j=0}^{n} \frac{\partial^2 E_j}{\partial \theta^2} \gamma_j. \\
\end{equation}

\noindent Thus, computing gradients and performing error mitigation within the larger context of automatic differentiation will in theory commute.

When we apply RE in practice, each $E_j$ in Eq. (\ref{eq:richardson_extrap})
has been computed independently based on $N_j$ shots. What is obtained, then, is an estimate of $\mathcal{E}$, $\hat{\mathcal{E}}$, computed as a linear combination of estimates $\hat{E}_j$,
\begin{equation}
    \hat{\mathcal{E}} = \sum_{j=0}^{n} \hat{E}_j \gamma_j.
\end{equation}

\noindent The variance of the mitigated value, $\sigma^2$, can be computed in terms of the
variances of the individual estimates. Usually the variance of a sum of terms
involves covariance terms. However, it is reasonable to assume measurements of the expectation values are independent such that all
covariances are 0. This yields
\begin{equation}
    \sigma^2 = \sum_{j=0}^{n} \sigma^2_j \gamma_j^2. \label{eq:var_initial}
  \end{equation}

In the broader context of a quantum algorithm, such as the one considered here,
we are computing the expectation value of not only a
single observable, but of a full Hamiltonian. Let
\begin{equation}
 H = \sum_{i} c_i P_i,
\end{equation}

\noindent be the Hamiltonian of interest, where $P_i$ are $n$-qubit
Pauli operators. The error-mitigated estimate of the energy, $\hat{\mathcal{H}}$,
is given by
\begin{equation}
 \hat{\mathcal{H}} = \sum_{i} c_i \hat{\mathcal{P}}_i.
\end{equation}
\noindent If each estimate $\hat{\mathcal{P}}_i$ has variance $\sigma^2_{P_i}$
according to Eq. \ref{eq:var_initial}, then
the total variance is
\begin{equation}
 \sigma_H^2 = \sum_{i} c^2_i \sigma^2_{P_i} = \sum_i  c_i^2 \sum_{j=0}^n (\sigma^2_{P_i})_{j}\gamma_j^2.
\end{equation}

Unless the number of shots is adjusted, the variance of error-mitigated estimates
is generally higher than those of the unmitigated values, a fact which has been remarked upon in a number of previous works \cite{krebsbach2022optimization, wang2021error, qermit}. If ``too much" error mitigation is applied (e.g., if $n$ is too large), the results may be worse than an unmitigated value unless compensated for. The quality of mitigation can be quantified using an absolute mitigation error  \cite{qermit}, which is reported for these simulations in Appendix B.

The increase in variance is further compounded by the fact that we are computing mitigated values of first and second derivatives, each of which comprises a sum of multiple shifted terms. For example, for the gradient of a single expectation value,
\begin{equation}
    \frac{\partial \hat{\mathcal{E}}}{\partial \theta} = \sum_{j=0}^{n} \frac{\partial \hat{E}_j}{\partial \theta} \gamma_j = \sum_{j=0}^{n} \frac{1}{2} \left(\hat{E}_j (\theta+\pi/2) - \hat{E}_j(\theta - \pi/2) \right) \gamma_j.
\end{equation}

\noindent The variance of the estimate is
\begin{equation}
    (\sigma^2)^\prime = \sum_{j=0}^{n} \frac{1}{4} (\sigma^2_{j+} + \sigma^2_{j-}) \gamma_j^2 \approx \sum_{j=0}^{n} \frac{1}{2} (\sigma^2_{j+}) \gamma_j^2
\end{equation}

\noindent where $\sigma^2_{j\pm}$ represents the variance of the $j^{th}$ term at a shifted value of $\pm \pi/2$, and it is assumed that, due to symmetry, these variances will be equal (but note that they are not necessarily equal to the variance at the unshifted value) \cite{mari2021estimating}. Thus, for the derivative of the expectation value of a full Hamiltonian, a mitigated estimate is obtained with variance
\begin{equation}
 (\sigma_H^2)^\prime = \sum_{i} c^2_i (\sigma^2_{P_i})^\prime = \sum_i  c_i^2 \sum_{j=0}^{n} \frac{1}{2} (\sigma^2_{P_i})_{j+} \gamma_j^2.
\end{equation}

\noindent A similar expression can be derived for Hessians. In both cases there is a dependence on the parameter-shift rule used, and in particular on the number of terms it contains. The increased variance has consequences in our problem due to where they are used: in the solution of the response equations, Eq. (\ref{eq:first_order_response}), whose results are propagated into the computations of second energy derivative and fidelity susceptibility. Even small deviation can result in the system of equations being poorly conditioned, leading to suboptimal solutions.

\subsection{Software implementation}

 While we note that the open-source error-mitigation library Mitiq \cite{mitiq} has a great variety of error mitigation functionality, mitigation was implemented using the differentiable quantum transforms framework in PennyLane \cite{dimatteo2022transforms} in order to preserve differentiability. Transforms are composable
 ``metaprograms" that modify the behaviour of quantum functions and circuits.
In particular, we leverage the concept of a \emph{batch transform}: such transforms take quantum functions as input, and return a set of transformed functions, and a processing function that acts on the results of those functions.  This is shown below in pseudocode:

\begin{lstlisting}[language=python]
transformed_functions, processing_fn = batch_transform(function)
results = execute(transformed_functions)
processing_fn(results)
\end{lstlisting}

Expansion of a Hamiltonian into Pauli terms, computation of gradients and Hessians, and zero-noise extrapolation can all be expressed as batch transforms.  For example, the gradient computation as described in Sect.~\ref{subsec:quantum_gradients} is a batch transform: a quantum circuit is passed as input, and the transform returns two quantum circuits as output (one for each shifted parameter), and a function that combines the executed results (Eq.~(\ref{eq:param_shift})). Furthermore, if the extrapolation procedure is implemented in an autodifferentiable way, the error-mitigated results can be input to subsequent computations without compromising differentiability. A pseudocode implementation of the core function to compute mitigated gradients is shown below to highlight key elements of batch transform functionality such as composition and usage of the reconstruction functions. The full implementation is included in our GitHub repository \cite{ourgithub}.

\begin{lstlisting}[language=python]
import pennylane as qml

# A quantum device
dev = qml.device(...)

# A quantum circuit that computes the expectation of a Hamiltonian
# that we would like to obtain an error-mitigated gradient for
circuit = ...

# Some configuration options for error mitigation; passed to ZNE routine
mitigation_config = {...}

# Split into constituent terms so mitigation is applied to each one separately
ham_circuits, ham_fn = qml.transforms.hamiltonian_expand(circuit, group=False)

# Map gradient computation batch transform over each circuit from the previous step.
grad_circuits, grad_fn = qml.transforms.map_batch_transform(
    qml.gradients.param_shift, ham_circuits
)

# Map each gradient circuit into a set of circuits that computes the error-mitigated
# value. This set of circuits is then executed by the quantum device.
zne_circuits, zne_fn = qml.transforms.map_batch_transform(
    partial(zne, mitigation_config), grad_circuits,
)

execution_results = qml.execute(zne_circuits, dev, gradient_fn=None)

# The processing functions keep track of what should happen to each
# value during the transform, and are used to recombine the results
error_mitigated_expval_gradient = ham_fn(grad_fn(zne_fn(execution_results)))

\end{lstlisting}

The \texttt{mitigation\_config} parameter in the pseudocode described above
provides a great deal of choice in how error mitigation is performed. For instance, how many
noise scale factors are used, what their values are, and how noise gets added
based on them. This work considers two ways of adding noise: unitary folding \cite{giurgicatiron2020digital} and CNOT folding \cite{he2020zeronoise}, each of which are implemented as quantum function transforms \cite{dimatteo2022transforms}. In unitary folding, given an input circuit $U(\theta)$ and an integer noise scale factor $\lambda_j$, the circuit is transformed into $(U(\theta) U^\dagger(\theta))^{(\lambda_j-1)/2} U(\theta) $, i.e., the circuit is run forwards and backwards $(\lambda_j-1)/2$ times. In CNOT folding $(\lambda_j-1)/2$ pairs of CNOTs are added after each existing CNOT.

To verify the mitigation pipeline, it was first applied to computation of ground state energy, the second energy derivative, and the fidelity susceptibility by running the circuit at the optimal variational parameters with no noise, as well as with shot noise. Following this, the simulated device noise model of \texttt{ibmq\_manila} was applied, with the circuits adjusted for the coupling map prior to any transforms being applied. Simulations were run on a desktop machine, but due to the extensive run time of noisy simulations, execution was parallelized over the different values of $r$.

We experimented with sequences of noise scale factors with varying length. Given $\lambda_j = 2j + 1$, simulations were run for sequences of scale factors $[\lambda_0, \ldots, \lambda_n]$ for $n = 0$ (i.e., noisy simulation only) up to $n=4$ (i.e., a max $\lambda_j$ of $9$) for $L=4$, and up to $n=3$ for $L=6$. No notable difference was observed between CNOT folding and unitary folding. In both cases, as the number of folds increases, so does the variance in results; however, due to the nature of the quantities being computed, this variance was strongly amplified, most notably for the cases of max $n=3$ and $n=4$ (up to 3 and 4 folds, respectively).

\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{figs/l4_d2e.pdf}
    \caption{Second energy derivative per site for the 4-spin system. Results are plotted along the horizontal axis for $r \in [0.5, 0.6, \ldots 1.4]$ but offset for readability.  The numbers in the legend indicate which scale factors were used in the mitigation process (e.g., ``1,3,5" means mitigation was performed using 3 points for the extrapolation: the unmodified circuit, adding one fold, and adding two folds). For each mitigation method, 100 independent trials were run; the error bars correspond to the standard deviation of these results. As expected, a greater number of points in the extrapolation leads to higher variance in the potential result.}
    \label{fig:l4_d2e}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{figs/l4_fs.pdf}
    \caption{Fidelity susceptibility per site for the 4-spin system. For each mitigation method, 100 independent trials were run; the error bars correspond to the standard deviation of these results. Greater variation is observed, compared to the second energy derivative of Fig. \ref{fig:l4_d2e}, due to a double dependence on the derivatives of variational parameters with respect to the Hamiltonian parameter $r$.}
    \label{fig:l4_fs}
\end{figure}

Results are shown here for up to the $n=3$ case for $L=4$, and up to $n=2$ for $L=6$; the full set of numerical results is available for perusal on GitHub \cite{ourgithub}.
The $L=4$ ($L=6$) results are shown in Fig. \ref{fig:l4_d2e} and \ref{fig:l4_fs} (Figs. \ref{fig:l6_d2e} and \ref{fig:l6_fs}). To obtain the error-mitigated estimate 8192 shots were used for \emph{each circuit} that was run. The entire computation was run 100 times for each quantity of interest; the error bars on the plots show the mean and standard deviation of the results. We observe the clear increase in variance which occurs from performing error mitigation as the number of points used in the extrapolation is increased. Note that, when results are averaged over all error-mitigated attempts, we do in many cases obtain a mitigated value that improves over the unmitigated one; however this requires a substantial amount of quantum resources. This is further discussed in Appendix B, where plots are provided for the absolute mitigation error over the distribution of results.

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{figs/l6_d2e.pdf}
    \caption{Second energy derivative per site for the 6-spin system. Error bars correspond to the standard deviation of 100 independent trials. Compared to results in the 4-spin case, the variance is higher for both types of folding.}
    \label{fig:l6_d2e}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{figs/l6_fs.pdf}
    \caption{Fidelity susceptibility per site for the 6-spin system. Error bars correspond to the standard deviation of 100 independent trials. Similar to the 4-spin case, the variance in results is higher than that of the second energy derivative. In the case of unitary folding at $r=1.4$, one trial resulted in a markedly large condition number for the matrix inversion involved in solution of the response equations, leading to significantly less accurate results. We found such behaviour to be common for the 6-spin case when pushing beyond the use of three scale factors for both types of folding.}
    \label{fig:l6_fs}
\end{figure}

In some trials, a significant variation was observed, most notably for simulations with a larger number of folds and scale factors. It was found that the reason for this was due to the dependence of both the second energy derivative and the fidelity susceptibility on the quantities $\partial \theta_i / \partial r$. Error mitigation was applied to the computation of each element of the Hessian matrix and gradient of the expectation value of $H_1$ required as per the response equations (Eq. (\ref{eq:first_order_response})). Their solution involves a linear inversion, but when the Hessian is noisy the resulting matrix may be poorly conditioned, leading to very large values of the $\partial \theta_i / \partial r$. These are then used directly in computation of the second energy derivative, and \emph{twice} in the computation of fidelity susceptibility (so its variance is even higher). Obtaining accurate values of these parameters is thus critical.

\section{Summary}

In this work an end-to-end implementation of quantum automatic differentiation for a problem in condensed matter physics was presented. The use of automatic differentiation enables the computation of quantities of interest (second energy derivative and fidelity susceptibility) using effectively a single quantum circuit, without knowledge of the full spectrum of the system. It is hoped that this demonstration will motivate further exploration into how these tools can be applied to other physical problems.

The methods were successfully applied here to small quantum systems with correct results, up to deviations due to statistical and simulated hardware noise. However, the computation of error-mitigated gradients with respect to Hamiltonian parameters was found to be particularly sensitive to such noise. A more general and systematic study should be performed to analyze the theory underlying error-mitigation of energy derivatives, as it has implications for many other problems in areas such as quantum chemistry. This could include testing different methods for computing the matrix inverse when solving the response equations (e.g., through a pseudoinverse / singular value decomposition), or alternative strategies for evaluating the derivatives themselves \cite{obrien2019calculating, Azad2022}. Different (non-ZNE) error mitigation techniques should also be analyzed, and it would be of interest to determine the extent to which the mitigation strategy itself could be optimized, using, e.g., methods developed in \cite{krebsbach2022optimization} or \cite{vaqem}. Finally, as noted in \cite{vaqem}, results of error mitigation applied to noisy simulations may not be representative of those obtained using an actual quantum processor, and so, if feasible, future work should include execution on real devices.

\section{Acknowledgements}

ODM thanks Josh Izaac and Tom Bromley for helpful discussions, and acknowledges funding provided by NSERC grant no. RGPIN-2022-04609. TRIUMF receives federal funding via a contribution agreement with the National Research Council of Canada. We acknowledge the use of the IBM Qiskit framework for
this work. The views expressed are those of the authors,
and do not reflect the official policy or position of IBM
or the IBM Quantum team

\section*{Appendix A: Overlap algorithms}

\begin{figure}[p]
\centering
\includegraphics[scale=0.4,trim={4cm 1cm 6cm 1cm},clip=true,angle=-90]{./figs/hadamard.pdf}
\caption{Circuit for calculating the overlap of the wave functions instantiated by the
unitary operators $U_i$ and $U_f$ using the Hadamard test. The gate R is R${_y}(-\pi/2)$ or
R${_x}(\pi/2)$ for the real and imaginary parts respectively.}
\label{hadam}

\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/hadFS.pdf}
\caption{Fidelity susceptibility of the transverse field Ising model as a function of the
Hamiltonian parameter $r$ for four spins comparing the Hadamard test  (squares)with the overlap
circuit Fig. \ref{uudag} (circles). The solid line is the exact result.}
\label{hadFS}
\end{figure}

The Hadamard test is a standard way to calculate the overlap between two states.
The circuit for this test is shown in generic form in Fig. \ref{hadam}. The advantage of
the Hadamard test is that it can provide separately the real and imaginary parts
of the overlap. The disadvantages, compared to the circuit in Fig. \ref{uudag}
are that it uses an additional qubit and, more importantly, it requires
controlled versions of the unitary operators that create the states. This leads
to a large increase in the number of multi-qubit gates resulting in greater
sensitivity to hardware noise. This is illustrated in Fig. \ref{hadFS} where the fidelity
susceptibility as a function of the Hamiltonian parameter $r$ for 4 spins calculated
using the Hadamard test, and is compared to that calculated with overlap from
Fig. \ref{uudag}. The simulation is done using the \texttt{ibmq\_manila} gate noise model.

\begin{figure}[p]
\centering
\includegraphics[scale=0.4,trim={3cm 1cm 3cm 1cm},clip=true,angle=-90]{./figs/swap.pdf}
\caption{Circuit for calculating the squared overlap of the wave functions instantiated by the
unitary operators $U_i$ and $U_f$ using the swap test. }
\label{swap}

\centering
\includegraphics[scale=0.45,trim={1cm 1cm 1cm 1cm},clip=true,angle=0]{./figs/swapFS.pdf}
\caption{Fidelity susceptibility of the transverse field Ising model as a function of the Hamiltonian parameter $r$ for
four spins using the swap test Fig. \ref{swap} (squares) compared to the overlap circuit
Fig. \ref{uudag} (circles). Results of using alternate versions of the swap test discussed
in \cite{Cincio_2018} are shown using up-triangles  (Ancilla-Based Algorithm, Ref. \cite{Cincio_2018}, Fig. 5)
and down-triangles (Bell-Based Algorithm, Ref. \cite{Cincio_2018}, Fig. 6).
The solid line is the exact result.}
\label{swapFS}
\end{figure}

The swap test is another way to compute the squared overlap of two states. The generic circuit is
shown in Fig. \ref{swap}. The advantage of the swap test over the Hadamard test is that the unitary
operators creating the states are not controlled. However, the number of qubits required to contain
the states is doubled. Instead of controlled unitaries, controlled swaps are needed which is a
disadvantage compared to the overlap circuit in Fig. \ref{uudag}.

In Ref. \cite{Cincio_2018} Cincio \emph{et al.} discuss ways to improve the swap test.
Fig. \ref{swapFS} compares
the fidelity susceptibility as a function of the Hamiltonian parameter
$r$ for 4 spins calculated with the \texttt{ibmq\_manila} gate noise model using
the Ancilla-Based Algorithm (ABA, Ref. \cite{Cincio_2018}, Fig. 5) and the
Bell-Based Algorithm (BBA, Ref. \cite{Cincio_2018}, Fig. 6) with results using a
simple unimproved swap, Fig. \ref{swap}, and the overlap circuit in Fig. \ref{uudag}.
The algorithms from Ref. \cite{Cincio_2018} provide an improvement but, as one might expect,
the circuit Fig. \ref{uudag} which does not introduce controlled gates, other than those contained
in the state-creating unitaries, shows the least sensitivity to hardware noise.

\section*{Appendix B: Additional results}
\label{appendix:more_results}

In this section additional numerical results are presented for the 4- and 6-spin systems. In particular, we show results for mitigated estimation of the ground state energy of both systems, and the absolute mitigation error \cite{qermit} (which, for the unmitigated case, is simply the error).

Note that error mitigation has the best overall performance for the estimation of the ground state energy; this is not surprising given that there is no dependence on the results on the solution of a noisy linear system. As 100 trials were performed at each $r$, there are two errors to plot: the absolute mitigation error for each trial, averaged over all trials; or the absolute mitigation error of the average of all results together. For the ground state energy of the $L=4$ system (Fig. \ref{fig:l4_gse_error}) one sees that a small amount of mitigation generally improves on the results of any individual trial. When considering all the mitigation results together, all mitigation methods see improvement. As each trial is independent, this provides insight into what could be expected in the limit of an infinite number of shots (in this case, obtained from running 819200 shots for each term in the Hamiltonian, for each circuit). For $L=6$ the benefits are more clear. This is not unexpected: given the small size of the 4-spin circuit, the increase in statistical noise that occurs when mitigating is worse than simply using the results as-is given the noise levels in today's processors. On the other hand for the 6-spin case, where there are five CNOTs in the circuit, mitigation is helpful.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l4_gse.pdf}
    \caption{Error-mitigated ground state energies of the 4-spin system. Mean and standard deviation of 100 independent trials are shown for each value of $r$.}
    \label{fig:l4_gse}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l4_gse_error.pdf}
    \caption{Mitigation error for the ground state energy of the 4-spin system. }
    \label{fig:l4_gse_error}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l6_gse.pdf}
    \caption{Error-mitigated ground state energies of the 6-spin system. Mean and standard deviation of 100 independent trials are shown for each value of $r$.}
    \label{fig:l6_gse}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l6_gse_error.pdf}
    \caption{Mitigation error for the ground state energy of the 6-spin system. }
    \label{fig:l6_gse_error}
\end{figure}

For the second energy derivative (Fig. \ref{fig:l4_d2e_error}) and fidelity susceptibility (Fig. \ref{fig:l4_fs_error} in the $L=4$ case), the benefits of mitigation in an average trial are less clear, but again in the limit of an infinite number of shots one can expect an advantage, at least up to a maximum of 3 folds; running with 4 folds produces worse results for both metrics, even for the 4-spin system (this data is provided in the GitHub repository). Running with two different scale factors and then extrapolating back produces the best results for both system sizes; while for the 4-spin case there is again no real improvement noted, there is a clear improvement in the $L=6$ case.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l4_d2e_error.pdf}
    \caption{Mitigation error for the second energy derivative of the 4-spin system. }
    \label{fig:l4_d2e_error}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l4_fs_error.pdf}
    \caption{Mitigation error for the fidelity susceptibility of the 4-spin system. }
    \label{fig:l4_fs_error}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l6_d2e_error.pdf}
    \caption{Mitigation error for the second energy derivative of the 6-spin system. }
    \label{fig:l6_d2e_error}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{figs/l6_fs_error.pdf}
    \caption{Mitigation error for the fidelity susceptibility of the 6-spin system. }
    \label{fig:l6_fs_error}
\end{figure}

\bibliographystyle{h-physrev4}


\section{INTRODUCTION}

Image filtering such as smoothing an image reduces noise, blurred images can be rectified. There are broadly two types of algorithms- linear\cite{kailath1974view,satpathy2010image} and non-linear\cite{addabbo2019classification}. Linear filter is achieved through convolution and Fourier multiplication whereas Non-linear filter cannot be achieved through any of these. Its output is not the linear function of its input thus, its result varies in a non-intuitive manner. Here, the following image shows how the median filter\cite{juneja2009improved}  enhances the images by reducing the noise and smoothing.
\begin{figure}[!htbp]
\begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=.5\textwidth]{median1}
         \caption{}
         \label{fig1a}
     \end{subfigure}
      \begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=.5\textwidth]{median2}
         \caption{}
         \label{fig1b}
     \end{subfigure}
     \caption{(a)Original image , (b) Filtered image}

\end{figure}

In order to do further processing like image segmentation, edge detection etc. noise should be eliminated. Median filter is the most effective non-linear filtering algorithm to detect and remove salt and pepper noise. Median filter retains the edges of the while removing the noise that's why it is used most widely thus, apart from trampling the noise up to 5\% to 60\%, it also preserves the image details. The noisy pixels are evaluated and labelled as noisy pixels and the switching based median filter is applied to other pixels which are not noisy.

Bilateral filter \cite{paris2007gentle} is a type of non-liner filter, it reduces noise by smoothing and preserves edges of the images. It takes weighted sum of the pixels which are nearby of each pixel and replaces the intensity of the pixels with the average of that weighted sum. The neighbourhood pixels are identified through their locations which are relative to the input pixel. Algorithms \cite{chandel2013image} used for linear filtering are: Box blur, Gaussian, bilateral and Hann window.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=10cm]{name}
    \caption{Types of Filtering Algorithms}
\end{figure}

In box blur \cite{zimmer2009blur}, an image with 9x9 pixel values, can be considered having a 3x3 neighbourhood values.  smoothing is achieved by averaging the neighbourhood pixel values of the particular pixel in the output image. In this way the pixels with higher intensity gets converted into a lower valued pixel and vice-versa thus, balancing the image pixels.
Gaussian filter is based on the equation of a Gaussian, it can be used to generate a kernel.
A kernel is a matrix (small) which is a convoluted matrix used in various image filtering techniques for embossing, smoothing, sharpening, blurring etc.
Gaussian blur can be applied for smoothing or filtering operations.

\subsection{Working Example- Box Blur:}

Let's take a 9*9 image as an input for the example. Here in the input image, 90 is the peak value and also some of the holes are 0 which are black point. Considering a 3*3 neighbourhood kernel around the pixel of which average values will be taken. Similarly, average values around all the pixels will be calculated by moving the kernel frame all over the input image to be able to generate newer versions of the output image. For example, if sum around a particular pixel is 360 then the value in corresponding pixel of output image will be 40 since, there are 9 pixels around the input pixel's kernel. Thus, in this manner all the values of the output matrix will be filled.

\begin{figure}[!htbp]
\begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{box-2}
         \caption{}
         \label{fig2a}
     \end{subfigure}
      \begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{gray}
         \caption{}
         \label{fig2b}
     \end{subfigure}
     \caption{(a)Input and  Output image , (b) Input and Output image in shades of gray}

\end{figure}

It can be observed from the image that the place where zero was surrounded by 90s, in the output image is replaced by a much higher value and thus smoothing intensity of the pixel. Also 90 which was surrounded by 0's is now reduced to 10 producing a smoothen output.

\section{Image filtering algorithms }
\subsection{Non-linear}
\subsubsection{Median Filter \cite{juneja2009improved}}
Before proceeding for filtering the image, noise should be detected properly in order to rectify the noise and smoother the image. Following assumptions are made for error detection i) in an image which is noise-free edges separate the smoothly varying areas, ii) the pixel which is noisy in the image has uneven intensity i.e. either very low or very high in comparison of the neighborhood pixels.
During the impulse noise detection procedure, two image sequences are generated. One is the series of gray-scale images,  $x^{(0)}_{(i,j)} , x^{(1)}_{(i,j)} ,\dots,x^{(n)}_{(i,j)}$.
Here, the first image $x^{(0)}_{(i,j)}$ is noisy image, ${(i,j)}$ is the position of pixels in the image where $1 \leq i \leq M, 1 \leq j \leq N$. Here, N and M represent direction of pixels in vertical and horizontal direction respectively and $x^{(n)}_{(i,j)}$ is the image after nth iteration.

Noise detection procedure:
\begin{enumerate}
  \item	Assume a window $(2W+1) x (2W+1)$ around $x^{(n-1)}_{(i,j)}$  i.e.
      \[ x^{(i+k,j+l)}_{(n-1)}\]  where $-W<=k<=W, -W<=l<=W$ and $W>=1$.

  \item	Finding windows, median value $m^{(i,j)}_{(n-1)}$  ,
        \[m^{(i,j)}_{(n-1)}   = median (x^{(i+k,j+l)}_{(n-1)})\]
  \item	The absolute difference of $x^{(n-1)}_{(i,j)}$   and  $m^{(i,j)}_{(n-1)}$  ,
\[f^{(i,j)}_{(n)} =
\begin{cases}
      {f ^{(i,j)}_{(n-1)},   if | x^{(n-1)}_{(i,j)}  -   m^{(i,j)}_{(n-1)} | <T}\\
      {1, otherwise}
    \end{cases} \]
\item	If noise is detected in (i,j)th  then the value of x(n)(I,j)   will be changed as,
\[x^{(n)}_{(i,j)} =
\begin{cases}
     m^{(i,j)}_{(n-1)}  ,    if f^{(i,j)}_{(n)} \neq f^{(i,j)}_{(n-1)} \\
     x(n-1)(i,j)  ,    if f^{(i,j)}_{(n)} \neq f^{(i,j)}_{(n-1)}
    \end{cases} \]

\end{enumerate}
All these steps needed to be repeated for p times. Where p can be 2,3,4…. and so on.
After the pth iteration, $x^{(p)}_{(i,j)}$  and  $f^{(p)}_{(i,j)}$    are two images formed. But for noise filtering process, the binary flag image $f^{(p)}_{(i,j)}$ is only required. Sun and Neuvo introduced this median based noise-based detection and later Wang and Zhang (PSM) modified it in a progressive manner. This method or approach is distinct from PSM as in this method various adaptive filtering approach as it focuses and improves the filter's filtering performance.

\subsubsection{Bilateral Filtering}
   It is a simple and non-iterative method. Bilateral filtering\cite{weiss2006fast} preserves the edges while smoothing the noisy image by virtue of combination of non-linear nearby imagee value. Based on photometric similarity and geometric closeness. In this method, no phantom colours along the edges of the output images are produced.
    Bilateral filter is similar to Gaussian convolution as it is the average of the pixels but the bilateral filter considers the variations in intensities thus, preserves the edges of the image. Bilateral filtering can be justified as when the two pixels are near to each other and they are similar in terms of the photometric range along with the occupancies of the spatial locations which are nearby.
\subsection{Gaussian \textit{vs.} Bilateral}
The heat equation or isotropic diffusion has a Gaussian correlation. It applies uniformly to the picture, thereby encompassing the borders and image content. The Gaussian and Bilateral filters are nearly identical, except that the Gaussian inflects through a modulated function that determines the distance between the center and neighboring pixels. There is no filtering at the object's boundaries with Gaussian filtering, and the filter is applied evenly throughout the picture. On the other hand, a bilateral filter prevents blurring between objects while yet being successful in evenly eliminating noise and so giving a superior outcome.

\section{Applications}
\begin{enumerate}

\item \textit{Image Enhancement:} Image filtering helps in enhancing the quality of the image pixels. It comprises of the processes of altering the images in terms of pixel values (blurring, smoothing etc.) irrespective of the image type i.e. photochemical, digital or even illustration, this processing is known as photo retouching in case of the analog image enhancing. Various softwares are also used which encompasses three dimensional modelers, vector and raster graphic editors. These softwares are primary tools through which an image can be enhanced and transformed as desired\cite{jain1989fundamentals}. \\
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=8cm]{imen}
    \caption{Noisy Image and Filtered Image}
\end{figure}

\item \textit{Denoising:} A ten function is used as a range weight, while a square box function is used as a spatial weight. The bilateral filter, unlike Gaussian blur, delivers sharper and superior results by focusing on maintaining the object outlines. It has many applications in fields like movie restoration and medical imaging.
\begin{figure}[!htbp]
\begin{subfigure}[b]{0.33\textwidth}
         \centering
         \includegraphics[width=.8\textwidth]{den-input}
         \caption{}
         \label{fig3a}
     \end{subfigure}%
      \begin{subfigure}[b]{0.33\textwidth}
         \centering
         \includegraphics[width=.8\textwidth]{den-noisy}
         \caption{}
         \label{fig3b}
     \end{subfigure}%
     \begin{subfigure}[b]{0.33\textwidth}
         \centering
         \includegraphics[width=.8\textwidth]{den-noisy}
         \caption{}
         \label{fig3c}
     \end{subfigure}
     \caption{(a)Input image , (b) Noisy image, (c) Surface blur filter}

\end{figure}

\item \textit{Image Compression:}
Compression can lead to partial data loss i.e., it leads to lossy compression. But in the field of medical image processing \cite{schalkoff1989digital}, clipart, comics, and technical drawings, lossless image data is preferred. Compression artifacts can be introduced during lossy compression methods when short-bit rates are used. The compression, which is lossy but doesn't create unobtrusiveness, can be called visually lossless.

\item \textit{Texture and illumination separation:} Texture and illumination separation is done in image based modelling. The uneven texture in images is removed by the bilateral filter successfully and also preserves the cessation which are caused by geometrical and brightness changes. This method is inspired from the fact that in an image, variations in brightness are much higher as compared to that in texture patterns\cite{katsaggelos2012digital,satpathy2010adaptive,satpathy2010image}.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=8cm]{restore}
    \caption{Noisy Image and Restored Image}
\end{figure}

\item \textit{Character Recognition:} It is a common procedure to digitize printed compositions with the ultimate goal that they can be altered, looked at and stored even more firmly in machine procedures, For instance, webbased machine interpretation, content todiscourse, extraction of important data and content mining. OCR is a knoll of insight, example, and PC vision research.Early shapes needed to be computerized with each character's images and worked on one text style at a time. There are some appealing techniques for copying engineered yields, particularly after the first filtered sheet, including sections, images and other unprinted segments\cite{mollah2011design}.

\item \textit{Signature Verification:} A computerized mark is a science scheme to speak to an advanced correspondence's authenticity.A lawful computerized mark bears the beneficiary's motivation to think that a perceived sender made the message to such an extent that the sender could not refuse to send the message with nonrepudiation and confirmation and that the message was not changed in exchange.Advanced marks are usually used to define programming, budget correspondence, and in other circumstances where impersonation or alteration is essential.

\item \textit{Biometrics:}
To depict the last class of biometrics, some researchers have instituted the word behavior measurements \cite{nagwanshi2012b,nagwanshi2013d,nagwanshi2018s,nagwanshi2019c}. The enmity of Fingerprint Verification (FVC) is an intercontinental challenge focused on evaluating the programming of distinctive mark confirmation. A subgroup of distinctive mark impressions was provided to enlisted people with different sensors allowing the parameters of distinct calculations to be changed. Individuals were spoken to offer enlist and match executable records of their calculations ; the evaluation was carried out at the courtesies of the coordinators using the submitted compelling documents on an appropriate database acquired as the preparation set with the indistinguishable sensors \cite{nagwanshi2020estimation,basheer2021fesd,nagwanshi2012introduction}.

\item \textit{Object Recognition:} Item discovery is a PC innovation recognized with PC vision and image preparation that manages to see depictions of semantic objects of a class, such as individuals, constructions or cars in computerized recordings and images. Well-inquired regions of discovery of article integrate face location and recognition of walkers. Item recognition is guaranteed in many areas of PC vision, similar to picture recovery and video observation.

\begin{figure}[!htbp]
\begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{objrec}
         \caption{}
         \label{fig2a}
     \end{subfigure}
      \begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{facedet}
         \caption{}
         \label{fig2b}
     \end{subfigure}
     \caption{(a)Images in Object Recognition , (b) Face Images.}
\end{figure}

\item \textit{Face Detection:}
Face identification is subject to PC technology that constructs nonsensical (sophisticated) photos of the dimensions and regions of human faces \cite{ahmad2013image,satpathy2010image}.
\item \textit{Medical Applications:} In particular, we attach importance to restorative imagery as the critical use of image handling. Despite the fact that imaging can be performed for medicinal subtleties of expelled organs and tissues, such occasions are not usually referred to as therapeutic imaging, but to some extent they are a piece of pathology.As a reprimand and in its broadest context, it is a piece of natural imaging and includes radiology that utilizes Xbeam radiography imaging advances, appealing reverberation imaging, ultrasonographic or ultrasound therapeutic material, endoscopy, elastography, imaging, thermography, restorative and nuclear prescription helpful imaging techniques as positron outflow tomography\cite{Raja21,mahmood19,patra20}.
\end{enumerate}

\section{Conclusions}
Image filtering has a wide range of applications including edge detection, removing noise, sharpening and smoothing. Image filtering has brought an evolutionary change in the field of image processing.  A filter can be characterized as a kernel, which can be defined as a small array applied to every pixel and their adjacent neighbors in the given image. Here, in this paper we have briefly explained about the image filtering and compared various algorithms. We have also briefly discussed about its various applications such as image enhancement, image compression and restoration, denoising and in medical field. When it is needed to reduce noise while preserving the peaks and edges of the image, median and bilateral filter is recommended. But, if only the peaks of the image are needed to be preserved and effect on the edge is not to be bothered then gaussian filter is sufficient as it requires less computation comparatively. Whereas if only noise reduction is concerned and peaks and edges are not in consideration, box filter will suffice. It is recommended keeping the kernel size less than or equal to 7 for any algorithm for optimal computational cost.

\section{INTRODUCTION}

Image filtering such as smoothing an image reduces noise, blurred images can be rectified. There are broadly two types of algorithms- linear\cite{kailath1974view,satpathy2010image} and non-linear\cite{addabbo2019classification}. Linear filter is achieved through convolution and Fourier multiplication whereas Non-linear filter cannot be achieved through any of these. Its output is not the linear function of its input thus, its result varies in a non-intuitive manner. Here, the following image shows how the median filter\cite{juneja2009improved}  enhances the images by reducing the noise and smoothing.
\begin{figure}[!htbp]
\begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=.5\textwidth]{median1}
         \caption{}
         \label{fig1a}
     \end{subfigure}
      \begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=.5\textwidth]{median2}
         \caption{}
         \label{fig1b}
     \end{subfigure}
     \caption{(a)Original image , (b) Filtered image}

\end{figure}

In order to do further processing like image segmentation, edge detection etc. noise should be eliminated. Median filter is the most effective non-linear filtering algorithm to detect and remove salt and pepper noise. Median filter retains the edges of the while removing the noise that's why it is used most widely thus, apart from trampling the noise up to 5\% to 60\%, it also preserves the image details. The noisy pixels are evaluated and labelled as noisy pixels and the switching based median filter is applied to other pixels which are not noisy.

Bilateral filter \cite{paris2007gentle} is a type of non-liner filter, it reduces noise by smoothing and preserves edges of the images. It takes weighted sum of the pixels which are nearby of each pixel and replaces the intensity of the pixels with the average of that weighted sum. The neighbourhood pixels are identified through their locations which are relative to the input pixel. Algorithms \cite{chandel2013image} used for linear filtering are: Box blur, Gaussian, bilateral and Hann window.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=10cm]{name}
    \caption{Types of Filtering Algorithms}
\end{figure}

In box blur \cite{zimmer2009blur}, an image with 9x9 pixel values, can be considered having a 3x3 neighbourhood values.  smoothing is achieved by averaging the neighbourhood pixel values of the particular pixel in the output image. In this way the pixels with higher intensity gets converted into a lower valued pixel and vice-versa thus, balancing the image pixels.
Gaussian filter is based on the equation of a Gaussian, it can be used to generate a kernel.
A kernel is a matrix (small) which is a convoluted matrix used in various image filtering techniques for embossing, smoothing, sharpening, blurring etc.
Gaussian blur can be applied for smoothing or filtering operations.

\subsection{Working Example- Box Blur:}

Let's take a 9*9 image as an input for the example. Here in the input image, 90 is the peak value and also some of the holes are 0 which are black point. Considering a 3*3 neighbourhood kernel around the pixel of which average values will be taken. Similarly, average values around all the pixels will be calculated by moving the kernel frame all over the input image to be able to generate newer versions of the output image. For example, if sum around a particular pixel is 360 then the value in corresponding pixel of output image will be 40 since, there are 9 pixels around the input pixel's kernel. Thus, in this manner all the values of the output matrix will be filled.

\begin{figure}[!htbp]
\begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{box-2}
         \caption{}
         \label{fig2a}
     \end{subfigure}
      \begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{gray}
         \caption{}
         \label{fig2b}
     \end{subfigure}
     \caption{(a)Input and  Output image , (b) Input and Output image in shades of gray}

\end{figure}

It can be observed from the image that the place where zero was surrounded by 90s, in the output image is replaced by a much higher value and thus smoothing intensity of the pixel. Also 90 which was surrounded by 0's is now reduced to 10 producing a smoothen output.

\section{Image filtering algorithms }
\subsection{Non-linear}
\subsubsection{Median Filter \cite{juneja2009improved}}
Before proceeding for filtering the image, noise should be detected properly in order to rectify the noise and smoother the image. Following assumptions are made for error detection i) in an image which is noise-free edges separate the smoothly varying areas, ii) the pixel which is noisy in the image has uneven intensity i.e. either very low or very high in comparison of the neighborhood pixels.
During the impulse noise detection procedure, two image sequences are generated. One is the series of gray-scale images,  $x^{(0)}_{(i,j)} , x^{(1)}_{(i,j)} ,\dots,x^{(n)}_{(i,j)}$.
Here, the first image $x^{(0)}_{(i,j)}$ is noisy image, ${(i,j)}$ is the position of pixels in the image where $1 \leq i \leq M, 1 \leq j \leq N$. Here, N and M represent direction of pixels in vertical and horizontal direction respectively and $x^{(n)}_{(i,j)}$ is the image after nth iteration.

Noise detection procedure:
\begin{enumerate}
  \item	Assume a window $(2W+1) x (2W+1)$ around $x^{(n-1)}_{(i,j)}$  i.e.
      \[ x^{(i+k,j+l)}_{(n-1)}\]  where $-W<=k<=W, -W<=l<=W$ and $W>=1$.

  \item	Finding windows, median value $m^{(i,j)}_{(n-1)}$  ,
        \[m^{(i,j)}_{(n-1)}   = median (x^{(i+k,j+l)}_{(n-1)})\]
  \item	The absolute difference of $x^{(n-1)}_{(i,j)}$   and  $m^{(i,j)}_{(n-1)}$  ,
\[f^{(i,j)}_{(n)} =
\begin{cases}
      {f ^{(i,j)}_{(n-1)},   if | x^{(n-1)}_{(i,j)}  -   m^{(i,j)}_{(n-1)} | <T}\\
      {1, otherwise}
    \end{cases} \]
\item	If noise is detected in (i,j)th  then the value of x(n)(I,j)   will be changed as,
\[x^{(n)}_{(i,j)} =
\begin{cases}
     m^{(i,j)}_{(n-1)}  ,    if f^{(i,j)}_{(n)} \neq f^{(i,j)}_{(n-1)} \\
     x(n-1)(i,j)  ,    if f^{(i,j)}_{(n)} \neq f^{(i,j)}_{(n-1)}
    \end{cases} \]

\end{enumerate}
All these steps needed to be repeated for p times. Where p can be 2,3,4…. and so on.
After the pth iteration, $x^{(p)}_{(i,j)}$  and  $f^{(p)}_{(i,j)}$    are two images formed. But for noise filtering process, the binary flag image $f^{(p)}_{(i,j)}$ is only required. Sun and Neuvo introduced this median based noise-based detection and later Wang and Zhang (PSM) modified it in a progressive manner. This method or approach is distinct from PSM as in this method various adaptive filtering approach as it focuses and improves the filter's filtering performance.

\subsubsection{Bilateral Filtering}
   It is a simple and non-iterative method. Bilateral filtering\cite{weiss2006fast} preserves the edges while smoothing the noisy image by virtue of combination of non-linear nearby imagee value. Based on photometric similarity and geometric closeness. In this method, no phantom colours along the edges of the output images are produced.
    Bilateral filter is similar to Gaussian convolution as it is the average of the pixels but the bilateral filter considers the variations in intensities thus, preserves the edges of the image. Bilateral filtering can be justified as when the two pixels are near to each other and they are similar in terms of the photometric range along with the occupancies of the spatial locations which are nearby.
\subsection{Gaussian \textit{vs.} Bilateral}
The heat equation or isotropic diffusion has a Gaussian correlation. It applies uniformly to the picture, thereby encompassing the borders and image content. The Gaussian and Bilateral filters are nearly identical, except that the Gaussian inflects through a modulated function that determines the distance between the center and neighboring pixels. There is no filtering at the object's boundaries with Gaussian filtering, and the filter is applied evenly throughout the picture. On the other hand, a bilateral filter prevents blurring between objects while yet being successful in evenly eliminating noise and so giving a superior outcome.

\section{Applications}
\begin{enumerate}

\item \textit{Image Enhancement:} Image filtering helps in enhancing the quality of the image pixels. It comprises of the processes of altering the images in terms of pixel values (blurring, smoothing etc.) irrespective of the image type i.e. photochemical, digital or even illustration, this processing is known as photo retouching in case of the analog image enhancing. Various softwares are also used which encompasses three dimensional modelers, vector and raster graphic editors. These softwares are primary tools through which an image can be enhanced and transformed as desired\cite{jain1989fundamentals}. \\
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=8cm]{imen}
    \caption{Noisy Image and Filtered Image}
\end{figure}

\item \textit{Denoising:} A ten function is used as a range weight, while a square box function is used as a spatial weight. The bilateral filter, unlike Gaussian blur, delivers sharper and superior results by focusing on maintaining the object outlines. It has many applications in fields like movie restoration and medical imaging.
\begin{figure}[!htbp]
\begin{subfigure}[b]{0.33\textwidth}
         \centering
         \includegraphics[width=.8\textwidth]{den-input}
         \caption{}
         \label{fig3a}
     \end{subfigure}%
      \begin{subfigure}[b]{0.33\textwidth}
         \centering
         \includegraphics[width=.8\textwidth]{den-noisy}
         \caption{}
         \label{fig3b}
     \end{subfigure}%
     \begin{subfigure}[b]{0.33\textwidth}
         \centering
         \includegraphics[width=.8\textwidth]{den-noisy}
         \caption{}
         \label{fig3c}
     \end{subfigure}
     \caption{(a)Input image , (b) Noisy image, (c) Surface blur filter}

\end{figure}

\item \textit{Image Compression:}
Compression can lead to partial data loss i.e., it leads to lossy compression. But in the field of medical image processing \cite{schalkoff1989digital}, clipart, comics, and technical drawings, lossless image data is preferred. Compression artifacts can be introduced during lossy compression methods when short-bit rates are used. The compression, which is lossy but doesn't create unobtrusiveness, can be called visually lossless.

\item \textit{Texture and illumination separation:} Texture and illumination separation is done in image based modelling. The uneven texture in images is removed by the bilateral filter successfully and also preserves the cessation which are caused by geometrical and brightness changes. This method is inspired from the fact that in an image, variations in brightness are much higher as compared to that in texture patterns\cite{katsaggelos2012digital,satpathy2010adaptive,satpathy2010image}.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=8cm]{restore}
    \caption{Noisy Image and Restored Image}
\end{figure}

\item \textit{Character Recognition:} It is a common procedure to digitize printed compositions with the ultimate goal that they can be altered, looked at and stored even more firmly in machine procedures, For instance, webbased machine interpretation, content todiscourse, extraction of important data and content mining. OCR is a knoll of insight, example, and PC vision research.Early shapes needed to be computerized with each character's images and worked on one text style at a time. There are some appealing techniques for copying engineered yields, particularly after the first filtered sheet, including sections, images and other unprinted segments\cite{mollah2011design}.

\item \textit{Signature Verification:} A computerized mark is a science scheme to speak to an advanced correspondence's authenticity.A lawful computerized mark bears the beneficiary's motivation to think that a perceived sender made the message to such an extent that the sender could not refuse to send the message with nonrepudiation and confirmation and that the message was not changed in exchange.Advanced marks are usually used to define programming, budget correspondence, and in other circumstances where impersonation or alteration is essential.

\item \textit{Biometrics:}
To depict the last class of biometrics, some researchers have instituted the word behavior measurements \cite{nagwanshi2012b,nagwanshi2013d,nagwanshi2018s,nagwanshi2019c}. The enmity of Fingerprint Verification (FVC) is an intercontinental challenge focused on evaluating the programming of distinctive mark confirmation. A subgroup of distinctive mark impressions was provided to enlisted people with different sensors allowing the parameters of distinct calculations to be changed. Individuals were spoken to offer enlist and match executable records of their calculations ; the evaluation was carried out at the courtesies of the coordinators using the submitted compelling documents on an appropriate database acquired as the preparation set with the indistinguishable sensors \cite{nagwanshi2020estimation,basheer2021fesd,nagwanshi2012introduction}.

\item \textit{Object Recognition:} Item discovery is a PC innovation recognized with PC vision and image preparation that manages to see depictions of semantic objects of a class, such as individuals, constructions or cars in computerized recordings and images. Well-inquired regions of discovery of article integrate face location and recognition of walkers. Item recognition is guaranteed in many areas of PC vision, similar to picture recovery and video observation.

\begin{figure}[!htbp]
\begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{objrec}
         \caption{}
         \label{fig2a}
     \end{subfigure}
      \begin{subfigure}[b]{0.5\textwidth}
         \centering
         \includegraphics[width=\textwidth]{facedet}
         \caption{}
         \label{fig2b}
     \end{subfigure}
     \caption{(a)Images in Object Recognition , (b) Face Images.}
\end{figure}

\item \textit{Face Detection:}
Face identification is subject to PC technology that constructs nonsensical (sophisticated) photos of the dimensions and regions of human faces \cite{ahmad2013image,satpathy2010image}.
\item \textit{Medical Applications:} In particular, we attach importance to restorative imagery as the critical use of image handling. Despite the fact that imaging can be performed for medicinal subtleties of expelled organs and tissues, such occasions are not usually referred to as therapeutic imaging, but to some extent they are a piece of pathology.As a reprimand and in its broadest context, it is a piece of natural imaging and includes radiology that utilizes Xbeam radiography imaging advances, appealing reverberation imaging, ultrasonographic or ultrasound therapeutic material, endoscopy, elastography, imaging, thermography, restorative and nuclear prescription helpful imaging techniques as positron outflow tomography\cite{Raja21,mahmood19,patra20}.
\end{enumerate}

\section{Conclusions}
Image filtering has a wide range of applications including edge detection, removing noise, sharpening and smoothing. Image filtering has brought an evolutionary change in the field of image processing.  A filter can be characterized as a kernel, which can be defined as a small array applied to every pixel and their adjacent neighbors in the given image. Here, in this paper we have briefly explained about the image filtering and compared various algorithms. We have also briefly discussed about its various applications such as image enhancement, image compression and restoration, denoising and in medical field. When it is needed to reduce noise while preserving the peaks and edges of the image, median and bilateral filter is recommended. But, if only the peaks of the image are needed to be preserved and effect on the edge is not to be bothered then gaussian filter is sufficient as it requires less computation comparatively. Whereas if only noise reduction is concerned and peaks and edges are not in consideration, box filter will suffice. It is recommended keeping the kernel size less than or equal to 7 for any algorithm for optimal computational cost.