# CT Material Decomposition using
Spectral Diffusion Posterior Sampling

Xiao Jiang, Grace J. Gang, and J. Webster Stayman Xiao Jiang and J. Webster
Stayman are with the Department of Biomedical Engineering, Johns Hopkins
University, MD 21205 USA (e-mail<EMAIL_ADDRESS>web.stayman@jhu.edu).
Grace J. Gang is with the Department of Radiology, University of Pennsylvania.
PA 19104 USA (e-mail: grace.j.gang@pennmedicine.upenn.edu).

###### Abstract

In this work, we introduce a new deep learning approach based on diffusion
posterior sampling (DPS) to perform material decomposition from spectral CT
measurements. This approach combines sophisticated prior knowledge from
unsupervised training with a rigorous physical model of the measurements. A
faster and more stable variant is proposed that uses a “jumpstarted” process
to reduce the number of time steps required in the reverse process and a
gradient approximation to reduce the computational cost. Performance is
investigated for two spectral CT systems: dual-kVp and dual-layer detector CT.
On both systems, DPS achieves high Structure Similarity Index Metric
Measure(SSIM) with only $10\%$ of iterations as used in the model-based
material decomposition(MBMD). Jumpstarted DPS (JSDPS) further reduces
computational time by over $85\%$ and achieves the highest accuracy, the
lowest uncertainty, and the lowest computational costs compared to classic DPS
and MBMD. The results demonstrate the potential of JSDPS for providing
relatively fast and accurate material decomposition based on spectral CT data.

###### Index Terms:

Spectral CT, material decomposition, diffusion model.

##  I Introduction

Spectral CT has enabled a number of applications including density estimation,
virtual monochromatic imaging, and contrast agent enhancement. Many of these
applications rely on material decomposition, which reconstructs basis material
density maps from the spectral projections[1]. Material decomposition is an
ill-conditioned nonlinear inverse problem without an explicit solution.
Existing material decomposition algorithms can be categorized into three main
types: analytical decomposition[2], iterative/model-based decomposition[3],
and learning-based decomposition[4].

Deep learning methods have been extensively used in medical image formation
including spectral CT[4, 5, 6, 7]. Such approaches leverage prior knowledge
learned from large datasets and can generally surpass classic approaches in
terms of image quality. However, many deep learning methods do not directly
leverage a physical model of the measurements, raising concerns about
robustness, network hallucinations, etc. Some methods integrate a physical
model as part of network training to improve the data consistency[8, 9].
However, these trained networks are tailored to specific system models,
requiring network retraining for system changes in protocol, technique,
device, etc.

Diffusion Posterior Sampling (DPS)[10] has established a new framework to
integrate a learned prior and physics model. Specifically, DPS starts with
unsupervised training of a score-based generative model (SGM) to capture the
target domain distribution, then application of a reverse process to estimate
image parameters - alternating between SGM reverse sampling to drive the image
towards the target distribution and model-based updates to improve the data
consistency with measurements. As a result, the final output adheres to both
the prior distribution and the actual measurements. A major advantage of DPS
is that its network training is not specific to any physical model, allowing
for application across different imaging systems. Recently, we have applied
DPS to nonlinear CT reconstruction[11], demonstrating its effectiveness for
both low-mA and sparse-view single-energy CT reconstruction.

In this work, the original DPS framework is expanded to Spectral DPS (SDPS)
specifically for material decomposition in spectral CT. To further enhance
this method, we propose a “jumpstarted” sampling strategy with gradient
approximation, which significantly stabilizes the sampling process and reduces
computational costs. The performance of both SDPS and its jumpstarted variant
(JSDPS) are demonstrated in dual-kVp[12] and dual-layer CT[13].

##  II Methodology

###  II-A Spectral Deep Posterior Sampling

####  II-A1 Score-based Generative Model (SGM) for Inverse Problems

A SGM[14] defines a forward process which continuously perturbs the target-
domain sample with time-dependent noise. The corresponding reverse process
generates new samples by inverting the perturbation process. Both forward and
reverse processes are described by a stochastic differential equation (SDE).
Specifically, the SDEs for a Denoising Diffusion Probabilistic Model
(DDPM)[15] have the following form at time $t$:

$Forward:\mathrm{d}\textbf{x}=-\frac{\beta_{t}}{2}\textbf{x}\mathrm{d}t+\sqrt{\beta_{t}}\mathrm{d}\textbf{w}$
(1a)
$Reverse:\mathrm{d}\textbf{x}=[-\frac{\beta_{t}}{2}\textbf{x}-\beta_{t}\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{x}_{t}})]\mathrm{d}t+\sqrt{\beta_{t}}\mathrm{d}\textbf{w}$
(1b)

where $\beta_{t}$ is the time-dependent noise variance and
$\mathrm{d}\textbf{w}$ is the standard Wiener process. The unknown score
function $\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{x}_{t}})$ is
approximated by a deep neural network $\textbf{s}_{\theta}(\textbf{x}_{t},t)$.
A SGM is also capable of posterior sampling from conditional distribution
$p(\textbf{x}|\textbf{y})$:

$\mathrm{d}\textbf{x}=[-\frac{\beta_{t}}{2}\textbf{x}-\beta_{t}\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{x}_{t}|\textbf{y}})]\mathrm{d}t+\sqrt{\beta_{t}}\mathrm{d}\textbf{w}$
(2)

In the context of CT, y and x are the projection measurements and the image
volume, respectively. Leveraging Bayes rule $p(\textbf{x}|\textbf{y})\propto
p(\textbf{x})p(\textbf{y}|\textbf{x})$, the reverse sampling (2) may be
reformulated as:

$\begin{aligned}
\mathrm{d}\textbf{x}=&[-\frac{\beta_{t}}{2}\textbf{x}-\beta_{t}\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{x}_{t}})]\mathrm{d}t+\sqrt{\beta_{t}}\mathrm{d}\textbf{w}\\\
&-\beta_{t}\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{y}|\textbf{x}_{t}})\mathrm{d}t\end{aligned}\
$ (3)

The conditional distribution $p_{t}({\textbf{y}|\textbf{x}_{t}})$ can be
approximated[10] as

$\ p_{t}(\textbf{y}|\textbf{x}_{t})\approx
p_{t}(\textbf{y}|\hat{\textbf{x}}_{0}),\text{where}\
\hat{\textbf{x}}_{0}=\frac{1}{\sqrt{\overline{\alpha}_{t}}}(\textbf{x}_{t}+(1-\overline{\alpha}_{t}))\textbf{s}_{\theta}(\textbf{x}_{t},t)$
(4)

Based on this approximation, the deep posterior sampling process could be
expressed as[11]:

$\displaystyle\mathrm{d}\textbf{x}=$
$\displaystyle[-\frac{\beta_{t}}{2}\textbf{x}-\beta_{t}\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{x}_{t}})]\mathrm{d}t+\sqrt{\beta_{t}}\mathrm{d}\textbf{w}$
(5)
$\displaystyle-\beta_{t}\nabla_{\hat{\textbf{x}}_{0}}\mathrm{log}p_{t}({\textbf{y}|\hat{\textbf{x}}_{0}})\nabla_{\textbf{x}_{t}}{\hat{\textbf{x}}_{0}}\mathrm{d}t$

Note that the first two term is exactly same as the unconditional DDPM
sampling (1b), and $p_{t}(\textbf{y}|\hat{\textbf{x}}_{0})$ is the likelihood
term for the measurements. Therefore, the reverse sampling (5) integrates the
prior information captured by DDPM and the physical measurement model provided
by the likelihood function.

The forward process is discretized into $T$ time steps [15]:

$\displaystyle\textbf{x}_{t}=\sqrt{\bar{\alpha}_{t}}\textbf{x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\boldsymbol{\epsilon}$
(6) where $\displaystyle\ \bar{\alpha}_{t}=\prod_{i=1}^{t}(1-\beta_{t}),\
\boldsymbol{\epsilon}\sim\mathcal{N}(0,\boldsymbol{I}),\ t=1,2...T$

Since $\boldsymbol{\epsilon}=-\sqrt{1-\bar{\alpha}_{t}}\
\mathrm{log}p_{t}({\textbf{x}_{t}})$, DDPM trains a network
$\boldsymbol{\epsilon}_{\theta}(\textbf{x},t)$ to predict the noise
$\boldsymbol{\epsilon}$. That is, network parameters, $\theta$, are estimated

$\theta^{*}=\text{argmin}\
\mathbb{E}_{\textbf{x}_{0}}\mathbb{E}_{\boldsymbol{\epsilon},t}\|\boldsymbol{\epsilon}_{\theta}(\textbf{x}_{t},t)-\boldsymbol{\epsilon}\|_{2}^{2}$
(7)

Using the trained network and physical model, the image $\textbf{x}_{0}$ may
be reconstructed by solving the reverse SDE[15, 16].

####  II-A2 Diffusion Posterior Sampling for Material Decomposition

This work aims to apply the DPS for the material decomposition problem, which
estimates basis material densities from the spectral CT measurements[1]. We
adopt a general spectral CT model[13] where measurements y are assumed to
follow a multi-variate Gaussian distribution with covariance K:
$\textbf{y}\sim\mathcal{N}(\overline{\textbf{y}},\textbf{K})$ and mean:

$\overline{\textbf{y}}=\textbf{BS}\exp({-\textbf{QAx}})$ (8)

where x are densities for each basis volume, Q are mass attenuation
coefficients for each basis, A represents (channel-specific) projection
matrices, S are channel-specific spectral response, and B captures overall
gain effects. Substituting Eq.(8) into Eq.(5), we have the spectral DPS (SDPS)
for material decomposition:

$\displaystyle\mathrm{d}\textbf{x}=$
$\displaystyle[-\frac{\beta_{t}}{2}\textbf{x}-\beta_{t}\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{x}_{t}})]\mathrm{d}t+\sqrt{\beta_{t}}\mathrm{d}\textbf{w}$
(9)
$\displaystyle-\beta_{t}\nabla_{\hat{\textbf{x}}_{0}}\|\textbf{BS}\exp({-\textbf{QA}\hat{\textbf{x}}_{0}})-\textbf{y}\|_{\textbf{K}^{-1}}^{2}\nabla_{\textbf{x}_{t}}{\hat{\textbf{x}}_{0}}\mathrm{d}t$

The SDPS pseudo-code is shown below and illustrated in Fig. 1.

Algorithm 1 Spectral Diffusion Posterior Sampling (SPDS)

1:$T$: diffusion steps

2:$\eta_{t}$: step size

3:y: spectral CT measurement

4:$\textbf{x}_{T}\sim\mathcal{N}~{}(0,\boldsymbol{I})$

5:for $t=T$ to $1$ do:

6: $\textbf{z}\sim\mathcal{N}~{}(0,\boldsymbol{I})$

7:
$\hat{\textbf{x}}_{0}=\frac{1}{\sqrt{\bar{\alpha}_{t}}}(\textbf{x}_{t}-\sqrt{1-\bar{\alpha}_{t}}\boldsymbol{\epsilon}_{\theta}(\textbf{x}_{t},t))$

8:
$\textbf{x}_{t-1}^{\prime}=\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}\textbf{x}_{t}+\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}}\hat{\textbf{x}}_{0}+\sigma_{t}\textbf{z}$

9:
$\textbf{x}_{t-1}=\textbf{x}_{t-1}^{\prime}-\eta_{t}\nabla_{\textbf{x}_{t}}\left\|\textbf{BS}\exp({-\textbf{QA}\hat{\textbf{x}}_{0}})-\textbf{y}\right\|_{\textbf{K}^{-1}}^{2}\nabla_{\textbf{x}_{t}}{\hat{\textbf{x}}_{0}}$

10:end for

Figure 1: SDPS Workflow: Basis material images are concatenated to form
$\textbf{x}_{0}$. The forward process progressively perturbs images with
noise, and a noise prediction network is trained. The reverse process uses
DDPM sampling to progressively denoise the image and uses a likelihood-based
update to enforce data consistency.

###  II-B Strategies for Fast and Stable SPDS

####  II-B1 Instability of SDPS

The SDPS framework has the advantage of combining a learned prior and a
physics model. However, we have observed that the stability of SDPS is highly
dependent on the step size, $\eta_{t}$, and that even with careful tuning,
decompositions present large variations over posterior samples and
hallucinations. Careful design of a step size scheduler is required to balance
the prior information and measurements, particularly for the earlier steps in
the reverse process. In early steps, $\textbf{x}_{t}$ is far from the
solution. While large step sizes may be favorable to enhance data consistency,
this risks disrupting diffusion updates - making it challenging to pick
optimal step sizes.

####  II-B2 Jumpstarted Sampling

For many imaging problems, an initial estimate is readily available through
fast computation. For example, in spectral CT, projection-domain or image-
domain material decomposition [17, 18] can be used as the first-pass estimate,
which we denote as $\hat{\textbf{x}}_{0}^{f}$. Forward diffusion can be
directly performed on $\hat{\textbf{x}}_{0}^{f}$ as
$\hat{\textbf{x}}_{t}^{f}=\sqrt{\bar{\alpha}}_{t}\hat{\textbf{x}}_{0}^{f}+(\sqrt{1-\bar{\alpha}_{t}})\boldsymbol{\epsilon}$.
The difference between the distribution of $\hat{\textbf{x}}_{t}$ and
$\hat{\textbf{x}}_{t}^{f}$ can be quantified by the KL divergence:

$D_{KL}(p(\textbf{x}_{t}|\textbf{x}_{0})\parallel
p(\hat{\textbf{x}}_{t}^{f}|\hat{\textbf{x}}_{0}^{f}))=\frac{\bar{\alpha}_{t}}{2(1-\bar{\alpha}_{t})}\|\hat{\textbf{x}}_{0}^{f}-\textbf{x}_{0}\|_{2}^{2}.$
(10)

In forward diffusion, ${\bar{\alpha}_{t}}/{(1-\bar{\alpha}_{t})}$
progressively diminishes, ultimately converging towards zero. Consequently,
Eq.(10) implies that we can expect for sufficiently large $t$, the difference
between the distribution of $\textbf{x}_{t}$ and $\hat{\textbf{x}}_{t}^{f}$ is
small enough to be ignored, then the reverse sampling can start from
$\hat{\textbf{x}}_{0}^{f}$ instead of pure noise. In effect, we can skip early
time steps by using an approximate solution to which an appropriate amount of
noise has been added - yielding a more stabilized decomposition by using an
initialization closer the solution as well as a faster solution having
“jumpstarted” over many early time steps.

####  II-B3 Gradient Approximation

SDPS computes the data fidelity gradient term via chain rule:
$\nabla_{\textbf{x}_{t}}\mathrm{log}p_{t}({\textbf{y}|\hat{\textbf{x}}_{0}})=\nabla_{\hat{\textbf{x}}_{0}}\mathrm{log}p_{t}({\textbf{y}|\hat{\textbf{x}}_{0}})\nabla_{\textbf{x}_{t}}{\hat{\textbf{x}}_{0}}$.
The $\nabla_{\textbf{x}_{t}}{\hat{\textbf{x}}_{0}}$ term may be expanded as:

$\nabla_{\textbf{x}_{t}}{\hat{\textbf{x}}_{0}}=\frac{1}{\sqrt{\bar{\alpha}_{t}}}(1-\sqrt{1-\alpha_{t}}\nabla_{\textbf{x}_{t}}\boldsymbol{\epsilon}_{\theta}(\textbf{x}_{t},t))\vspace{-.1in}$
(11)

The computation of the Jacobian
$\nabla_{\textbf{x}_{t}}\boldsymbol{\epsilon}_{\theta}(\textbf{x}_{t},t)$ is
both time- and memory-intensive, particularly for high dimensional images and
deep neural networks. We have experimentally observed that the Jacobian is
well-approximated by a diagonal matrix for each time step, and the diagonal
elements are of limited range. Therefore, we approximate the Jacobian by a
constant (which could be further absorbed into the step size). Incorporating
both the jumpstarted sampling and the gradient approximation, we propose the
jumpstarted spectral DPS (JSDPS), which is outlined in Algorithm 2 below.

Algorithm 2 Jumpstarted Spectral DPS (JSDPS)

1:$T,T^{\prime}$: training steps, sampling steps, $T^{\prime}\ll T$

2:$\eta_{t}$: step size

3:y: spectral CT measurement

4:$\hat{\textbf{x}}_{0}^{f}$: image-domain decomposition

5:$\boldsymbol{\epsilon}\sim\mathcal{N}~{}(0,I)$

6:$\textbf{x}_{T^{\prime}}=\hat{\textbf{x}}_{0}^{T^{\prime}}=\sqrt{\bar{\alpha}_{T^{\prime}}}\textbf{x}_{0}+\sqrt{1-\bar{\alpha}_{T^{\prime}}}\boldsymbol{\epsilon}$

7:for $t=T^{\prime}$ to $1$ do:

8: $\textbf{z}\sim\mathcal{N}~{}(0,\boldsymbol{I})$

9:
$\hat{\textbf{x}}_{0}=\frac{1}{\sqrt{\bar{\alpha}_{t}}}(\textbf{x}_{t}-\sqrt{1-\bar{\alpha}_{t}}\boldsymbol{\epsilon}_{\theta}(\textbf{x}_{t},t))$

10:
$\textbf{x}_{t-1}^{\prime}=\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}\textbf{x}_{t}+\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}}\hat{\textbf{x}}_{0}+\sigma_{t}\textbf{z}$

11:
$\textbf{x}_{t-1}=\textbf{x}_{t-1}^{\prime}-\eta_{t}\nabla_{\hat{\textbf{x}}_{0}}\left\|\textbf{BS}\exp({-\textbf{QA}\hat{\textbf{x}}_{0}})-\textbf{y}\right\|_{\textbf{K}^{-1}}^{2}$

12:end for

###  II-C Unconditional Multi-Material Generation

####  II-C1 Dataset Generation

This work focuses on two material decomposition with water and calcium as
bases. To build spectral CT and material density dataset, we collected 32000
clinical chest CT slices from the CT Lymph Nodes dataset[19]. Those images are
pre-processed to convert Hounsfield Units to attenuation coefficients and
remove patient beds. An upper bound value (2000HU) was also applied to remove
large attenuation values from metal. Two soft threshold functions were applied
to approximate water and calcium densities (unit: $g/ml$):

$\rho_{w}(\mu)=\begin{cases}k_{w}\mu,&\text{if }\mu\leq\mu_{w}\\\
k_{w}\mu_{w}-k_{wc}(\mu-\mu_{w}),&\text{if }\mu_{w}<\mu<\mu_{c}\\\
0,&otherwise\end{cases}$ (12a)
$\rho_{c}(\mu)=\begin{cases}k_{c}(\mu-\mu_{c})+k_{cw}(\mu_{c}-\mu_{w}),&\text{if
}\mu\geq\mu_{c}\\\ k_{cw}(\mu-\mu_{w}),&\text{if }\mu_{w}<\mu<\mu_{c}\\\
0,&otherwise\end{cases}$ (12b)

In this work, $k_{w}$, $k_{wc}$, $k_{cw}$, $k_{c}$ were empirically set to
$5.18$, $-8.77$, $5.69$, $2.12~{}g/cm^{2}$, while $\mu_{w}$ and $\mu_{c}$ were
$0.22~{}cm^{-1}$ and $0.35~{}cm^{-1}$. The water and calcium densities were
used as ground truth to simulate spectral measurements.

####  II-C2 Network Training

DDPM[15] is employed as the SGM framework, with the Residual Unet as the
backbone network. Paired water and calcium images are concatenated to form the
2 channel image $\textbf{x}_{0}$. The discretized diffusion process uses
$T=1000$ time steps with a linear variance scheduler from $\beta_{1}=1e^{-4}$
to $\beta_{1000}=0.02$. Training and validation datasets contain 25000 slices
and 5000 slices, respectively. The loss function (7) is minimized by the Adam
optimizer with a batch size of $16$ and a learning rate of $10^{-4}$. Training
terminated after 200 epochs.

###  II-D Conditional Posterior Sampling for Different Spectral Systems

The unconditional DDPM training captures the material bases distribution
without specification of the spectral CT device. The trained model can be used
in material decomposition for arbitrary spectral CT systems. Here we
investigate the SDPS/JSDPS performance on a simulated dual-layer CT [13] and
dual-kVp CT system [12]. The dual-layer CT uses a detector with $300\mu m$ CsI
and $600\mu m$ CsI scintillators in the top and bottom layers, respectively
(with a 5mm gap between layers). The x-ray tube operates at 120 kVp. The dual-
kVp system uses single-layer detector with $600\mu m$ CsI. The tube voltage
alternates between 80 kVp and 120 kVp every other view. For both systems, 800
projections were simulated with Poisson noise equivalent to an exposure of
0.05 mAs/view. Voxel size and detector pixel size were set to $0.8~{}mm$ and
$1.0~{}mm$, respectively.

###  II-E Evaluation

We first evaluated the performance of JSDPS with and without gradient
approximation. A standard image-domain decomposition from filtered
backprojection reconstructions, i.e., the initialization of the JSDPS
algorithm, was also shown for comparison. Second, we compared the performance
of several material decomposition algorithms including model-based material
decomposition (MBMD) [3], SDPS, and JSDPS. MBMD was formulated as the
following optimization problem:

$\hat{\textbf{x}}=\text{argmin}\|\textbf{BS}\exp({-\textbf{QAx}})-\textbf{y}\|_{\textbf{K}^{-1}}^{2}+\lambda_{w}\textbf{R}(\textbf{x}_{w})+\lambda_{c}\textbf{R}(\textbf{x}_{c})$
(13)

where R is a quadratic gradient roughness penalty. Regularization strengths,
$\lambda_{w},\lambda_{c}$, for water and calcium were set to $10^{-4}$ and
$4\times 10^{-4}$, respectively. We used 10000 iterations of separable
paraboloidal surrogates updates, and final change in loss is $\sim 0.01\%$
from the previous iteration. SDPS was implemented according to Algorithm 1,
and the step size scheduler follows [10]:
$\eta_{t}=\eta/\|\textbf{BS}\exp({-\textbf{QAx}})-\textbf{y}\|_{\textbf{K}^{-1}}$.
We swept $\eta$ from $0.1$ to $10$, and the optimal $\eta$ was determined as
the one with minimal MSE from ground truth. JSDPS used a constant step
scheduler with $T^{\prime}=150,\eta_{t}=0.02$. Performance was evaluated for
computation time (calculated for eight parallel decompositions), and image
quality - quantified by Structural Similarity Index Measure(SSIM) and Peak
Signal-to-Noise Ratio (PSNR).

Additionally, to compare variability and biases in the posterior samples
provided by SDPS and JSDPS, we compute sample bias and standard deviation over
an ensemble of 16 outputs:

$\begin{matrix}\text{Bias}=\|\mathbb{E}\\{\hat{\textbf{x}}\\}-\textbf{x}\|_{2}&\text{Std}=\|\mathbb{E}\\{(\hat{\textbf{x}}-\mathbb{E}\\{\hat{\textbf{x}}\\})^{2}\\}\|_{2}\end{matrix}$
(14)

##  III Results

Figure 2: Decomposed water (top, W/L: $1.2/0.6~{}g/ml$) and Calcium (bottom,
W/L: $0.05/0.1~{}g/ml$) images. Left to right: Image-domain decomposition,
JSDPS without gradient approximation, JSDPS with gradient approximation,
ground truth. Computational time for eight outputs is provided in the bottom
left corner. Figure 3: Decomposed water (top, W/L: $1.2/0.6~{}g/ml$) and
Calcium (bottom, W/L: $0.05/0.1~{}g/ml$) for (left to right): MBMD, SDPS,
JSDPS, ground truth.

###  III-A Illustration of JSDPS variants

Fig.2 displays JSDPS results as compared with standard image-domain
decomposition (initialization) and the ground truth. Image-domain
decomposition shows significant image noise and streaking. JSDPS without
gradient approximation effectively captures most of the high-contrast
structures, including smaller pulmonary vessels. However, as seen in the
zoomed-in ROIs, low-contrast soft tissue hallucinations are present. In the
contrast enhanced areas of the calcium image, such as the aorta, deformations
are also observed. By using the gradient approximation, the reverse sampling
achieves $20.2\%$ and $29.3\%$ computation time decrease for the dual-layer
and dual-kVp system, respectively. Gradient approximation also effectively
mitigates the hallucinations and enhances the consistency with ground truth
for vascular structures.

Figure 4: Bias and standard deviation map for SDPS and JSDPS. Here we only
show the results of dual-kVp system, which is similar to that of the dual-
layer system.

###  III-B Comparison of SDPS, JSDPS, and MBMD

####  III-B1 Image Quality

Figure 3 compares the performance of SDPS, JSDPS, and MBMD in water and
calcium bases for the dual layer and dual kV systems. Four zoomed in areas are
shown in the (1) lung, (2) heart (iodine-filled aorta in the calcium basis),
(3) low-contrast soft tissue, and (4) spine. Comparing across all three
algorithms, MBMD produces images with the highest noise (evident in ROIs 2 and
3) at the lowest spatial resolution (evident in ROI 1). The SDPS results have
low noise and high spatial resolution, but are prone to hallucinations
observed as erroneous anatomy in all four ROIs. While SDPS produces results
that are more closely aligned with the prior distribution, the instability of
the algorithm is problematic for the low fidelity (low dose) data simulated in
this work. The JSDPS algorithm outperforms both MBMD and SDPS from both visual
observations of anatomical structures faithful to the ground truth and
according to SSIM and PSNR. Good decomposition results were obtained for both
the dual layer and dual kV system, demonstrating generalizability to different
spectral systems.

Figure 4 shows the bias and standard deviation maps computed from an ensemble
of 16 individual output samples. SDPS exhibits large variability and bias
around edge features, suggesting substantial image uncertainty. In contrast,
JSDPS displays markedly reduced bias and standard deviation, demonstrating the
capability to effectively stabilize the sampling process.

####  III-B2 Computational Cost

Table I summarizes the computational cost of different material decomposition
algorithms for eight output estimates. MBMD is the slowest due to the large
number of iterations (10000) required for convergence. SDPS substantially
reduces the number of iterations (1000), thereby shortening the computation
time. Leveraging the jumpstarted sampling strategy and the gradient
approximation, JSDPS requires only 150 iterations, which results in an
additional $86.83\%$ and $88.07\%$ reduction in computation time for dual-
layer and dual kVp CT, respectively. Furthermore, a $23.16\%$ and $42.54\%$
memory saving was also achieved compared to SDPS for the two spectral systems
due to the absence of Jacobian computation.

TABLE I: Computational cost of material decomposition algorithms | Dual layer | Dual kVp
---|---|---
| Time(s) | Memory(GB) | Time(s) | Memory(GB)
MBMD | 6406.3 | 11.3 | 2802.5 | 4.8
SDPS | 728.2 | 17.7 | 478.6 | 13.4
JSDPS | 95.9 | 13.6 | 57.1 | 7.7

##  IV Conclusion and Discussion

Spectral CT material decomposition is a challenging and ill-conditioned
problem, often exhibiting excessive noise and slow convergence. The novel SDPS
framework is designed to tackle these issues - combining a sophisticated
learned SGM prior and a physical model, to achieve low image noise while
enhancing image quality. We additionally adopted jumpstarted sampling and
gradient approximation strategies which were found to be effective in
stabilizing performance as well as reducing computation time and memory
consumption. In evaluations on both dual-layer and dual-kVp systems, JSDPS
achieves the highest accuracy, the lowest uncertainty, and the lowest
computational costs compared to SDPS and MBMD. This work demonstrated that
JSDPS is a promising approach for spectral CT decomposition.

## Acknowledgments

This work is supported, in part, by NIH grant R01EB030494.

## References

  * [1] Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstruction for spectral ct,” _IEEE transactions on medical imaging_ , vol. 33, no. 8, pp. 1614–1626, 2014.
  * [2] X. Jiang, C. Fang, P. Hu, H. Cui, L. Zhu, and Y. Yang, “Fast and effective single-scan dual-energy cone-beam ct reconstruction and decomposition denoising based on dual-energy vectorization,” _Medical physics_ , vol. 48, no. 9, pp. 4843–4856, 2021.
  * [3] S. Tilley, W. Zbijewski, and J. W. Stayman, “Model-based material decomposition with a penalized nonlinear least-squares ct reconstruction algorithm,” _Physics in Medicine & Biology_, vol. 64, no. 3, p. 035005, 2019.
  * [4] H. Gong, S. Tao, K. Rajendran, W. Zhou, C. H. McCollough, and S. Leng, “Deep-learning-based direct inversion for material decomposition,” _Medical physics_ , vol. 47, no. 12, pp. 6294–6309, 2020.
  * [5] W. Wu, D. Hu, C. Niu, L. V. Broeke, A. P. Butler, P. Cao, J. Atlas, A. Chernoglazov, V. Vardhanabhuti, and G. Wang, “Deep learning based spectral ct imaging,” _Neural Networks_ , vol. 144, pp. 342–358, 2021.
  * [6] J. F. Abascal, N. Ducros, V. Pronina, S. Rit, P.-A. Rodesch, T. Broussaud, S. Bussod, P. C. Douek, A. Hauptmann, S. Arridge _et al._ , “Material decomposition in spectral ct using deep learning: a sim2real transfer approach,” _IEEE Access_ , vol. 9, pp. 25 632–25 647, 2021.
  * [7] J. Zhu, T. Su, X. Zhang, J. Yang, D. Mi, Y. Zhang, X. Gao, H. Zheng, D. Liang, and Y. Ge, “Feasibility study of three-material decomposition in dual-energy cone-beam ct imaging with deep learning,” _Physics in Medicine & Biology_, vol. 67, no. 14, p. 145012, 2022.
  * [8] W. Fang, D. Wu, K. Kim, M. K. Kalra, R. Singh, L. Li, and Q. Li, “Iterative material decomposition for spectral ct using self-supervised noise2noise prior,” _Physics in Medicine & Biology_, vol. 66, no. 15, p. 155013, 2021.
  * [9] A. Eguizabal, O. Öktem, and M. U. Persson, “Deep learning for material decomposition in photon-counting ct,” _arXiv preprint arXiv:2208.03360_ , 2022.
  * [10] H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, and J. C. Ye, “Diffusion posterior sampling for general noisy inverse problems,” _arXiv preprint arXiv:2209.14687_ , 2022.
  * [11] S. Li, M. Tivnan, Y. Shen, and J. W. Stayman, “Diffusion posterior sampling for nonlinear ct reconstruction,” _arXiv preprint arXiv:2312.01464_ , 2023.
  * [12] R. Cassetta, M. Lehmann, M. Haytmyradov, R. Patel, A. Wang, L. Cortesi, D. Morf, D. Seghers, M. Surucu, H. Mostafavi _et al._ , “Fast-switching dual energy cone beam computed tomography using the on-board imager of a commercial linear accelerator,” _Physics in Medicine & Biology_, vol. 65, no. 1, p. 015013, 2020.
  * [13] W. Wang, Y. Ma, M. Tivnan, J. Li, G. J. Gang, W. Zbijewski, M. Lu, J. Zhang, J. Star-Lack, R. E. Colbeth _et al._ , “High-resolution model-based material decomposition in dual-layer flat-panel cbct,” _Medical physics_ , vol. 48, no. 10, pp. 6375–6387, 2021.
  * [14] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” _arXiv preprint arXiv:2011.13456_ , 2020.
  * [15] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” _Advances in neural information processing systems_ , vol. 33, pp. 6840–6851, 2020.
  * [16] J. Song, C. Meng, and S. Ermon, “Denoising diffusion implicit models,” _arXiv preprint arXiv:2010.02502_ , 2020.
  * [17] P. Stenner, T. Berkus, and M. Kachelriess, “Empirical dual energy calibration (edec) for cone-beam computed tomography,” _Medical physics_ , vol. 34, no. 9, pp. 3630–3641, 2007.
  * [18] P. R. Mendonça, P. Lamb, and D. V. Sahani, “A flexible method for multi-material decomposition of dual-energy ct images,” _IEEE transactions on medical imaging_ , vol. 33, no. 1, pp. 99–116, 2013.
  * [19] H. R. Roth, L. Lu, A. Seff, K. M. Cherry, J. Hoffman, S. Wang, J. Liu, E. Turkbey, and R. M. Summers, “A new 2.5 d representation for lymph node detection using random sets of deep convolutional neural network observations,” in _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014: 17th International Conference, Boston, MA, USA, September 14-18, 2014, Proceedings, Part I 17_. Springer, 2014, pp. 520–527.


# Controlization Schemes Based on Orthogonal Arrays

Anirban Chowdhury Ewout van den Berg Pawel Wocjan
IBM Quantum, IBM Thomas J. Watson Research Center, Yorktown Heights, NY, USA

###### Abstract

Realizing controlled operations is fundamental to the design and execution of
quantum algorithms. In quantum simulation and learning of quantum many-body
systems, an important subroutine consists of implementing a controlled
Hamiltonian time-evolution. Given only black-box access to the uncontrolled
evolution $e^{-iHt}$, _controlizing_ it, i.e., implementing
$\mathrm{ctrl}(e^{-iHt})=|0\rangle\langle 0|\otimes I+|1\rangle\langle
1|\otimes e^{-iHt}$ is non-trivial. A recent paper introduced a controlization
scheme in the context of designing algorithms for transforming unknown
Hamiltonian dynamics [32]. The main idea behind the scheme is to intersperse
the uncontrolled evolution with suitable operations such that the overall
dynamics approximates the desired controlled evolution. The scheme in Ref.
[32], although efficient, uses operations randomly sampled from an
exponentially large set. In the present work, we show that more efficient
controlization schemes can be constructed with the help of orthogonal arrays
for unknown 2-local Hamiltonians. This construction can also be generalized to
$k$-local Hamiltonians. Moreover, our controlization schemes based on
orthogonal arrays can take advantage of the interaction graph’s structure and
be made more efficient.

##  1 Introduction

The task of efficiently realizing controlled operations plays a fundamental
role in the design and execution of quantum algorithms. Algorithmic primitives
such as quantum phase estimation [31], quantum amplitude estimation [3],
linear-combination-of-unitaries (LCU) [13, 26], and quantum singular value
transformation (QSVT) [28, 18] require implementing a _controlled_ unitary
$U$, an operation of the form $\mathrm{ctrl}(U)=|0\rangle\langle 0|\otimes
I+|1\rangle\langle 1|\otimes U$. In applications related to quantum
simulation, the unitary $U$ is often the time-evolution under a Hamiltonian
$H$ for time $t$, i.e., $U=e^{-iHt}$. Quantum phase estimation, for example,
is used with controlled time-evolution to determine spectral properties of
Hamiltonians. Block-encoding techniques are used to manipulate controlled
time-evolution into implementing functions of $H$ – this idea is foundational
to quantum algorithms for preparing both ground and Gibbs states of quantum
many-body systems [34, 35, 12, 15, 39, 14], solving linear systems of
equations [19, 11], and many more [28, 18].

In most of these examples, the Hamiltonian $H$ is known beforehand. Thus there
is an explicit quantum circuit that realizes the time-evolution and turning
that into a controlled operation is straightforward. There exist problems, for
instance, in quantum meteorology [33] and quantum learning [43, 29, 1], where
the Hamiltonian may not be known a priori. However, we may still wish to use
the quantum algorithmic primitives mentioned above to design efficient
algorithms. This raises the question: “Can we _controlize_ an unknown
Hamiltonian evolution? That is, given the ability to perform $\exp{(-iHt)}$ as
a black-box, can we implement the operation
$\mathrm{ctrl}\left(\exp{(-iHt)}\right)$?”

Controlization of unknown quantum dynamics has been studied in a number of
different settings [21, 16, 32]. Janzing [21] showed how a controlization
protocol can be used to perform phase estimation of unknown quantum
Hamiltonians. Closely related are also the problems of reversing and fast-
forwarding Hamiltonian dynamics [30, 42]. While our focus here is on
Hamiltonian dynamics, it should be noted that controlizing arbitrary unitary
operations has also been investigated extensively [2, 10, 36, 16]. More
recently, controlization has been used as a primitive in algorithms for
transforming unknown Hamiltonian dynamics to implement functions of
Hamiltonians [32].

The controlization protocol for unknown Hamiltonian dynamics in Ref. [32]
proceeds as follows. It invokes the black-box to implement the Hamiltonian
dynamics $\exp(-iH\tau)\in\mathcal{U}(\mathcal{H})$ for different times
$\tau>0$ and intersperses them with different control operations from a finite
subset of $\mathcal{U}(\mathbb{C}^{2}\otimes\mathcal{H})$ such that the
resulting evolution approximates the unitary

$\displaystyle\mathrm{ctrl}\big{(}\exp(-iHt)\big{)}$
$\displaystyle=|0\rangle\langle 0|\otimes I+|1\rangle\langle
1|\otimes\exp(-iHt)\in\mathcal{U}(\mathbb{C}^{2}\otimes\mathcal{H})$ (1)

for some desired value of $t$. More precisely, the unknown Hamiltonian $H$
acts on $n$ qubits, that is, $\mathcal{H}=(\mathbb{C}^{2})^{\otimes n}$ and
the control operations are controlled $n$-fold tensor products of the Pauli
matrices $\\{I,X,Y,Z\\}$. In particular, the operations are randomly sampled
from an exponentially large set.

The purpose of this work is to present efficient controlization schemes for
the situation where the quantum system consists of $n$ qudits, that is,
$\mathcal{H}=(\mathbb{C}^{d})^{\otimes n}$, and the unknown system Hamiltonian
$H\in\mathcal{L}(\mathcal{H})$ is assumed to be $k$-local with $k$ being a
constant. Our construction of controlization schemes is based on the
combinatorial concept of orthogonal arrays [20] and improves upon some aspects
of the controlization schemes presented in [32]. It significantly reduces the
size of the set of operations used to intersperse the black-box Hamiltonian
dynamics. The main idea underlying our construction is that there is a close
connection between controlization and decoupling schemes and that the latter
can be constructed from orthogonal arrays.

##  2 Controlization of known Hamiltonian dynamics

Before considering the case of unknown Hamiltonian dynamics, let us illustrate
how a decoupling scheme for a Hamiltonian $H$ for which the terms are
explicitly known can be naturally extended to a controlization scheme.

For simplicity assume initially that $H=P$, where $P$ is a Pauli operator.
Choosing a Pauli $Q$ that anti-commutes with $P$ enables us to decouple, i.e.,
to switch off the time evolution, since

$P+QPQ^{{\dagger}}=P-P=0$

and thus

$(e^{-iPt/2})(Qe^{-iPt/2}Q^{{\dagger}})=(e^{-iPt/2})(e^{-iQPQ^{{\dagger}}t/2})=(e^{-iPt/2})(e^{+iPt/2})=I.$

Omitting $Q$ and $Q^{{\dagger}}$ in the product above yields

$(e^{-iPt/2})(e^{-iPt/2})=e^{-iPt}$

This allows us to implement the conditional time evolution with $P$ as shown
in the circuit below:

Since we can always choose $Q$ to be a weight-one Pauli, the advantage of this
approach is that we replace a conditional multi-qubit Pauli exponentiation
gate by one that merely requires two conditional single-qubit Pauli
operations.

Time evolution for a known Hamiltonian $H=\sum_{i}\alpha_{i}H_{i}$ consisting
of a weighted sum of Pauli operators can be implemented by means of
Trotterization or schemes such as QDRIFT [8, 24]. The latter is a
probabilistic simulation technique that often outperforms deterministic
simulation techniques.

A conditional time evolution is then obtained by splitting the time evolution
of each individual term in two pieces and conjugate one of them by an operator
that negates the term. Various circuit simplifications can be made. For
instance, if the Pauli terms $A$ and $B$ for successive exponentiation terms
$e^{iAt}$ and $e^{iBt}$, have overlapping support we can always find a unit-
weight Pauli $Q$ on that support that anti-commutes with both terms. This
allows us to write:

The controlled $Q$ and $Q^{{\dagger}}$ gates originally present between
$e^{-iAt/2}$ and $e^{-iBt/2}$ cancel and are therefore omitted. Grouping the
Pauli terms into commuting sets and appropriately choosing the ordering of the
terms allows additional simplification along the same lines.

##  3 Controlization of unknown Hamiltonian dynamics based on decoupling
schemes

The goal of dynamical decoupling is to effectively stop the system from
evolving by interspersing the time evolution according to the system
Hamiltonian $H$ by suitable control operations. We refer the reader to [45,
46, 44] for an introduction to a combinatorial approach for constructing
Hamiltonian simulation schemes based on average Hamiltonian theory. Decoupling
and time-reversal are two special cases of these schemes. We make use of the
combinatorial decoupling schemes presented in these works to construct our
controlization schemes.

In the current section we do not assume any particular structure of the
quantum system $\mathcal{H}$ and its Hamiltonian
$H\in\mathcal{L}(\mathcal{H})$. We only assume the Hamiltonian $H$ to be
traceless throughout. We also assume that the control operations in the
control group $\mathcal{G}$ – a finite subgroup of $\mathcal{U}(\mathcal{H})$
– can be implemented instantaneously (this is often referred to as bang-bang
control).111We refer the reader to the paper [7, 5] and the references therein
for Hamiltonian simulation based on bounded-strength controls. In the next
section we consider quantum systems consisting of qudits and local
Hamiltonians describing the couplings among the qudits.

To illustrate the idea of Hamiltonian simulation based on average Hamiltonian
theory, let us consider a simple example that uses three control operations
$V_{1}$, $V_{2}$, $V_{3}\in\mathcal{G}$ and invokes the black-box to implement
the dynamics for two times $\tau_{1}$ and $\tau_{2}$. We can rewrite the
resulting unitary as follows:

$\displaystyle V_{3}\,e^{-iH\tau_{2}}\,V_{2}\,e^{-iH\tau_{1}}\,V_{1}$ (2)
$\displaystyle=(V_{3}V_{2}V_{1})\,(V_{2}V_{1})^{\dagger}\,e^{-iH\tau_{2}}\,(V_{2}V_{1})\,V_{1}^{\dagger}\,e^{-iH\tau_{1}}\,V_{1}$
(3)
$\displaystyle=e^{-iU_{2}^{\dagger}HU_{2}\tau_{2}}\,e^{-iU_{1}^{\dagger}HU_{1}\tau_{1}}$
(4)

The operators $U_{1}$ and $U_{2}$ are given by

$\displaystyle U_{1}$ $\displaystyle=V_{1}$ (5) $\displaystyle U_{2}$
$\displaystyle=V_{2}V_{1}$ (6)

assuming that the final control operation $V_{3}$ is chosen such that
$V_{3}V_{2}V_{1}=I$ (the latter means that the control sequence is cyclic).

Obviously, the above method of rewriting the resulting time evolution
generalizes to an arbitrary number of control operations. We only consider
cyclic control schemes. We have:

$\displaystyle U_{1}$ $\displaystyle=V_{1}$ (7) $\displaystyle U_{2}$
$\displaystyle=V_{2}V_{1}$ (8) $\displaystyle\,\,\,\vdots$ $\displaystyle
U_{N}$ $\displaystyle=V_{N}\cdots V_{2}V_{1}$ (9)

We always chose $V_{N+1}=(V_{N}\cdots V_{2}V_{1})^{\dagger}$ to make the
scheme cyclic.

Using the simple Trotter approximation, we obtain

$\displaystyle\prod_{j=1}^{N}e^{-iU_{j}^{\dagger}HU_{j}\tau_{j}}$
$\displaystyle=e^{-iU_{N}^{\dagger}HU_{N}\tau_{N}}\,\cdots\,e^{-iU_{2}^{\dagger}HU_{2}\tau_{2}}\,e^{-iU_{1}^{\dagger}HU_{1}\tau_{1}}$
(10) $\displaystyle\approx e^{-i\tilde{H}},$ (11)

where the average Hamiltonian $\tilde{H}$ is given by the sum

$\displaystyle\tilde{H}$
$\displaystyle=\tau_{1}U_{1}HU_{1}^{\dagger}+\tau_{2}U_{2}HU_{2}^{\dagger}+\ldots+\tau_{N}U_{N}HU_{N}^{\dagger}.$
(12)

Thus, after the control cycle the resulting time evolution is approximately as
if the system had evolved under the average Hamiltonian $\tilde{H}$. This is
referred to as average Hamiltonian theory.

A better approximation can be achieved according to the formula:

$\displaystyle\prod_{j=N}^{1}e^{-iU_{j}^{\dagger}HU_{j}\tau_{j}/2}\prod_{j=1}^{N}e^{-iU_{j}^{\dagger}HU_{j}\tau_{j}/2}$
(13)

Take note that the two products above multiply the exponentials
$e^{-iU_{j}^{\dagger}Hu_{j}\tau_{j}/2}$ in opposite orders and use times
$\tau_{j}/2$. This formula is a special case of the Lie-Trotter-Suzuki
formulas [40, 41]. Observe that these higher-order formulas do not increase
the overall time of the control cycle, but only increase the number of control
operations.

The average Hamiltonian in eq. (12) serves as the starting point for these
higher-order formulas. This is why it is important to express the desired
target Hamiltonian $\tilde{H}$ as a weighted sum containing a small number $N$
of conjugates $U_{1}^{\dagger}HU_{1},\ldots,U_{N}^{\dagger}HU_{N}$ and having
a small overall time $\sum_{j=1}^{N}\tau_{j}$.

######  Definition 1 (Simulation, decoupling, and time reversal schemes).

Let $H,\tilde{H}\in\mathcal{L}(\mathcal{H})$ be two Hamiltonians. We say that
a scheme $\mathcal{S}=(\tau_{1},U_{1};\ldots;\tau_{N},U_{N})$ is a simulation
scheme for the system Hamiltonian $H$ and the target Hamiltonian $\tilde{H}$
if

$\displaystyle\mathcal{S}(H)$
$\displaystyle=\sum_{j=1}^{N}\tau_{j}U_{j}HU_{j}^{\dagger}=\tilde{H}.$ (14)

There are two important special cases. For $\tilde{H}=0_{\mathcal{H}}$, where
$0_{\mathcal{H}}\in\mathcal{L}(\mathcal{H})$ denotes the zero operator, we
call the scheme a _decoupling_ scheme222For decoupling schemes, we always use
$\sum_{j=1}^{N}\tau_{j}=1$. and use $\mathcal{D}$ to denote it. For
$\tilde{H}=-H$, we call the scheme a _time reversal_ scheme and use
$\mathcal{R}$ to denote it.

######  Remark 1 (Decoupling $\Rightarrow$ time reversal).

Before discussing in detail how decoupling schemes can be extended to
controlization schemes, let us briefly describe the connection to time
reversal schemes. Let
$\mathcal{D}=(\tau_{1},U_{1};\tau_{2},U_{2};\ldots;\tau,U_{N})$ be a
decoupling scheme. We may assume w.l.o.g. that $U_{1}=I_{\mathcal{H}}$. We
also assume that $\sum_{j}t_{j}=1$ as explained in the footnote below. Then,
we have

$\displaystyle\tau_{1}H+\tau_{2}U_{2}HU_{2}^{\dagger}+\ldots\tau_{N}U_{N}HU_{N}^{\dagger}$
$\displaystyle=0,$ (15)

which is equivalent to

$\displaystyle\frac{\tau_{2}}{\tau_{1}}\,U_{2}HU_{2}^{\dagger}+\ldots+\frac{\tau_{2}}{\tau_{1}}\
U_{N}HU_{N}^{\dagger}$ $\displaystyle=-H.$ (16)

Thus, the decoupling scheme
$\mathcal{D}=(\tau_{1},U_{1}=I;\tau_{2},U_{2};\ldots;\tau_{N},U_{N})$ gives
rise to the time reversal scheme
$\mathcal{R}=(\tau_{2}/\tau_{1},U_{2};\ldots;\tau_{N}/\tau_{N},U_{N})$.
Observe that the time needed to realize time reversal is
$(1-\tau_{1})/\tau_{1}$ which will be greater than
$1=\tau_{1}+\ldots+\tau_{N}$ in contrast to the situation for decoupling. This
slow-down for time-reversal is discussed in [23] and more generally for
Hamiltonian simulation in [45].

Before we state the theorem connecting controlization and decoupling, we need
to introduce the following definition. Let $M\in\mathcal{L}(\mathcal{H})$ be
an arbitrary operator. Then, $\Lambda(M)$ is the controlled operator defined
to be

$\displaystyle\Lambda(M)$ $\displaystyle=|0\rangle\langle 0|\otimes
M+|1\rangle\langle 1|\otimes
I_{\mathcal{H}}\in\mathcal{L}(\mathbb{C}^{2}\otimes\mathcal{H}),$ (17)

where $I_{\mathcal{H}}\in\mathcal{L}(\mathcal{H})$ denotes the identity
operator.

######  Theorem 1 (Decoupling $\Rightarrow$ controlization).

Let $\mathcal{D}=(\tau_{1},U_{1};\ldots;\tau_{N},U_{N})$ be a decoupling
scheme for some Hamiltonian $H\in\mathcal{L}(\mathcal{H})$. Then, the scheme
$\mathcal{C}=\Lambda(\mathcal{D})=(\tau_{1},\Lambda(U_{1});\ldots;\tau_{N},\Lambda(U_{N}))$
is a simulation scheme for the pair

$\displaystyle I_{2}\otimes H\mbox{ and }|1\rangle\langle 1|\otimes H$ (18)

of Hamiltonians in $\mathcal{L}(\mathbb{C}^{2}\otimes\mathcal{H})$. For this
reason, $\mathcal{C}$ is a controlization scheme that makes it possible to
approximately implement the controlled time evolution according to $H$, that
is,

$\displaystyle|0\rangle\langle 0|\otimes I_{\mathcal{H}}+|1\rangle\langle
1|\otimes\exp(-iHt).$ (19)

###### Proof.

We adjoin a single qubit $\mathbb{C}^{2}$ to the quantum system $\mathcal{H}$.
The overall Hamiltonian of the joint system $\mathbb{C}^{2}\otimes\mathcal{H}$
is of the form

$\displaystyle I_{2}\otimes H,$ (20)

that is, there is no coupling between $\mathbb{C}^{2}$ and $\mathcal{H}$ and
the control qubit $\mathbb{C}^{2}$ does not have any non-trivial internal
dynamics. We have

$\displaystyle\Lambda(U_{j})$ $\displaystyle=|0\rangle\langle 0|\otimes
U_{j}+|1\rangle\langle 1|\otimes I_{\mathcal{H}}$ (21)

as the control operations in $\mathcal{C}$ for $k=1,\ldots,N$. We obtain

$\displaystyle\mathcal{C}(I_{2}\otimes H)$
$\displaystyle=\sum_{j=1}^{N}\tau_{j}\Lambda(U_{j})(I_{2}\otimes
H)\Lambda(U_{j})^{\dagger}$ (22) $\displaystyle=|0\rangle\langle
0|\otimes\sum_{j=1}^{N}\tau_{j}U_{j}HU_{j}^{\dagger}+|1\rangle\langle
1|\otimes\sum_{j=1}^{N}\tau_{j}H$ (23) $\displaystyle=|0\rangle\langle
0|\otimes\mathcal{D}(H)+|1\rangle\langle 1|\otimes H$ (24)
$\displaystyle=|1\rangle\langle 1|\otimes H,$ (25)

where we used the fact that $\mathcal{D}(H)=0_{\mathcal{H}}$ since
$\mathcal{D}$ is a decoupling scheme for $H$ and $\sum_{j=1}^{N}\tau_{j}=1$ by
convention for decoupling schemes. Finally, observe that time evolution
according to the target Hamiltonian $|1\rangle\langle 1|\otimes H$ gives rise
to the desired controlled time evolution in eq. (19). ∎

######  Remark 2.

The idea of converting a decoupling scheme into a controlization scheme is
sketched but not pursued in [21] because a different model is considered in
that work. Let us briefly explain the main differences. Observe that in our
theorem the operations $U_{j}$ of the decoupling scheme become controlled
operations of the form $\Lambda(U_{j})=|0\rangle\langle 0|\otimes
U_{j}+|1\rangle\langle 1|\otimes I$ of the controlization scheme. Such
controlled operations are specifically disallowed in that work. For instance,
when $\mathcal{H}=(\mathbb{C}^{2})^{\otimes n}$ and the control unitaries
$U_{j}$ are tensor products of single qubit gates, the implementation of
$\Lambda(U_{j})$ would require two-qubit gates. While bang-bang control of
single qubit operations is often a valid assumption, bang-bang control of two-
qubit operations is more problematic.

Therefore, that work considers the following model.333Both to simplify the
presentation and to facilitate the comparison with our model, we describe here
the model of [21] as having a _single_ qubit as control, whereas that work
considers more generally a _multi-qubit_ register. The joint quantum system is
$\mathbb{C}^{2}\otimes\mathcal{H}$, where $\mathbb{C}^{2}$ is the ancilla
qubit and $\mathcal{H}$ is an $n$-qubit system. The former is called the
control qubit and the latter the target register. The starting Hamiltonian
$H_{\mathrm{init}}$ is a pair-interaction Hamiltonian of the form

$\displaystyle H_{\mathrm{init}}$ $\displaystyle=H_{\mathbb{C}^{2}}\otimes
I_{\mathcal{H}}+H_{\mathbb{C}^{2},\mathcal{H}}+I_{2}\otimes H.$ (26)

The term $H_{\mathbb{C}^{2}}$ is the Hamiltonian of the control qubit and the
term $H_{\mathbb{C}^{2},\mathcal{H}}$ specifies the pair-interactions444The
pair-interactions between the control qubit and the qubits of the target
register need to be non-trivial to enable one to construct a simulation scheme
based on average Hamiltonian theory. between the control qubit and the qubits
of the target register. Both terms are assumed to be known. The term $H$ is
the Hamiltonian of the target register. It is an unknown pair-interaction
Hamiltonian, and the task is to controlize it. More precisely, that work
constructs a simulation scheme $\mathcal{C}$ using only single qubit
operations such that

$\displaystyle\mathcal{S}(H_{\mathrm{init}})$ $\displaystyle=\sigma_{z}\otimes
H.$ (27)

Note that having $\sigma_{z}\otimes H$ is equivalent to having
$|1\rangle\langle 1|\otimes H$.

##  4 Decoupling schemes based on orthogonal arrays

The previous section established the connection between controlization and
decoupling schemes without making any assumptions about the quantum system
$\mathcal{H}$ and the Hamiltonian $H\in\mathcal{L}(H)$.

In this section we consider $k$-local Hamiltonians acting on $n$ qudits, that
is, $\mathcal{H}=(\mathbb{C}^{d})^{\otimes n}$. We review some known
constructions of decoupling schemes based on the combinatorial concept of
orthogonal arrays, which then give rise to controlization schemes via Theorem
1. Orthogonal arrays, which were first used for decoupling two-local qubit
Hamiltonians in [38], are defined as follows [20]:

######  Definition 2 (Orthogonal arrays).

An $N\times n$ array $M$ with entries from a finite alphabet $\mathcal{A}$ is
an $OA(N,n,s,t)$ orthogonal array with $s=|\mathcal{A}|$ levels and strength
$t$ if and only if each $N\times t$ subarray contains each $t$-tuple of
$\mathcal{A}$ as a row exactly $\lambda=N/s^{t}$ times.

The following result is known for the situation where the quantum system is
$\mathcal{H}=(\mathbb{C}^{d})^{\otimes n}$, that is, it consists of $n$ qudits
and the system Hamiltonian $H\in\mathcal{L}(\mathcal{H})$ is an arbitrary
$2$-local operator, that is,

$\displaystyle H$ $\displaystyle=\sum_{k\ell}H_{k\ell},$ (28)

where $H_{k\ell}$ are (traceless) Hermitian operators acting only the $k$th
and $\ell$th qudits for $1\leq k<\ell\leq n$:

######  Theorem 2 (OA $\Rightarrow$ decoupling scheme, [45], Theorem 1).

Let $\mathcal{A}$ be the finite alphabet $\\{1,\ldots,d^{2}\\}$. Then, any
orthogonal array with parameters $OA(N,n,d^{2},2)$ over $\mathcal{A}$ can be
used to define a decoupling scheme $\mathcal{D}$ that annihilates any
$2$-local Hamiltonian on $n$ qudits. The number of local operations used in
this scheme is given by $N$.

The unitaries $U_{j}$ for $j=1,\ldots,N$ in the OA-based decoupling schemes
are local operations of the form $U_{j}\in\mathcal{G}^{\otimes n}$, where
$\mathcal{G}$ is a finite subgroup of the unitary group $\mathcal{U}(d)$. The
times $\tau_{j}$ are all equal to $\frac{1}{N}$. For the sake of completeness,
let us briefly explain the correspondence between the entries of orthogonal
array and the control operations. We may always assume that the first unitary
$U_{1}$ is the identity operator. Next, define the $d^{2}$ many generalized
Pauli matrices

$\displaystyle X^{a}Z^{b}\in\mathcal{G}$ (29)

for $a,b\in\\{0,\ldots,d-1\\}$, where

$\displaystyle X$ $\displaystyle=\sum_{x=0}^{d-1}|x+1\mbox{ mod
}d\rangle\langle x|$ (30) $\displaystyle Z$
$\displaystyle=\mathrm{diag}(\omega^{0},\omega^{1}\ldots,\omega^{d-1}),$ (31)

where $\omega=\exp(2\pi i/d)$ is a $d$th root of unity. Denote these
generalized Pauli matrices by $P_{1},\ldots,P_{d^{2}}$. This collection is a
special case of a so-called nice error basis555A collection of $d^{2}$ unitary
matrices $P_{1},\ldots,P_{d^{2}}\in\mathcal{U}(d)$ that are orthogonal with
respect to to the trace inner product $\langle
A,B\rangle=\mathrm{Tr}(A^{\dagger}B)/d$ form a so-called unitary error basis.
Any unitary error basis gives rise to a depolarizing channel on
$\mathbb{C}^{d}$. It is called nice error basis if the matrices are indexed by
a finite group and the multiplication of the matrices corresponds to the
multiplication of the group elements up to scalar factors. for
$\mathbb{C}^{d}$ [25, 27]. It can be shown that all $|\mathcal{A}|^{2}$ tensor
products of two generalized Pauli matrices give rise to the depolarizing
channel $\Phi$ on $\mathbb{C}^{d}\times\mathbb{C}^{d}$. More precisely, the
depolarizing channel is given by

$\displaystyle\Phi(\bullet)$
$\displaystyle=\frac{1}{|\mathcal{A}|^{2}}\sum_{(j,j^{\prime})\in\mathcal{A}\times\mathcal{A}}(P_{j}\otimes
P_{j^{\prime}})\bullet(P_{j}\otimes P_{j^{\prime}})^{\dagger}.$ (32)

In particular, any traceless two-qudit Hamiltonian is annihilated by $\Phi$.
This annihilation property can be “lifted” with help of the orthogonal array
as follows. The $j$th control operation $U_{j}$ is the $n$-fold tensor product
of generalized Pauli matrices defined to be

$\displaystyle U_{j}$ $\displaystyle=P_{m_{j1}}\otimes
P_{m_{j2}}\otimes\cdots\otimes P_{m_{jn}}\in\mathcal{G}^{\otimes n},$ (33)

where the entries $(m_{j1},m_{j2},\ldots,m_{jn})\in\mathcal{A}^{n}$ correspond
to the $j$th row of the orthogonal array. It can be shown that the scheme
simultaneously defines a depolarizing channel on any pair of qudits. This is
due to the defining property of orthogonal arrays: all pairs in
$\mathcal{A}\times\mathcal{A}$ appear exactly $\lambda$ times for any two
columns of the orthogonal array.

1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 3 | 3 | 3 | 3 | 4 | 4 | 4 | 4
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 2 | 3 | 4 | 1 | 2 | 3 | 4 | 1 | 2 | 3 | 4 | 1 | 2 | 3 | 4
1 | 2 | 3 | 4 | 4 | 3 | 2 | 1 | 2 | 1 | 4 | 3 | 3 | 4 | 1 | 2
1 | 2 | 3 | 4 | 2 | 1 | 4 | 3 | 3 | 4 | 1 | 2 | 4 | 3 | 2 | 1
1 | 2 | 3 | 4 | 3 | 4 | 1 | 2 | 4 | 3 | 2 | 1 | 2 | 1 | 4 | 3
Figure 1: Orthogonal array $OA(16,5)$ with multiplicity $1$. The table shows
the transpose of the OA.

The construction of a decoupling scheme with the above approach requires an
orthogonal array of appropriate size. For qubits, we consider the alphabet
$\mathcal{A}=\\{1,2,3,4\\}$ corresponding to the Pauli matrices
$\\{I,X,Y,Z\\}$, which gives $s=|\mathcal{A}|=4$. Figure 1 shows an example
orthogonal array, which can be used to decouple an arbitrary pair-interaction
Hamiltonian on up to $n=5$ qubits with $N=16$ control operations (columns of
the orthogonal array can always be removed whenever $n>t$). Explicit
constructions of orthogonal arrays are given in [20], which allows the
creation of the following arrays:

######  Theorem 3 ([20, Theorem 3.20]).

If $s$ is a prime power then an $OA(s^{n},(s^{n}-1)/(s-1),s,2)$ exists
whenever $n\geq 2$.

######  Theorem 4 ([20, Theorem 6.40]).

If $s$ is a power of a prime and $n\geq 2$, then an orthogonal array
$OA(2s^{n},2(s^{n}-1)/(s-1)-1,s,2)$ can be obtained by using difference
schemes.

Applying these theorems in our setting gives orthogonal arrays of the size
shown in Table 1. Note that most, if not all, constructions in [20] require
$s=d^{2}$, and therefore $d$ itself, to be a prime power. Tables of select
orthogonal arrays are available; see for instance [4, 9, 20].

$\gamma=1$ | $OA(16,5,4,2)$ | Follows from Theorem 3
---|---|---
$\gamma=2$ | $OA(32,9,4,2)$ | Follows from Theorem 4
$\gamma=4$ | $OA(64,21,4,2)$ | Follows from Theorem 3
$\gamma=8$ | $OA(128,41,4,2)$ | Follows from Theorem 4
$\gamma=16$ | $OA(256,85,4,2$ | Follows from Theorem 3
$\gamma=32$ | $OA(512,169,4,2)$ | Follows from Theorem 4
Table 1: Example orthogonal array parameters for used in Pauli ($s=4$)
decoupling schemes for two-local Hamiltonians ($t=2$). For instance , based on
the second entry, an arbitrary pair-interaction Hamiltonian on $n=9$ qubits
can be decoupled using $N=32$ control operations, where each control operation
is tensor product of Pauli operators.

We now consider the situation when the qudits are not all coupled to each
other, that is, the interaction graph of the Hamiltonian is assumed to be a
non-complete known graph. We color the vertices of this graph by assigning
each vertex one of a number of different colors. We need the coloring to be
proper meaning that connected vertices always receive different colors. The
chromatic number is $\chi$ is the minimum number of colors required to achieve
a proper coloring. The chromatic number $\chi$ of a partially coupled qudit
Hamiltonian can be significantly smaller than $n$, which is the chromatic
number of the complete graph. This observation enables us to construct more
efficient decoupling schemes since there are no constraints on the control
sequences between qudits with the same color [22]. Thus, it suffices to
construct decoupling scheme of only a fully coupled $\chi$-qudit (and not a
fully coupled $n$-qudit) system, and apply identical control operations to
qudits of the same color. To summarize, the chromatic number becomes the
effective number of qudits when selecting an orthogonal array.

We mention that it is also possible to handle more general scenarios with
orthogonal arrays. The first generalization encompasses Hamiltonians that have
higher locality than $2$. In this case, one has to rely on orthogonal arrays
of strength equal to the locality parameter of the Hamiltonian. The second
generalization includes an uncontrollable bath that couples to groups of
qudits and the task of decoupling now asks to to remove these bath-system
couplings as well. For simplicity, we restricted our attention to the case of
pair-interactions between the qudits of the system without any bath-system
couplings and refer the reader to [45, 46, 44, 37] for these more general
cases.

##  5 Improved controlization

We now discuss how the OA-based constructions can give improved controlization
schemes. This sections considers only qubit Hamiltonians in order to easily
compare with the controlization schemes in Ref. [32]. The latter work
effectively uses a decoupling scheme that averages over conjugations with all
tensor products of Pauli matrices over $n$ qubits. This is convenient because
it gives a controlization scheme for _any_ $n$-qubit Hamiltonian $H$ without
us having to know anything about $H$. A drawback to this is that the target
Hamiltonian $\ket{1}\bra{1}\otimes H$ generating the controlled-evolution ends
up being a sum of $O(4^{n})$ terms. Thus, the controlization scheme cannot be
implemented by a simple Trotter-Suzuki product formula, as we would need to
implement $4^{n}$ terms for each Trotter step. Ref. [32] overcomes this by
using a randomized product formula, specifically the QDRIFT protocol by
Campbell [8]. The QDRIFT algorithm samples terms from the Hamiltonian with
probabilities according to their weights in the decomposition and implements a
randomized product-formula which converges to the target evolution. The
overall cost scales polynomially with the sum of the weights of the individual
terms instead of the total number of terms. Since the total weight is $O(1)$
in the controlization scheme of [32], QDRIFT enables an efficient
implementation.

Our controlization scheme can be implemented _deterministically_ without using
QDRIFT for local Hamiltonians. In general, the OA-based controlization schemes
can be more efficient when we have access to more structural information about
$H$. For an unknown $n$-qubit Hamiltonian $H$ that is $k$-local, one has to
find an $\mathrm{OA}_{\lambda}(N,n,4,k)$ with $N=\lambda 4^{k}$ such that the
multiplicity $\lambda$ is as small as possible. It is shown in [5] that using
certain BCH-codes in the construction of orthogonal arrays from error
correcting codes in [20, Theorem 4.6] yields $OA(N,n,4,k)$ such that
$N=O(n^{k-1})$. In addition to relying on such asymptotic results, one should
always search the literature or use a software package (e.g., [17]) to
identify an OA with the smallest possible $N$ for the concrete application at
hand. The OA-based decoupling scheme then gives a decomposition of the
effective Hamiltonian into $N$ terms. With this, a simple first-order Trotter
formula can implement controlization using $N$ operations to intersperse the
unknown Hamiltonian time evolution.

Our OA-based controlization scheme can also be realized _stochastically_ using
QDRIFT. The scaling of the gate complexity remains the same as in [32] because
it depends only on the sum of the coefficients, which is still $O(1)$.
However, there is an advantage in that we now only have to sample from $N$
terms, which can potentially make the implementation simpler since for most
situations – when the locality parameter $k$ is much smaller than the number
of qubits – the OA parameter $N$ is going to be significantly smaller than
$4^{n}$. Recall that $4^{n}$ is the size of the set of control operations in
Ref. [32] since it samples $n$-fold tensor products uniformly at random.

Finally, let us remark that $n$-qubit $2$-local Hamiltonians with an
interaction graph of degree at most $\Delta$, can always be decoupled and
controlized using an orthogonal array with $N=O(\Delta)$. Each of the $N$
controlled operations can be written as a product of $n$ two-qubit gates.

##  6 Some future research directions

Recall the correspondence between the $U_{j}$ and $V_{j}$ matrices for
decoupling schemes.

$\displaystyle V_{1}$ $\displaystyle=I$ (34) $\displaystyle V_{j}$
$\displaystyle=U_{j}U_{j-1}^{\dagger}\mbox{ for }j=2,\ldots,N$ (35)
$\displaystyle V_{N+1}$ $\displaystyle=U_{N}^{\dagger}$ (36)

The control operations $V_{j}$ that need to be implemented for the decoupling
scheme are of the form $P_{1j}\otimes P_{2j}\otimes\ldots\otimes
P_{nj}\in\mathcal{G}^{\otimes n}$, where $\mathcal{G}$ is a finite subgroup of
$\mathcal{U}(d)$. Let us assume that for each time step $j$ the local
operations can be implemented in parallel. When we convert the decoupling
scheme into a controlization scheme, however, we have to implement the control
operations

$\displaystyle\Lambda(U_{j})$ $\displaystyle=|0\rangle\langle 0|\otimes
P_{1j}\otimes P_{2j}\otimes\ldots\otimes P_{nj}+|1\rangle\langle 1|\otimes
I\otimes I\otimes\ldots\otimes I.$ (37)

Now, it becomes clear that the weight of $U_{j}$ (that is, the number of non-
identity components) impacts resources needed to realize $\Lambda(U_{j})$.
Thus, one possible question is how to construct decoupling schemes such that
the control operations $V_{j}$ have small weights. Here it may be helpful to
consider orthogonal arrays constructing from error correcting codes as in [5].

Another question related to this is how to realize the control gates
$\Lambda(P_{kj})$ acting on the control qubit and the $k$th qudit efficiently
since long-range gates may not always be directly available. Methods such as
those in [6] may help realize these gates.

### Acknowledgements

We would like to thank Patrick Rall for helpful discussions.

## References

  * AAKS [20] Anurag Anshu, Srinivasan Arunachalam, Tomotaka Kuwahara, and Mehdi Soleimanifar.  Sample-efficient learning of interacting quantum systems.  Nature Physics, 17:931 – 935, 2020.
  * AFCB [14] Mateus Araújo, Adrien Feix, Fabio Costa, and Časlav Brukner.  Quantum circuits cannot control unknown operations.  NEW JOURNAL OF PHYSICS, 22-Sep, Vol. 16(arXiv:1309.7976), 2014.
  * BHMT [02] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  In Quantum computation and information, volume 305 of Contemporary Mathematics, pages 53–74. AMS, 2002.
  * BJL [99] Thomas Beth, Deiter Jungnickel, and Hanfried Lenz.  Design Theory.  Encyclopedia of Mathematics and its Applications. Cambridge University Press, 2 edition, 1999.
  * BRW [16] Adam D. Bookatz, Martin Rötteler, and Pawel Wocjan.  Improved bounded-strength decoupling schemes for local Hamiltonians.  IEEE Transactions on Information Theory, 62(5):2881–2894, 2016\.
  * BTW+ [23] Elisa Bäumer, Vinay Tripathi, Derek S. Wang, Patrick Rall, Edward H. Chen, Swarnadeep Majumder, Alireza Seif, and Zlatko K. Minev.  Efficient long-range entanglement using dynamic circuits, 2023.
  * BWV [14] Adam D. Bookatz, Pawel Wocjan, and Lorenza Viola.  Hamiltonian quantum simulation with bounded-strength controls.  New Journal of Physics, 16(4):045021, apr 2014.
  * Cam [19] Earl Campbell.  Random compiler for fast hamiltonian simulation.  Phys. Rev. Lett., 123:070503, Aug 2019.
  * CD [06] Charles J. Colbourn and Jeffrey H. Dinitz, editors.  Handbook of Combinatorial Designs.  Chapman and Hall/CRC, 2006.
  * CE [16] Giulio Chiribella and Daniel Ebler.  Optimal quantum networks and one-shot entropies.  New Journal of Physics, 18(9):093053, sep 2016.
  * CKS [17] Andrew M. Childs, Robin Kothari, and Rolando D. Somma.  Quantum linear systems algorithm with exponentially improved dependence on precision.  SIAM J. Comp., 46:1920, 2017.
  * CS [17] Anirban Narayan Chowdhury and Rolando D. Somma.  Quantum algorithms for Gibbs sampling and hitting-time estimation.  Quantum Information & Computation, 17(1–2):41–64, feb 2017\.
  * CW [12] Andrew M. Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  Quantum Information & Computation, 12(11–12):901–924, nov 2012\.
  * DDTL [24] Zhiyan Ding, Yulong Dong, Yu Tong, and Lin Lin.  Robust ground-state energy estimation under depolarizing noise, 2024.
  * DLT [22] Yulong Dong, Lin Lin, and Yu Tong.  Ground-state preparation and energy estimation on early fault-tolerant quantum computers via quantum eigenvalue transformation of unitary matrices.  PRX Quantum, 3:040305, Oct 2022.
  * DNSM [21] Qingxiuxiong Dong, Shojun Nakayama, Akihito Soeda, and Mio Murao.  Controlled quantum operations and combs, and their applications to universal controllization of divisible unitary operations, 2021.
  * EV [19] Pieter Thijs Eendebak and Alan Roberto Vazquez.  OApackage: A Python package for generation and analysis of orthogonal arrays, optimal designs and conference designs.  Journal of Open Source Software, 4(34):1097, 2019.
  * GSLW [19] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics.  In Proc. of the 51st Annual ACM SIGACT Symp. Theor. Comp., STOC 2019, page 193–204, New York, NY, USA, 2019. Association for Computing Machinery.
  * HHL [09] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Phys. Rev. Lett., 103:150502, Oct 2009.
  * HSS [99] A. S. Hedayat, Neil J. A. Sloane, and John Stufken.  Orthogonal Arrays Theory and Applications.  Springer, 1999.
  * Jan [02] Dominik Janzing.  Quantum algorithm for measuring the energy of $n$ qubits with unknown pair-interactions.  Quantum Information and Computation, 2(3):198–207, 2002.
  * JK [99] Jonathan A. Jones and Emanuel Knill.  Efficient refocusing of one-spin and two-spin interactions for NMR quantum computation.  Journal of Magnetic Resonance, 141(2):322–325, 1999.
  * JWB [02] Dominik Janzing, Pawel Wocjan, and Thomas Beth.  Complexity of decoupling and time reversal for n spins with pair interactions: Arrow of time in quantum control.  Phys. Rev. A, 66:042311, Oct 2002.
  * KGR [23] Oriel Kiss, Michele Grossi, and Alessandro Roggero.  Importance sampling for stochastic quantum simulations.  Quantum, 7:977, April 2023.
  * Kni [96] Emanuel Knill.  Non-binary unitary error bases and quantum codes, 1996.
  * Kot [14] Robin Kothari.  Efficient algorithms in quantum query complexity.  PhD thesis, University of Waterloo, 2014.
  * KR [02] Andreas A. Klappenecker and Martin Rötteler.  Beyond stabilizer codes I: Nice error bases.  IEEE Transactions on Information Theory, 48(8):2392–2395, 2002\.
  * LC [19] Guang Hao Low and Isaac L. Chuang.  Hamiltonian simulation by qubitization.  Quantum, 3:163, 2019.
  * LZH [20] Zhi Li, Liujun Zou, and Timothy H. Hsieh.  Hamiltonian tomography via quantum quench.  Phys. Rev. Lett., 124:160502, Apr 2020.
  * Nav [18] Miguel Navascués.  Resetting uncontrolled quantum systems.  Phys. Rev. X, 8:031008, Jul 2018.
  * NC [00] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2000.
  * OKTM [24] Tatsuki Odake, Hlér Kristjánsson, Philip Taranto, and Mio Murao.  Universal algorithm for transforming Hamiltonian eigenvalues, 2024.
  * PB [14] Shengshi Pang and Todd A. Brun.  Quantum metrology for a general Hamiltonian parameter.  Physical Review A, 90:022117, 2014.
  * [34] David Poulin and Pawel Wocjan.  Preparing ground states of quantum many-body systems on a quantum computer.  Phys. Rev. Lett., 102:130503, 2009.
  * [35] David Poulin and Pawel Wocjan.  Sampling from the thermal quantum Gibbs state and evaluating partition functions with a quantum computer.  Phys. Rev. Lett., 103:220502, 2009.
  * QDS+ [19] Marco Túlio Quintino, Qingxiuxiong Dong, Atsushi Shimbo, Akihito Soeda, and Mio Murao.  Reversing unknown quantum transformations: Universal quantum circuit for inverting general unitary operations.  Phys. Rev. Lett., 123:210502, Nov 2019.
  * RW [06] Martin Rötteler and Pawel Wocjan.  Equivalence of decoupling schemes and orthogonal arrays.  IEEE Transactions on Information Theory, 52(9):4171–4181, 2006\.
  * SM [01] Marcus Stollsteimer and Günter Mahler.  Suppression of arbitrary internal coupling in a quantum register.  Phys. Rev. A, 64:052301, Oct 2001.
  * STCA [23] Thais L. Silva, Márcio M. Taddei, Stefano Carrazza, and Leandro Aolita.  Fragmented imaginary-time evolution for early-stage quantum signal processors.  Scientific Reports, 13(1):18258, 2023.
  * Suz [90] Masuo Suzuki.  Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations.  Physics Letters A, 146(6):319–323, 1990.
  * Suz [91] Masuo Suzuki.  General theory of fractal path integrals with applications to many-body theories and statistical physics.  Journal of Mathematical Physics, 32(2):400–407, February 1991.
  * TDN [20] David Trillo, Benjamin Dive, and Miguel Navascués.  Translating Uncontrolled Systems in Time.  Quantum, 4:374, December 2020.
  * WGFC [14] Nathan Wiebe, Christopher Granade, Christopher Ferrie, and D. G. Cory.  Hamiltonian learning and certification using quantum resources.  Phys. Rev. Lett., 112:190501, May 2014.
  * Woc [03] Pawel Wocjan.  Computational Power of Hamiltonians in Quantum Computing.  PhD thesis, University of Karlsruhe, 2003.  https://schneider128k.github.io/dissertation.pdf.
  * [45] Pawel Wocjan, Martin Rötteler, Dominik Janzing, and Thomas Beth.  Simulating hamiltonians in quantum networks: Efficient schemes and complexity bounds.  Phys. Rev. A, 65:042309, Mar 2002.
  * [46] Pawel Wocjan, Martin Rötteler, Dominik Janzing, and Thomas Beth.  Universal simulation of Hamiltonians using a finite set of control operations.  Quantum Information and Computation, 2(2):133–150, 2002.


# On Some Properties of the Trigamma Function

Kwara Nantomah1,∗, Gregory Abe-I-Kpeng2 and Sunday Sandow3 1Department of
Mathematics, School of Mathematical Sciences, C. K. Tedam University of
Technology and Applied Sciences, P. O. Box 24, Navrongo, Upper-East Region,
Ghana<EMAIL_ADDRESS>2Department of Industrial Mathematics, School
of Mathematical Sciences, C. K. Tedam University of Technology and Applied
Sciences, P. O. Box 24, Navrongo, Upper-East Region, Ghana.
<EMAIL_ADDRESS>3Department of Mathematics, School of Mathematical
Sciences, C. K. Tedam University of Technology and Applied Sciences, P. O. Box
24, Navrongo, Upper-East Region, Ghana<EMAIL_ADDRESS>
(Date: Received: xxxxxx; Revised: xxxxxx; Accepted: xxxxxx.
∗ Corresponding author)

###### Abstract.

In 1974, Gautschi proved an intriguing inequality involving the gamma function
$\Gamma$. Precisely, he proved that, for $z>0$, the harmonic mean of
$\Gamma(z)$ and $\Gamma(1/z)$ can never be less than 1. In 2017, Alzer and
Jameson extended this result to the digamma function $\psi$ by proving that,
for $z>0$, the harmonic mean of $\psi(z)$ and $\psi(1/z)$ can never be less
than $-\gamma$ where $\gamma$ is the Euler-Mascheroni constant. In this paper,
our goal is to extend the results to the trigamma function $\psi^{\prime}$. We
prove among other things that, for $z>0$, the harmonic mean of
$\psi^{\prime}(z)$ and $\psi^{\prime}(1/z)$ can never be greater than
$\pi^{2}/6$.

###### Key words and phrases:

Gamma function; digamma function; trigamma function; harmonic mean inequality

###### 2010 Mathematics Subject Classification:

33B15, 26A48, 26A51, 41A29

##  1\. Introduction

The classical gamma function which is an extension of the factorial function
is frequently defined as

$\Gamma(z)=\int_{0}^{\infty}r^{z-1}e^{-r}dr$

for $z>0$. Closely connected to the gamma function is the the digamma (or psi)
function is which is defined as

$\displaystyle\psi(z)=\frac{d}{dz}\ln\Gamma(z)$
$\displaystyle=-\gamma+\int_{0}^{\infty}\frac{e^{-r}-e^{-zr}}{1-e^{-r}}dr,$
(1)
$\displaystyle=\int_{0}^{\infty}\left(\frac{e^{-r}}{r}-\frac{e^{-zr}}{1-e^{-r}}\right)dr,$
(2) $\displaystyle=-\gamma-\frac{1}{z}+\sum_{n=1}^{\infty}\frac{z}{n(n+z)},$
(3)

where $\gamma$ is the Euler-Mascheroni constant. Derivatives of the digamma
function which are called polygamma functions are defined as

$\displaystyle\psi^{(c)}(z)$
$\displaystyle=(-1)^{c+1}\int_{0}^{\infty}\frac{r^{c}e^{-zr}}{1-e^{-r}}dr,$
(4) $\displaystyle=(-1)^{c+1}\sum_{n=0}^{\infty}\frac{c!}{(n+z)^{c+1}},$ (5)

for $z>0$ and $c\in\mathbb{N}$. The particular case $\psi^{\prime}(z)$ is what
is referred to as the trigamma function. Also, it is well known in the
literature that the integral

$\frac{c!}{z^{c+1}}=\int_{0}^{\infty}r^{c}e^{-zr}\,dr$ (6)

holds for $z>0$ and $c\in\mathbb{N}_{0}$.

In 1974, Gautschi [11] presented an elegant inequality involving the gamma
function. Precisely, he proved that, for $z>0$, the harmonic mean of
$\Gamma(z)$ and $\Gamma(1/z)$ is at least 1. That is,

$\frac{2\Gamma(z)\Gamma(1/z)}{\Gamma(z)+\Gamma(1/z)}\geq 1,$ (7)

for $z>0$ and with equality when $z=1$. As a direct consequence of (7), the
inequalities

$\Gamma(z)+\Gamma(1/z)\geq 2$ (8)

and

$\Gamma(z)\Gamma(1/z)\geq 1$ (9)

are obtained for $z>0$. Attributing to the importance of this inequality, some
refinements and extensions have been investigated [1, 2, 3, 4, 5, 6, 12, 13].

In 2017, Alzer and Jameson [8] established a striking companion of (7) which
involves the digamma function $\psi(z)$. They established that the inequality

$\frac{2\psi(z)\psi(1/z)}{\psi(z)+\psi(1/z)}\geq-\gamma$ (10)

holds for $z>0$ and with equality when $z=1$. Thereafter, Alzer [7] refined
(10) by proving that

$\frac{2\psi(z)\psi(1/z)}{\psi(z)+\psi(1/z)}\geq-\gamma\frac{2z}{z^{2}+1}$
(11)

holds for $z>0$ and with equality when $z=1$.

In 2018, Yin et al. [23] extended inequality (10) to the $k$-analogue of the
digamma function by establishing that

$\frac{2\psi_{k}(z)\psi_{k}(1/z)}{\psi_{k}(z)+\psi_{k}(1/z)}\geq\frac{\ln^{2}k+\gamma^{2}-2(\gamma+1)\ln
k}{k\left[\ln k+\psi(1/k)\right]}$ (12)

for $z>0$ and $\frac{1}{\sqrt[3]{3}}\leq k\leq 1$.

In 2020, Yildirim [22] improved on the inequality (12) by establishing that

$\frac{2\psi_{k}(z)\psi_{k}(1/z)}{\psi_{k}(z)+\psi_{k}(1/z)}\geq\psi_{k}(1)$
(13)

for $z>0$ and $k>0$. When $k=1$, inequalities (12) and (13) both return to
inequality (10).

In 2021, Bouali [10] extended inequalities (7) and (10) to the $q$-analogues
of the gamma and digamma functions by proving that

$\frac{2\Gamma_{q}(z)\Gamma_{q}(1/z)}{\Gamma_{q}(z)+\Gamma_{q}(1/z)}\geq 1$
(14)

for $z>0$ and

$\frac{2\psi_{q}(z)\psi_{q}(1/z)}{\psi_{q}(z)+\psi_{q}(1/z)}\geq\psi_{q}(1)$
(15)

for $z>0$ and $q\in(0,p_{0})$, where $p_{0}\simeq 3.239945$.

For similar results involving other special functions, one may refer to the
works [14, 15, 16, 17, 18, 19].

In the present investigation, our goal is to extend the results of Alzer and
Jameson [8] to the trigamma function $\psi^{\prime}$. Specifically, we prove
among other things that, for $z>0$, the harmonic mean of $\psi^{\prime}(z)$
and $\psi^{\prime}(1/z)$ can never be greater than $\pi^{2}/6$. We present our
results in Section 2. In order to establish our results, we require the
following preliminary definitions and lemmas.

######  Definition 1.1 ([20]).

A function $H:\mathcal{I}\subseteq\mathbb{R}^{+}\rightarrow\mathbb{R}$ is
referred to as GG-convex if

$H(x^{1-k}y^{k})\leq H(x)^{1-k}H(y)^{k}$ (16)

for all $x,y\in\mathcal{I}$ and $k\in[0,1]$. If the inequality in (16) is
reversed, then $H$ is said to be GG-concave.

######  Definition 1.2 ([20]).

A function $H:\mathcal{I}\subseteq\mathbb{R}^{+}\rightarrow\mathbb{R}$ is
referred to as GA-convex if

$H(x^{1-k}y^{k})\leq(1-k)H(x)+kH(y)$ (17)

for all $x,y\in\mathcal{I}$ and $k\in[0,1]$. If the inequality in (17) is
reversed, then $H$ is said to be GA-concave.

######  Lemma 1.3 ([20]).

A function $H:\mathcal{I}\subseteq\mathbb{R}^{+}\rightarrow\mathbb{R}$ is GG-
convex (or GG-concave) if and only if $\frac{zH^{\prime}(z)}{H(z)}$ is
increasing (or decreasing) on $\mathcal{I}$ respectively.

######  Lemma 1.4 ([24]).

A function $H:\mathcal{I}\subseteq\mathbb{R}^{+}\rightarrow\mathbb{R}$ is GA-
convex if and only if

$H^{\prime}(z)+zH^{\prime\prime}(z)\geq 0$ (18)

for all $z\in\mathcal{I}$. The function $H$ is said to be GA-concave if and
only if the inequality in (18) is reversed.

The following lemma is well known in the literature as the convolution theorem
for Laplace transforms.

######  Lemma 1.5.

Let $f(r)$ and $g(r)$ be any two functions with convolution $f\ast
g=\int_{0}^{r}f(r-s)g(s)\,ds$. Then the Laplace transform of the convolution
is given as

$\mathcal{L}\left\\{f\ast
g\right\\}=\mathcal{L}\left\\{f\right\\}\mathcal{L}\left\\{g\right\\}.$

In other words,

$\int_{0}^{\infty}\left[\int_{0}^{r}f(r-s)g(s)\,ds\right]e^{-zr}\,dr=\int_{0}^{\infty}f(r)e^{-zr}\,dr\int_{0}^{\infty}g(r)e^{-zr}\,dr.$
(19)

######  Lemma 1.6 ([21]).

Let $-\infty\leq u<v\leq\infty$ and $p$ and $q$ be continuous functions that
are differentiable on $(u,v)$, with $p(u+)=q(u+)=0$ or $p(v-)=q(v-)=0$.
Suppose that $q(z)$ and $q^{\prime}(z)$ are nonzero for all $z\in(u,v)$. If
$\frac{p^{\prime}(z)}{q^{\prime}(z)}$ is increasing (or decreasing) on
$(u,v)$, then $\frac{p(x)}{q(x)}$ is also increasing (or decreasing) on
$(u,v)$.

##  2\. Results

######  Theorem 2.1.

The function $\psi^{\prime}(z)$ is GG-convex on $(0,\infty)$. In other words,

$\psi^{\prime}(x^{1-k}y^{k})\leq\left[\psi^{\prime}(x)\right]^{1-k}\left[\psi^{\prime}(y)\right]^{k}$
(20)

is satisfied for $x>0$, $y>0$ and $k\in[0,1]$.

###### Proof.

As a result of Lemma 1.3, it suffices to show that the function
$z\frac{\psi^{\prime\prime}(z)}{\psi^{\prime}(z)}$ is increasing on
$(0,\infty)$ and this follows from Lemma 2 of [7]. ∎

######  Corollary 2.2.

The inequality

$\psi^{\prime}(z)\psi^{\prime}(1/z)\geq\left(\frac{\pi^{2}}{6}\right)^{2}$
(21)

holds for $z\in(0,\infty)$ and with equality when $z=1$.

###### Proof.

By letting $x=z$, $y=1/z$ and $k=\frac{1}{2}$ in (20), we obtain

$\sqrt{\psi^{\prime}(z)\psi^{\prime}(1/z)}\geq\psi^{\prime}(1)=\frac{\pi^{2}}{6}$

which gives the desired result. ∎

######  Lemma 2.3.

For $r>0$, we have

$0<\frac{re^{-r}}{1-e^{-r}}<1.$ (22)

###### Proof.

By direct computation, we obtain

$\mathcal{B}(r)=\frac{re^{-r}}{1-e^{-r}}=\frac{p_{1}(r)}{q_{1}(r)}$

where $p_{1}(r)=re^{-r}$, $q_{1}(r)=1-e^{-r}$ and $p_{1}(0+)=q_{1}(0+)=0$.
Then

$\frac{p^{\prime}_{1}(r)}{q^{\prime}_{1}(r)}=1-r$

and then

$\left(\frac{p^{\prime}_{1}(r)}{q^{\prime}_{1}(r)}\right)^{\prime}=-1<0.$

Thus, $\frac{p^{\prime}_{1}(r)}{q^{\prime}_{1}(r)}$ is decreasing and as a
result of Lemma 1.6, the function $\mathcal{B}(r)$ is also decreasing. Hence

$0=\lim_{r\to\infty}\mathcal{B}(r)<\mathcal{B}(r)<\lim_{r\to
0+}\mathcal{B}(r)=1$

which completes the proof. ∎

######  Theorem 2.4.

The function $\psi^{\prime}(z)$ is GA-convex on $(0,\infty)$. In other words,

$\psi^{\prime}(x^{1-k}y^{k})\leq(1-k)\psi^{\prime}(x)+k\psi^{\prime}(y)$ (23)

is satisfied for $x>0$, $y>0$ and $k\in[0,1]$.

###### Proof.

As a result of Lemma 1.4, it suffices to show that

$\phi(z)=\psi^{\prime\prime}(z)+z\psi^{\prime\prime\prime}(z)\geq 0$ (24)

for $z\in(0,\infty)$. By applying (4), (6) and Lemma 1.5, we obtain

$\displaystyle\frac{\phi(z)}{z}$
$\displaystyle=\frac{1}{z}\psi^{\prime\prime}(z)+\psi^{\prime\prime\prime}(z)$
$\displaystyle=-\int_{0}^{\infty}e^{-zr}dr\int_{0}^{\infty}\frac{r^{2}e^{-zr}}{1-e^{-r}}dr+\int_{0}^{\infty}\frac{r^{3}e^{-zr}}{1-e^{-r}}dr$
$\displaystyle=-\int_{0}^{\infty}\left[\int_{0}^{r}\frac{s^{2}}{1-e^{-s}}ds\right]e^{-zr}dr+\int_{0}^{\infty}\frac{r^{3}e^{-zr}}{1-e^{-r}}dr$
$\displaystyle=\int_{0}^{\infty}\mathcal{A}(r)e^{-zr}dr$

where

$\mathcal{A}(r)=\frac{r^{3}}{1-e^{-r}}-\int_{0}^{r}\frac{s^{2}}{1-e^{-s}}ds.$

Then by direct computations and as a result of (22), we have

$\displaystyle\mathcal{A}^{\prime}(r)$
$\displaystyle=\frac{3r^{2}}{1-e^{-r}}-\frac{r^{3}e^{-r}}{(1-e^{-r})^{2}}-\frac{r^{2}}{1-e^{-r}}$
$\displaystyle=\frac{r^{2}}{1-e^{-r}}\left[2-\frac{re^{-r}}{1-e^{-r}}\right]\geq
0.$

Hence $\mathcal{A}(r)$ is increasing and this implies that

$\mathcal{A}(r)\geq\lim_{r\to 0+}\mathcal{A}(r)=0.$

Therefore, $\phi(z)\geq 0$ which completes the proof. ∎

######  Remark 2.5.

Inequality (24) implies that the function $z\psi^{\prime\prime}(z)$ is
increasing.

######  Corollary 2.6.

The inequality

$\psi^{\prime}(z)+\psi^{\prime}(1/z)\geq\frac{\pi^{2}}{3}$ (25)

holds for $z\in(0,\infty)$ and with equality when $z=1$.

###### Proof.

By letting $x=z$, $y=1/z$ and $k=\frac{1}{2}$ in (23), we obtain

$\frac{\psi^{\prime}(z)}{2}+\frac{\psi^{\prime}(1/z)}{2}\geq\psi^{\prime}(1)=\frac{\pi^{2}}{6}$

which gives the desired result. ∎

######  Lemma 2.7 ([9]).

For $z>0$, the inequality

$\psi^{\prime}(z)\psi^{\prime\prime\prime}(z)-2\left[\psi^{\prime\prime}(z)\right]^{2}\leq
0$ (26)

is satisfied.

######  Lemma 2.8.

For $z>0$, the function

$F(z)=\frac{z\psi^{\prime\prime}(z)}{[\psi^{\prime}(z)]^{2}}$ (27)

is decreasing.

###### Proof.

By applying Lemma 2.7, we obtain

$\displaystyle[\psi^{\prime}(z)]^{3}F^{\prime}(z)$
$\displaystyle=\psi^{\prime}(z)\psi^{\prime\prime}(z)+z\psi^{\prime}(z)\psi^{\prime\prime\prime}(z)-2z[\psi^{\prime\prime}(z)]^{2}$
$\displaystyle=\psi^{\prime}(z)\psi^{\prime\prime}(z)+z\left[\psi^{\prime}(z)\psi^{\prime\prime\prime}(z)-2\left[\psi^{\prime\prime}(z)\right]^{2}\right]$
$\displaystyle<0.$

Hence $F^{\prime}(z)<0$ which completes the proof. ∎

######  Theorem 2.9.

For $z>0$, the inequality

$\frac{2\psi^{\prime}(z)\psi^{\prime}(1/z)}{\psi^{\prime}(z)+\psi^{\prime}(1/z)}\leq\frac{\pi^{2}}{6}$
(28)

holds and equality is attained if $z=1$.

###### Proof.

The case for $z=1$ is apparent. For this reason, we only prove the case for
$z\in(0,1)\cup(1,\infty)$. Let

$\mathcal{K}(z)=\frac{2\psi^{\prime}(z)\psi^{\prime}(1/z)}{\psi^{\prime}(z)+\psi^{\prime}(1/z)}\quad\text{and}\quad\beta(z)=\ln\mathcal{K}(z)$

for $z\in(0,1)\cup(1,\infty)$. Then direct calculations gives

$\beta^{\prime}(z)=\frac{\psi^{\prime\prime}(z)}{\psi^{\prime}(z)}-\frac{1}{z^{2}}\frac{\psi^{\prime\prime}(1/z)}{\psi^{\prime}(1/z)}-\frac{\psi^{\prime\prime}(z)-\frac{1}{z^{2}}\psi^{\prime\prime}(1/z)}{\psi^{\prime}(z)+\psi^{\prime}(1/z)}$

which implies that

$z\left[\psi^{\prime}(z)+\psi^{\prime}(1/z)\right]\beta^{\prime}(z)=z\frac{\psi^{\prime\prime}(z)}{\psi^{\prime}(z)}\psi^{\prime}(1/z)-\frac{1}{z}\frac{\psi^{\prime\prime}(1/z)}{\psi^{\prime}(1/z)}\psi^{\prime}(z).$

This further gives rise to

$\displaystyle
z\left[\frac{1}{\psi^{\prime}(z)}+\frac{1}{\psi^{\prime}(1/z)}\right]\beta^{\prime}(z)$
$\displaystyle=z\frac{\psi^{\prime\prime}(z)}{[\psi^{\prime}(z)]^{2}}-\frac{1}{z}\frac{\psi^{\prime\prime}(1/z)}{[\psi^{\prime}(1/z)]^{2}}$
$\displaystyle:=T(z).$

As a result of Lemma 2.8, we conclude that $T(z)>0$ if $z\in(0,1)$ and
$T(z)<0$ if $z\in(1,\infty)$. Thus, $\beta(z)$ is increasing on $(0,1)$ and
decreasing on $(1,\infty)$. Accordingly, $\mathcal{K}(z)$ is increasing on
$(0,1)$ and decreasing on $(1,\infty)$. Therefore, on both intervals, we have

$\mathcal{K}(z)<\lim_{z\to
1}\mathcal{K}(z)=\psi^{\prime}(1)=\frac{\pi^{2}}{6}$

completing the proof. ∎

## References

  * [1] H. Alzer, A harmonic mean inequality for the gamma function, J. Comp. Appl. Math., 87 (1997), 195-198.
  * [2] H. Alzer, Inequalities for the gamma function, Proc. Amer. Math. Soc., 128 (1999), 141-147.
  * [3] H. Alzer, On a gamma function inequality of Gautschi, Proc. Edinburgh Math. Soc., 45 (2002), 589-600.
  * [4] H. Alzer, On Gautschi’s harmonic mean inequality for the gamma function, J. Comp. Appl. Math., 157 (2003), 243-249.
  * [5] H. Alzer, Inequalities involving $\Gamma(x)$ and $\Gamma(1/x)$, J. Comp. Appl. Math., 192 (2006), 460-480.
  * [6] H. Alzer, Gamma function inequalities, Numer. Algor., 49 (2008), 53-84.
  * [7] H. Alzer, A Mean Value Inequality for the Digamma Function, Rendiconti Sem. Mat. Univ. Pol. Torino, 75(2)(2017), 19-25.
  * [8] H. Alzer and G. Jameson, A harmonic mean inequality for the digamma function and related results, Rend. Sem. Mat. Univ. Padova., 137 (2017), 203-209.
  * [9] H. Alzer and J. Wells, Inequalities for the Polygamma Functions, SIAM J. Math. Anal., 29(6)(1998), 1459-1466.
  * [10] M. Bouali, A harmonic mean inequality for the $q$-gamma and $q$-digamma functions, Filomat 35(12)(2021), 4105-4119.
  * [11] W. Gautschi, A harmonic mean inequality for the gamma function, SIAM J. Math. Anal., 5(1974), 278-281.
  * [12] C. Giordano and A. Laforgia, Inequalities and monotonicity properties for the gamma function, J. Comp. Appl. Math., 133 (2001), 387-396.
  * [13] G. J. O. Jameson and T. P. Jameson, An inequality for the gamma function conjectured by D. Kershaw, J. Math. Ineq., 6 (2012), 175-181.
  * [14] L. Matejicka, Proof of a Conjecture On Nielsen’s $\beta$-Function, Probl. Anal. Issues Anal., 8(26), (2019), 105-111.
  * [15] K. Nantomah, Certain Properties of the Nielsen’s $\beta$-Function, Bull. Int. Math. Virtual Inst., 9(2019), 263-269.
  * [16] K. Nantomah, Harmonic Mean Inequalities for Hyperbolic Functions, Earthline J. Math. Sci., 6(1)(2021), 117-129.
  * [17] K. Nantomah, A Harmonic Mean Inequality for the Exponential Integral Function, Int. J. Appl. Math., 34(4)(2021), 647-652.
  * [18] K. Nantomah, A Harmonic Mean Inequality Concerning the Generalized Exponential Integral Function, Adv. Math. Sci. J., 10(9)(2021), 3227-3231.
  * [19] K. Nantomah, Degenerate Exponential Integral Function and its Properties, Arab J. Math. Sci., Published Ahead-of-Print: https://doi.org/10.1108/AJMS-09-2021-0230.
  * [20] C. P. Niculescu, Convexity according to the geometric mean, Math. Inequal. Appl., 2(2)(2000), 155-167.
  * [21] I. Pinelis, L’hospital type Rules for Monotonicity, with Applications, J. Inequal. Pure Appl. Math., 3(1)(2002), Art No. 5, 5 pages.
  * [22] E. Yildirim, Monotonicity Properties on $k$-Digamma Function and its Related Inequalities, J. Math. Inequal., 14(1)(2020), 161-173.
  * [23] L. Yin, L-G. Huang, X-L. Lin and Y-L. Wang, Monotonicity, concavity, and inequalities related to the generalized digamma function, Adv. Difference Equ., (2018) 2018:246.
  * [24] X-M. Zhang, Y-M. Chu and X-H. Zhang, The Hermite-Hadamard Type Inequality of GA-Convex Functions and Its Application, J. Inequal. Appl., 2010(2010), Article ID 507560, 11 pages.