Optimizing dual-energy CT technique for iodine-based contrast-to-noise ratio
Fatma Terzioglu1, Emil Y. Sidky2, Jp Phillips2, Ingrid S. Reiser2, Guillaume
Bal3, and Xiaochuan Pan2
1Department of Mathematics, North Carolina State University, 2311 Stinson Dr.,
Raleigh, NC 27695, USA
3Departments of Statistics and Mathematics, The University of Chicago, 5747 S.
Ellis Ave., Chicago, IL 60637, USA
2Department of Radiology, The University of Chicago, 5841 S. Maryland Ave.,
Chicago, IL 60637, USA

Author to whom correspondence should be addressed. email<EMAIL_ADDRESS>

###### Abstract

Purpose: The goal of this study is to propose a systematic method for
determining the optimal x-ray tube settings/energy windows and fluence for
minimal noise and maximum CNR in material density images obtained from dual-
energy CT (DECT) scans by fixing the subject size and the total radiation
dose.
Methods: The noise propagation in the process of sinogram and image
reconstruction from DECT measurements is analyzed. The main objects of the
study are the pixel variances for the sinogram and monochromatic image and the
contrast-to-noise ratio (CNR), which were shown to depend on the Jacobian
matrix of the sinograms-to-DECT measurements map.

Analytic estimates for the sinogram and monochromatic image pixel variances
and the CNR as functions of tube potentials, fluence, and virtual
monochromatic image (VMI) energy are derived, and then used in a phantom
experiment as an objective function for optimizing the tube settings to
minimize the image noise and maximize the CNR.
Results: A non-trivial example that shows the existence of singular solutions
to the inversion of sinograms-to-DECT measurements map was presented.
Additionally, the optimal VMI energy for maximal CNR was determined. The
optimal energy VMI was found to be the least noisy monochromatic image
synthesized from the iodine and water density images, and it was shown that
using more general weights in combining the two images linearly does not
improve image quality.

When the x-ray beam filter material was fixed at 2mm of Aluminum and the
photon fluence for low and high kV scans were considered equal, the tube
potential pair of 60/120 kV led to the maximal CNR in the VMI formed at energy
55 KeV.
Conclusions: Optimizing DECT scan parameters to maximize the CNR can be done
in a systematic way. Also choosing the parameters that maximize the Jacobian
determinant over the sinogram domain would lead to more stable reconstructions
due to the reduced amplification of the measurement noise. Since the values of
the Jacobian determinant depend strongly on the imaging task, careful
consideration of all of the relevant factors is needed when implementing the
proposed framework.

###### Contents

  1. I. Introduction
  2. II. Methods
    1. II.A. Scan configuration, DECT technique, and noise simulation
    2. II.B. Jacobian of the sinogram-to-photon-counts transformation and uniqueness of reconstructions
    3. II.C. Noise propagation based on linearization of the inverse sinogram-to-photon-counts transformation
    4. II.D. Computation of pixel variances and the CNR of a monochromatic image
    5. II.E. Minimization of mean pixel variance
  3. III. Results
    1. III.A. Experimental setup
    2. III.B. Uniqueness of the reconstructed sinograms
    3. III.C. Optimization of the iodine CNR in virtual monochromatic images
  4. IV. Discussion and conclusion

##  I. Introduction

Dual-energy CT (DECT) systems enable the simultaneous acquisition of two
spectral measurements to identify different materials within the scanned
object. DECT has been demonstrated to outperform single-energy CT in terms of
image quality and contrast-to-noise ratio (CNR), allowing for reduced
radiation exposure and contrast agent concentration while maintaining image
quality 1, 2, 3, 4, 5.

The concept of DECT was introduced by Hounsfield in the 1970s 6, and the
mathematical framework for pre-reconstruction processing of DECT data was
developed by Alvarez and Macovski in a seminal paper in 1976 7. Their approach
to material decomposition is based on the assumption that energy-dependent
attenuation coefficients of chemical compounds can be approximated by a linear
combination of elemental mass attenuation maps weighted by the partial density
of each element. This reduces sinogram reconstruction to solving a nonlinear
system of equations for each sinogram value.

Due to the nonlinearity of the DECT measurement model, the uniqueness of
reconstructed sinograms from DECT measurements is not guaranteed. Levine first
provided an example of non-unique solutions in DECT for a material basis of
water and bone using spectra with three discrete photon energies 8. Alvarez
analyzed the non-invertibility of the sinogram to DECT measurements mapping by
studying the Jacobian determinant 9. In general, the nonvanishing of the
Jacobian determinant only guarantees local uniqueness (a.k.a. injectivity),
requiring additional constraints for global uniqueness. Bal and Terzioglu10
presented sufficient analytic criteria for the global injectivity of multi-
energy CT (MECT) measurement map for a given number of materials and equal
number of energy measurements. In the case of DECT, they showed that the
nonvanishing of the Jacobian determinant of the sinogram to DECT measurements
map throughout its domain is sufficient to ensure global uniqueness. They also
demonstrated how the choice of basis materials and x-ray spectra influences
the Jacobian determinant values and, consequently, the invertibility. In this
paper, we showcase the occurrence of nonuniqueness when the Jacobian
determinant vanishes within the rectangular region encompassing all possible
sinogram values, providing a clear example of two isolated solutions
corresponding to the same DECT measurement pair.

It was shown by Bal et al. 11 that the stability of the inversion of the
mapping of sinograms to MECT measurements is improved by choosing the x-ray
spectra that maximize Jacobian determinant values of the sinogram to DECT
measurements map. In this paper, we propose a systematic method for
determining optimal tube settings and fluence that minimize noise and maximize
CNR in material density images obtained from DECT scans while keeping subject
size and total radiation dose fixed. To achieve this, we consider a noise
model based on a compound Poisson process for the DECT measurements and
analyze the noise propagation from DECT measurements to material density
images by linearizing the inverse sinogram-to-photon-counts transformation. We
derive analytic expressions for sinogram and monochromatic image pixel
variances as functions of tube potentials, fluence, and virtual monochromatic
image (VMI) energy. These expressions, along with CNR, are used as an
objective function to optimize the tube settings and minimize image noise in a
phantom experiment. We determine the optimal VMI energy for maximal CNR and
also prove that the least noisy monochromatic image synthesized from iodine
and water density images corresponds to the optimal energy VMI. Consequently,
using more general weights to linearly combine the two images does not enhance
image quality.

The problem of improving iodine CNR in VMI has been previously addressed, but
our approach differs from existing works in terms of the analysis methods
employed. Specifically, the works by Leng et al. 2 and Tao et al. 12 focus on
applying denoising techniques during the reconstruction process. Yu et al. 13,
14 examine the effect of subject size and photon fluence on the image quality
of linearly mixed images generated from DECT scans using a dual-source CT
scanner, assuming a fixed x-ray source kilovoltage-peak (kV) setting and total
radiation dose. Michalak et al. 15 conduct a phantom study to empirically
determine the optimal photon energies for virtual mono-energetic imaging
across various phantom sizes. Dabli et al. 16 empirically determine optimal
tube potential settings that yield high image quality and accuracy for low
iodine concentration quantification. Ren et al. 17 analyze the conditioning of
spectral weights by employing singular value decomposition of the matrix
formed by sampling the intensity profile of each spectral weight over the
energy range.

##  II. Methods

In this section, we detail the approaches and procedures employed to optimize
Dual-Energy Computed Tomography (DECT) scan parameters for improved iodine-
based Contrast-to-Noise Ratio (CNR). We first present the physical model for
the DECT measurements, considering the well-established assumption of the
basis material decomposition 7, 18 and a noise model based on the compound
Poisson process 19. Next, we explore the role of the Jacobian matrix in
ensuring reconstruction uniqueness and examine noise propagation by
linearizing the inverse sinogram-to-photon-counts transformation. We then
derive analytical expressions for CNR, mean sinogram variance, and mean pixel
variance. We finally show that the smallest eigenvalue of the mean pixel
covariance matrix gives the minimum mean variance for VMI obtained by using
optimal VMI energy, identifying optimal energy levels for maximum image
quality.

###  II.A. Scan configuration, DECT technique, and noise simulation

The physics modeling for the dual-energy CT system includes available models
for Tungsten X-ray source spectra, low and high kV fluence, response of
energy-integrating detectors, and compound Poisson noise 19.

For iodine-based contrast imaging with a DECT system, we assume that the
scanned object is composed of iodine and water only. Consequently, the linear
attenuation coefficient is approximated by7

$\displaystyle\mu(E,y)\approx M_{1}(E)\rho_{1}(y)+M_{2}(E)\rho_{2}(y),$ (1)

where $M_{1}(E)$ and $M_{2}(E)$ denote material attenuation at energy $E$ for
iodine and water, respectively, and $\rho_{1}(y)$ and $\rho_{2}(y)$ are their
mass density at a spatial location $y$. While the mass densities are unknown
and need to be reconstructed, the energy dependent attenuation maps are known
_a priori_ , which are available at the NIST database20. To simplify the
notation, we write

$M(E)=\begin{bmatrix}M_{1}(E)\\\ M_{2}(E)\end{bmatrix}.$

For a given x-ray beam $l$, the x-ray transform (or sinogram) of mass density
$\rho_{j}$ of the $j$-th material is given by
$x_{j}(l)=\textstyle\int_{l}\rho_{j}(y)dy$. We let

$x=x_{l}=\begin{bmatrix}x_{1}(l)\\\ x_{2}(l)\end{bmatrix}.$

Let $S_{1}(E)$ and $S_{2}(E)$ be the x-ray energy spectra corresponding to low
and high energy x-ray tube potentials $tp_{1}$ and $tp_{2}$, respectively. The
x-ray spectra used in the experiments were known _a priori_ , which were
modeled, for given tube potentials, by using the publicly available Python
software toolkit SpekPy v2.021.

For an x-ray beam $l$, the number of photons incident on the detector with
energy $E$ per unit time corresponding to the $i$-th measurement is given by

$\displaystyle I_{i}(x_{l};E)=S_{i}(E)e^{-\int_{l}\mu(E,y)dy}\approx
S_{i}(E)e^{-M(E)\cdot x_{l}},\quad 1\leq i\leq 2,$ (2)

where $M(E)\cdot x_{l}=M_{1}(E)x_{1}(l)+M_{2}(E)x_{2}(l)$.

The total number of detected photons associated to line $l$ is then given by

$\displaystyle I_{i}(x_{l})=\int_{0}^{\infty}S_{i}(E)e^{-M(E)\cdot
x_{l}}D(E)dE,\quad 1\leq i\leq 2,$ (3)

where $D(E)$ is the energy dependence of the detector sensitivity. For energy
integrating detectors, $D(E)=\alpha E$ for some $\alpha>0$. 19

In this study, we consider the negative logarithm of the intensity
measurements

$g_{i}(x)=-\ln I_{i}(x),\quad i=1,2.$ (4)

We define

$I(x)=\begin{bmatrix}I_{1}(x)\\\ I_{2}(x)\end{bmatrix},\qquad
g(x)=\begin{bmatrix}g_{1}(x)\\\ g_{2}(x)\end{bmatrix}.$

We assume that the sinograms $x=x_{l}\in\mathcal{R}=[0,a_{1}]\times[0,a_{2}]$,
a rectangle in $\mathbb{R}^{2}$ with $a_{j}$ being the maximal attenuation
over all possible lines expected for the $j$-th material in a given imaging
task.

In the reconstruction of mass density maps from DECT measurements, we consider
a two step method given in the following diagram. For each line $l$,

$\displaystyle\begin{bmatrix}g_{1}\\\
g_{2}\end{bmatrix}\quad\xlongrightarrow[\text{method}]{\text{Newton's}}\quad\begin{bmatrix}x_{1}\\\
x_{2}\end{bmatrix}\quad\xlongrightarrow[\text{back-
projection}]{\text{Filtered}}\quad\begin{bmatrix}\rho_{1}\\\
\rho_{2}\end{bmatrix}.$ (5)

For $i=1,2$, the spectral measurements $I_{i}$ are assumed to be independent
random variables that follow a compound Poisson process 19. The covariance
matrix for the log-intensity measurements $g$ is then given by

$\displaystyle\Sigma_{g}(x)=\begin{bmatrix}\sigma_{g}(x)_{1}^{2}&0\\\
0&\sigma_{g}(x)_{2}^{2}\end{bmatrix},$ (6)

where

$\displaystyle\sigma_{g}(x)_{i}^{2}=\displaystyle\frac{\int_{0}^{\infty}D^{2}(E)S_{i}(E)e^{-M(E)\cdot
x}dE}{\Big{(}\int_{0}^{\infty}D(E)S_{i}(E)e^{-M(E)\cdot x}dE\Big{)}^{2}},$ (7)

is the variance of the $i$-th log-intensity measurement $g_{i}$.

We note that $\Sigma_{g}$ depends on the detector sensitivity $D(E)$ but is
independent of the factor $\alpha$ if $D(E)$ is replaced by $\alpha D(E)$.

###  II.B. Jacobian of the sinogram-to-photon-counts transformation and
uniqueness of reconstructions

Unlike the case of single energy CT, the reconstructions obtained from DECT
measurements may not always be unique, which is mainly due to the non-
linearity of the DECT measurement model with respect to sinogram values. For
nonlinear maps defined on convex domains, local constraints on the Jacobian of
the forward measurements provide sufficient criteria for the uniqueness of
reconstructions 22, 23, 24.

In our case, the map $x\to g(x)$ is smooth and its Jacobian matrix at point
$x\in\mathcal{R}$ is given by the matrix $J(x)$ with entries

$\displaystyle J_{ij}(x)=\frac{\partial g_{i}}{\partial
x_{j}}(x)=\frac{\int_{0}^{\infty}D(E)S_{i}(E)M_{j}(E)e^{-M(E)\cdot
x}dE}{\int_{0}^{\infty}D(E)S_{i}(E)e^{-M(E)\cdot x}dE},\qquad i,j=1,2.$ (8)

Based on the work of Gale and Nikaido22, Bal and Terzioglu 10 proved that if
the Jacobian determinant

$\displaystyle\det J(x)=J_{11}(x)J_{22}(x)-J_{12}(x)J_{21}(x)\neq 0,$ (9)

for all $x\in\mathcal{R}$, the map $x\to g(x)$ is globally injective. This
means that the reconstruction of $x$ from the knowledge of $g(x)$ (or $I(x)$)
is unique.

In DECT, the values of the Jacobian determinant depend on the material basis
and the x-ray spectra 10. For iodine-water material pair, we demonstrate the
dependence of the Jacobian determinant on the x-ray spectra in Figures 2 and 3
of section III.B. .

We also present in Fig. 4 an example scan protocol where the uniqueness does
not hold. That is, the Jacobian determinant vanishes inside the rectangle and
there exist two distinct sinogram values that are mapped by $g$ to the same
measurement pair.

###  II.C. Noise propagation based on linearization of the inverse sinogram-
to-photon-counts transformation

We now present an analysis of the noise propagation from DECT measurements to
the reconstructed sinograms by considering first order Taylor approximation to
inverse sinogram-to-photon-counts transformation.

For a given sinogram value $x_{l}\in\mathcal{R}$ that corresponds to a line
$l$, we let $g^{\eta}(x_{l})$ denote the noisy DECT measurement:

$\displaystyle g^{\eta}(x_{l})=g(x_{l})+\eta,$ (10)

where $\eta$ is the noise vector. Let $x_{l}^{\eta}$ be the sinogram value
that is reconstructed from $g^{\eta}(x_{l})$. Assuming that the noise is small
and considering a first order Taylor expansion, we have

$\displaystyle x_{l}^{\eta}\approx x_{l}+J^{-1}(x_{l})\eta.$ (11)

Under the linearization regime, the $2\times 2$ covariance matrix of the
reconstructed sinograms is given by

$\displaystyle\textstyle\Sigma_{l}=\Sigma_{x_{l}}=J^{-1}(x_{l})\Sigma_{g^{\eta}}(x_{l})J^{-t}(x_{l}),$
(12)

where $-t$ denotes the transpose of inverse matrix (see eg. Cowan, p. 21)25.
Here, the diagonal entries of (12) are the variances of the iodine and water
sinograms, and the off-diagonal entries are the covariance between them (see,
for example, Roessl et al. 26 for their explicit formulas).

We define the diagonal matrices iodine and water sinogram variances and
covariance as $C^{(ij)},i,j=1,2$. Therefore,

$\displaystyle(C^{(ij)})_{ll}=(\Sigma_{l})_{ij},\quad i,j=1,2.$ (13)

We also define a matrix $Q=(q_{ij})_{i,j=1,2}$ with diagonal entries being the
mean pixel variance of each material density map whereas the off-diagonal
entries are the mean covariance between them. Let $B=(b_{pl})_{1\leq p\leq
P,\;1\leq l\leq L}$ denote the filtered back-projection matrix. Then,

$\displaystyle q_{ij}=\frac{1}{P}{\rm{tr}}(BC^{(ij)}B^{t}),\quad j=1,2.$ (14)

It is well known that the acquisition of photons for each spectral measurement
may be shortened or lengthened to minimize reconstruction errors. Let $T$ be
the total time of acquisition and $0<\tau<T$ be the acquisition time of the
low energy measurement (so that $T-\tau$ is the time of acquisition in the
high kV setting). Since the variance of the measurement noise is proportional
to photon count, we observe that the covariance matrix of the measurement
error is divided by $\tau$ when $i=1$ and $T-\tau$ with $i=2$. We thus have a
modified covariance matrix for the reconstructed sinograms

$\displaystyle\textstyle\Sigma_{l}=J^{-1}(x_{l})F\Sigma_{g^{\eta}}(x_{l})J^{-t}(x_{l}),$
(15)

where

$F=\begin{bmatrix}\tau&0\\\ 0&T-\tau\end{bmatrix}.$

It is important to note that the matrix $\Sigma_{l}$ is always symmetric, and
is also positive definite provided that $\det J(x_{l})\neq 0$. We also observe
from eqs. (12) and (15) that the sinogram (co-)variances are inversely
proportional to the square of the Jacobian determinant. Therefore, optimizing
DECT scan parameters to maximize the minimum value of the Jacobian determinant
over the sinogram domain leads to more stable sinogram reconstructions due to
the reduced amplification of the measurement noise (also see Bal et al. 11).

###  II.D. Computation of pixel variances and the CNR of a monochromatic image

Once the sinograms $x_{1}$ and $x_{2}$ (for iodine and water, respectively)
are reconstructed from the DECT measurements using a nonlinear iterative
algorithm, for which we use Newton’s method in this study, they can be
combined linearly to obtain a monochromatic sinogram, that is

$\displaystyle x_{mono}=w_{1}x_{1}+w_{2}x_{2}=w^{t}x,$ (16)

where $w=[w_{1},w_{2}]^{t}$ is a unit vector in $\mathbb{R}^{2}$, i.e,
$\|w\|_{2}=1$.

We note that for virtual monochromatic images (VMI), one considers the
attenuation weights $w=\frac{M(E)}{\|M(E)\|_{2}}$ where the VMI energy $E$ is
chosen to maximize a given metric, e.g. the contrast-to-noise ratio (CNR) of a
region of interest (ROI) in a given imaging task.

We define the contrast-to-noise ratio of a signal as the difference between
the mean CT numbers of the signal and the background divided by the mean
standard deviation of the signal:

$\displaystyle{\rm{CNR}}(w,tp,\tau)=\frac{mean(CT\\#_{signal})-mean(CT\\#_{background})}{s(w,tp,\tau)_{signal}}.$
(17)

Here, we consider as the background the part of the image containing only
water.

In the following, we derive an analytic expression for the pixel variance
$s^{2}(w,tp,\tau)$.

Let $L$ denote the total number of lines (bins $\times$ views). If the
monochromatic sinograms (16) corresponding to different lines are independent,
then their covariance, which is denoted by $\Sigma_{mono}$ and is of size
$L\times L$, is a diagonal matrix of monochromatic variances

$\displaystyle\sigma^{2}_{l}(w,tp,\tau)=w^{t}\Sigma_{l}w.$ (18)

The pixel covariance matrix is then given by

$\displaystyle\Sigma_{y}$ $\displaystyle=B\Sigma_{mono}B^{t},$ (19)

where $B=(b_{pl})_{1\leq p\leq P,\;1\leq l\leq L}$ is the filtered back-
projection matrix.

By direct calculation, we obtain that

$\displaystyle(\Sigma_{y})_{ii}=\sum_{l=1}^{L}\sigma^{2}_{l}(w,tp,\tau)b_{pl}^{2},\quad
i=1,\dots,P.$ (20)

Therefore, the mean pixel variance of the monochromatic signal is given by

$\displaystyle s^{2}(w,tp,\tau)$
$\displaystyle=\frac{1}{P}\sum_{p=1}^{P}(\Sigma_{y})_{ii}=\frac{1}{P}\sum_{p=1}^{P}\sum_{l=1}^{L}\sigma^{2}_{l}(w,tp,\tau)b_{pl}^{2}$
$\displaystyle=\frac{1}{P}\sum_{l=1}^{L}\big{(}{\textstyle\sum}_{p=1}^{P}b_{pl}^{2}\big{)}\sigma^{2}_{l}(w,tp,\tau)=\frac{1}{P}\sum_{l=1}^{L}(b_{l}^{t}b_{l})\sigma^{2}_{l}(w,tp,\tau),$
(21)

where $b_{l}$ denote the $l$-th column of $B$. Now using eq. (18), we obtain
that

$\displaystyle s^{2}(w,tp,\tau)$
$\displaystyle=\frac{1}{P}\sum_{l=1}^{L}(b_{l}^{t}b_{l})w^{t}\Sigma_{l}w=w^{t}\left(\frac{1}{P}\sum_{l=1}^{L}b_{l}^{t}b_{l}\Sigma_{l}\right)w=w^{t}Qw,$
(22)

since

$\displaystyle
q_{ij}=\frac{1}{P}{\rm{tr}(BC^{(ij)}B^{t})}=\frac{1}{P}\sum_{l=1}^{L}\sum_{p=1}^{P}b_{pl}^{2}(\Sigma_{l})_{ij}=\frac{1}{P}\sum_{l=1}^{L}b_{l}^{t}b_{l}(\Sigma_{l})_{ij},\quad
j=1,2.$ (23)

###  II.E. Minimization of mean pixel variance

In the following, we show that the smallest eigenvalue of the mean pixel
covariance matrix gives the minimum mean variance for VMI obtained by using
optimal VMI energy, identifying optimal energy levels for maximum image
quality.

We first observe that if $\det J(x_{l})\neq 0$ for all $l$, then $\Sigma_{l}$
is positive definite for all $l$, and thus the matrix
$Q=\textstyle\frac{1}{P}\sum_{l=1}^{L}b_{l}^{t}b_{l}\Sigma_{l}$ is positive
definite. This implies that $s^{2}(w,tp,\tau)=w^{t}Qw$ is a positive definite
quadratic form.

When the matrix $Q$ is positive definite, by the spectral theorem, it is
orthogonally diagonalizable, that is the eigenvectors of $Q$ form an
orthogonal basis for $\mathbb{R}^{2}$. For $i=1,2$, we let
$Qu_{i}=\lambda_{i}u_{i}$ where $\|u_{i}\|_{2}=1$, and
$\lambda_{1}\geq\lambda_{2}>0$, i.e., $u_{i}$ is the unit eigenvector
associated to eigenvalue $\lambda_{i}$ (see eg. Horn and Johnson 27). Then,
for each tube potential pair $tp$ and fluence $\tau$,

$\displaystyle\min_{\|w\|_{2}=1}s^{2}(w,tp,\tau)=\min_{\|w\|_{2}=1}w^{t}Q(tp,\tau)w=\lambda_{2}(tp,\tau),$
(24)

where the minimizer is $u_{2}(tp,\tau)$. One can then minimize $\lambda_{2}$
over $tp$ and $\tau$ to find the optimal tube potential pair and the photon
fluence. Hence,

$\displaystyle(w^{*},tp^{*},\tau^{*})=\operatorname*{arg\,min}_{w,tp,\tau}s^{2}(w,tp,\tau),$
(25)

where

$\displaystyle(tp^{*},\tau^{*})=\operatorname*{arg\,min}_{tp,\tau}\lambda_{2}(tp,\tau),$
(26)

and $w^{*}$ is the unit eigenvector of $Q(tp^{*},\tau^{*})$ corresponding to
$\lambda_{2}(tp^{*},\tau^{*})$.

We observe that if the Jacobian determinant values increase, then the diagonal
entries of ($\Sigma_{l}$ and hence) $Q$ decrease. Since the smallest
eigenvalue of a $2\times 2$ matrix is always less than the smallest diagonal
entry (see eg. Horn and Johnson 27), scan parameters that maximize Jacobian
determinant values give the minimal smallest eigenvalue of $Q$ and hence it
reduces the mean pixel variance.

In general, when the mean pixel variance is minimized only over the set of the
attenuation weights, one expects to obtain a larger value in comparison to
generalized weights. However, for the iodine-water material pair, we
numerically observed that there is $E^{*}$ such that
$\frac{M(E^{*})}{\|M(E^{*})\|_{2}}=w^{*}$. This is mainly because
$M_{1}(E)>M_{2}(E)$ for all $E$ in the diagnostic energy range and the
components of $w^{*}$ have the same sign. We thus have

$\displaystyle\min_{0<E<E_{max}}s^{2}(E,tp,\tau)=\min_{\|w\|_{2}=1}s^{2}(w,tp,\tau)=\lambda_{2}(tp,\tau).$
(27)

We remark that this result is specific to iodine-water material pair and may
not hold for others.

##  III. Results

In this section, we present the results of our numerical experiments conducted
for the iodine-water material pair. We first explain our experimental setup.
We then demonstrate the dependence of uniqueness to the Jacobian determinant
values and showcase a non-unique reconstruction. We finally provide our
numerical results on the optimization of the VMI energy and iodine CNR in
virtual monochromatic images.

###  III.A. Experimental setup

In the numerical experiments, we used a $4\times 4\;cm^{2}$ phantom consisting
of a large water disk with inserted iodine-solution and calcium disk signals
(see Fig. 1). The centers and radii of each disk, and the concentration of
iodine solutions are given in table 1.

Figure 1: The plot of the circular water phantom with inserted iodine-solution and calcium disk signals. The red circle indicates the disk with the lowest iodine concentration, which is our region of interest. Material | Center (cm) | Radius (cm)
---|---|---
0.1% Iodine | (-1,1) | 0.3
0.2% Iodine | (1,1) | 0.3
0.5% Iodine | (1,-1) | 0.3
1% Iodine | (-1,1) | 0.3
Calcium | (-2,0) | 0.1
Calcium | (0,-2) | 0.2
Calcium | (2,0) | 0.3
Calcium | (0,2) | 0.4
Water | (0,0) | 3.6
Table 1: Phantom configuration.

For the phantom shown in Fig. 1, the DECT measurements were numerically
simulated for $512\times 512$ lines (bins $\times$ views) in fan beam
geometry. The distance from source to detector was 100 cm. For each line, 100
realizations were obtained by considering a compound Poisson process19. The
range of x-ray tube voltages was considered to be 30-150 kV. The x-ray beam
filter material was fixed at 2mm of Aluminum. The results were not
significantly affected by variations in the total number of photons and
fluence between low and high kV scans. As a result, these factors were
considered equal in the analysis.

###  III.B. Uniqueness of the reconstructed sinograms

As mentioned in section II.B., the uniqueness of reconstructions is ensured if
the Jacobian determinant is nonzero everywhere in the rectangle containing all
possible sinogram values (or pathlengths). Jacobian determinant values vary
according to the chosen tube potentials and the filters of the x-ray energy
spectra, where we fix the latter and focus on examining the effect of the
former.

According to the eqs. (12) and (15), the sinogram (co-)variances are inversely
proportional to the square of the Jacobian determinant. In Figure 2, we
present the plot of the quantity

$\displaystyle\min_{x\in\mathbb{R}}|\det J(x)|,$ (28)

This is the minimum of absolute values of Jacobian determinant over the
rectangular region $\mathcal{R}=[0,0.01]\times[0,7.2]$, which corresponds to
all possible pathlengths for iodine and water obtained from our phantom
(depicted in Figure 1). We caution that the Jacobian determinant also depends
on the tube settings, although we omit the specific notation for simplicity.

This plot reveals that lighter regions, which correspond to higher absolute
values of the Jacobian determinant, indicate tube potentials that lead to more
stable transmission-to-sinogram transformation. Therefore, the use of the tube
potentials in the region $50\leq tp_{1}\leq 80$ and $120\leq tp_{2}\leq 150$
yield more stable sinogram reconstructions.

Figure 2: The minimum of absolute value of Jacobian determinants over
pathlengths up to 0.01 and 7.2, for iodine and water respectively, as a
function of tube settings varying from 30 to 150 kV.

In Fig. 3, for a fixed high tube setting of $tp_{2}=120$ kV, we plot in black
the extremal Jacobian determinant values attained in the rectangle
$\mathcal{R}=[0,0.01]\times[0,7.2]$ cm2 for low tube setting varying from 30
to 90 kV. The red curve represents the Jacobian determinant of the linearized
map, which also corresponds to the Jacobian determinant values at zero
pathlength.

Figure 3: Plot of the Jacobian determinant for the map that transforms
iodine/water sinograms to dual-energy transmission at low and high kV. The
black curves show the range of determinant values over the expected tissue
path lengths and the red curve shows the determinant for the linearized
function. For the shown results, the high kV setting is 120, and the low kV
setting is on the x-axis. Maximizing determinant values leads to stable
transmission-to-sinogram transformation.

We observe in Fig. 3 that if the low tube potential is less than 45 kV, the
Jacobian determinant vanishes somewhere inside the rectangle, and hence the
uniqueness of sinogram reconstructions is not guaranteed. In fact, for the
tube potential pair $(tp_{1},tp_{2})=(35,120)$, for which the Jacobian
determinant becomes zero near the sinogram values $x=(x_{1},x_{2})=(4,6)$, we
encountered singular solutions to the inverse sinograms-to-DECT measurements
map. In Fig. 4, we present a plot of the level curves of the corresponding
log-intensity measurements as a function of the iodine and water sinogram
values. The red dot is where the Jacobian determinant is zero. The level curve
of $g_{2}=2.4$ (shown in dashed black) intersects the level curve of
$g_{1}=4.92$ (shown in blue) at two points (shown with blue dots). This means
that two different sinogram value pairs (shown with blue dots) lead to the
same measurement value of $(g_{1},g_{2})=(4.92,2.4)$.

Figure 4: Plot of the level curves of the log-intensity measurements as a
function of the water and iodine sinogram values demonstrating the existence
of singular solutions to the inverse sinograms-to-DECT measurements map. The
red dot is a point where the Jacobian vanishes. The dashed black curve
($g_{2}=2.4$) intersects the blue one ($g_{1}=4.92$) at two twice (shown with
blue dots), so both sinogram value pairs lead to the same log-intensity
measurement value of $(g_{1},g_{2})=(4.92,2.4)$.

###  III.C. Optimization of the iodine CNR in virtual monochromatic images

In this section, we present the results of the analysis conducted on the
Contrast-to-Noise Ratio (CNR) values of the 0.1% Iodine disk in the virtual
monochromatic image (VMI) which is enclosed by the red circle in Fig. 1. To
compare the analytical and empirical approaches, the CNR values were computed
using equation (17) analytically and through image reconstruction for 100
realizations. In the calculations, the 0.1% Iodine disk and the water only
parts of the phantom were considered as the signal and the background,
respectively.

Fig. 5 displays these CNR values for the low and high tube settings of 60/120
kV, with the blue curve representing the analytically computed values and the
red dots indicating the empirically obtained results. Notably, there is a
remarkable agreement between the empirical and analytical findings.
Furthermore, the maximum CNR for the 0.1% Iodine disk in the Virtual
Monochromatic Images (VMI) was observed at an energy level of 55 keV. Similar
results were obtained for the 60/140 kV setting. These findings provide
valuable insights into the optimal energy range for maximizing CNR in VMI.

Figure 5: For low/high tube settings of 60/120 kV, the CNR of the 0.1% Iodine
disk is shown for a noise level of 50,000 photons incident on each detector
pixel. Shown are the analytically computed CNR values with blue curve and the
empirically calculated values using noise realizations, with the red dots. The
agreement between theory and simulation validates the computational method. A
clear peak in the 0.1% Iodine disk CNR is seen in the VMI for an energy of 55
keV.

To further explore the influence of tube potentials on the CNR, we plot in
Fig. 6 the CNR values for various tube settings at this optimal VMI energy of
55 keV. Comparing the high tube potential settings of 120 kV and 140 kV, it is
observed that the former yields slightly higher CNR values. Consequently, the
tube potential pair of 60/120 kV emerges as the optimal choice for maximizing
the CNR of the 0.1% Iodine disk in VMI within the scope of this study. We
remark that, in view of Fig. 2, this pair of tube potentials lead to stable
transmission-to-sinogram transformation. This result highlights the
significance of selecting the appropriate kV settings to optimize the
visualization and diagnostic quality of the Iodine disk in VMI imaging.

Figure 6: Using the VMI energy that maximizes the 0.1% Iodine disk CNR, the
low and high kV settings are varied. Shown in the plot are only the high kV
settings of 120 and 140; the low kV setting is indicated on the x-axis. Of the
values computed, the 60/120 kV setting allows for maximum 0.1% Iodine CNR.

##  IV. Discussion and conclusion

In this study, we have introduced a novel approach for optimizing the DECT
scan protocol specifically for iodine-based CNR. Our methodology centers
around analyzing the propagation of noise from DECT measurements to material
density images. By linearizing the inverse sinogram-to-transmission
measurements map, we were able to derive analytic expressions for the mean
sinogram and pixel variances. These expressions serve as key components of our
framework, enabling us to systematically optimize the DECT scan parameters for
improved image quality in terms of iodine-based CNR. We also identified the
ideal VMI energy that maximizes CNR. We have shown that the VMI synthesized
from iodine and water density images with the optimal energy exhibits the
least amount of noise. As a result, our findings indicate that employing more
general weights to linearly combine these images does not provide any
substantial improvement in image quality.

Our DECT simulations were based on dual source DECT scanners. However, our
theoretical framework is applicable to dual-layer detectors and photon
counting detectors as well. In those cases, one needs to optimize the energy
window thresholds instead of the tube potentials. The photon fluence was taken
into account in our framework. However, we considered equal times for both
energy measurement in the numerical simulations as its influence on the
results was not significant. We remark that the results of this study are
specific to the phantom used as it directly affects the size of the sinogram
domain and consequently the range of Jacobian values observed. Our primary
objective in this research was to propose a framework for optimizing DECT scan
parameters to achieve maximum stability of reconstructions and contrast-to-
noise ratio (CNR). It is important to consider the specific imaging task at
hand and carefully assess all relevant factors when implementing this
framework.

In a future study, we plan to test our theoretical results using real DECT
data. We also plan to extend our analysis to the case of three or more
materials and energy measurements.

## Acknowledgment

F. Terzioglu’s work was supported in part by the NSF DMS grant 2206279. G.
Bal’s work was supported in part by the NSF DMS grants 2306411 and 1908736.
This work is also supported in part by NIH Grant Nos. R01-EB023968 and
R21-CA263660. The contents of this article are solely the responsibility of
the authors and do not necessarily represent the official views of the
National Institutes of Health.

##

  * 1 K. L. Grant, T. G. Flohr, B. Krauss, M. Sedlmair, C. Thomas, and B. Schmidt,  Assessment of an advanced image-based technique to calculate virtual monoenergetic computed tomographic images from a dual-energy examination to improve contrast-to-noise ratio in examinations using iodinated contrast media,  Investigative radiology 49, 586–592 (2014).
  * 2 S. Leng, L. Yu, J. G. Fletcher, and C. H. McCollough,  Maximizing iodine contrast-to-noise ratios in abdominal CT imaging through use of energy domain noise reduction and virtual monoenergetic dual-energy CT,  Radiology 276, 562 (2015).
  * 3 K. Kalisz, N. Rassouli, A. Dhanantwari, D. Jordan, and P. Rajiah,  Noise characteristics of virtual monoenergetic images from a novel detector-based spectral CT scanner,  European Journal of Radiology 98, 118–125 (2018).
  * 4 D. Sakabe, Y. Funama, K. Taguchi, T. Nakaura, D. Utsunomiya, S. Oda, M. Kidoh, Y. Nagayama, and Y. Yamashita,  Image quality characteristics for virtual monoenergetic images using dual-layer spectral detector CT: comparison with conventional tube-voltage images,  Physica Medica 49, 5–10 (2018).
  * 5 A. Tabari, M. S. Gee, R. Singh, R. Lim, K. Nimkin, A. Primak, B. Schmidt, and M. K. Kalra,  Reducing radiation dose and contrast medium volume with application of dual-energy CT in children and young adults,  American Journal of Roentgenology 214, 1199–1205 (2020).
  * 6 G. N. Hounsfield,  Computerized transverse axial scanning (tomography): Part 1. Description of system,  The British Journal of Radiology 46, 1016–1022 (1973).
  * 7 R. E. Alvarez and A. Macovski,  Energy-selective reconstructions in x-ray computerised tomography,  Physics in Medicine & Biology 21, 733 (1976).
  * 8 Z. H. Levine,  Nonuniqueness in dual-energy CT,  Medical physics 44, e202–e206 (2017).
  * 9 R. E. Alvarez,  Invertibility of the dual energy x-ray data transform,  Medical Physics 46, 93–103 (2019).
  * 10 G. Bal and F. Terzioglu,  Uniqueness criteria in multi-energy CT,  Inverse Problems 36, 065006 (2020).
  * 11 G. Bal, R. Gong, and F. Terzioglu,  An inversion algorithm for P-functions with applications to multi-energy CT,  Inverse Problems 38, 035011 (2022).
  * 12 S. Tao, K. Rajendran, W. Zhou, J. G. Fletcher, C. H. McCollough, and S. Leng,  Improving iodine contrast to noise ratio using virtual monoenergetic imaging and prior-knowledge-aware iterative denoising (mono-PKAID),  Physics in Medicine & Biology 64, 105014 (2019).
  * 13 L. Yu, A. N. Primak, X. Liu, and C. H. McCollough,  Image quality optimization and evaluation of linearly mixed images in dual-source, dual-energy CT,  Medical physics 36, 1019–1024 (2009).
  * 14 L. Yu, J. A. Christner, S. Leng, J. Wang, J. G. Fletcher, and C. H. McCollough,  Virtual monochromatic imaging in dual-source dual-energy CT: radiation dose and image quality,  Medical physics 38, 6371–6379 (2011).
  * 15 G. Michalak, J. Grimes, J. Fletcher, A. Halaweish, L. Yu, S. Leng, and C. McCollough,  Selection of optimal tube potential settings for dual-energy CT virtual mono-energetic imaging of iodine in the abdomen,  Abdominal Radiology 42, 2289–2296 (2017).
  * 16 D. Dabli, J. Frandon, A. Hamard, A. Belaouni, T. Addala, J.-P. Beregi, and J. Greffier,  Optimization of image quality and accuracy of low iodine concentration quantification as function of kVp pairs for abdominal imaging using dual-source CT: A phantom study,  Physica Medica 88, 285–292 (2021).
  * 17 Y. Ren, H. Xie, W. Long, X. Yang, and X. Tang,  On the Conditioning of Spectral Channelization (Energy Binning) and Its Impact on Multi-Material Decomposition Based Spectral Imaging in Photon-Counting CT,  IEEE Transactions on Biomedical Engineering 68, 2678–2688 (2021).
  * 18 J. F. Williamson, S. Li, S. Devic, B. R. Whiting, and F. A. Lerma,  On two-parameter models of photon cross sections: application to dual-energy CT imaging,  Medical physics 33, 4115–4129 (2006).
  * 19 B. R. Whiting,  Signal statistics in x-ray computed tomography,  in Medical Imaging 2002: Physics of Medical Imaging, volume 4682, pages 53–60, International Society for Optics and Photonics, 2002.
  * 20 J. H. Hubbell and S. M. Seltzer,  Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients 1 keV to 20 MeV for elements Z= 1 to 92 and 48 additional substances of dosimetric interest,  Technical report, National Inst. of Standards and Technology-PL, Gaithersburg, MD (United States), 1995.
  * 21 G. Poludniowski, A. Omar, R. Bujila, and P. Andreo,  SpekPy v2. 0—a software toolkit for modeling x-ray tube spectra,  Medical Physics 48, 3630–3637 (2021).
  * 22 D. Gale and H. Nikaido,  The Jacobian matrix and global univalence of mappings,  Mathematische Annalen 159, 81–93 (1965).
  * 23 C. Garcia and W. Zangwill,  On univalence and P-matrices, Report 7737,  Center for Mathematical Studies in Business and Economics, University of Chicago (1977).
  * 24 A. Mas-Colell,  Homeomorphisms of Compact, Convex Sets and the Jacobian Matrix,  SIAM Journal on Mathematical Analysis 10, 1105–1109 (1979).
  * 25 G. Cowan,  Statistical Data Analysis,  Oxford University Press, 1998.
  * 26 E. Roessl, A. Ziegler, and R. Proksa,  On the influence of noise correlations in measurement data on basis image noise in dual-energylike x-ray imaging,  Medical physics 34, 959–966 (2007).
  * 27 R. A. Horn and C. R. Johnson,  Matrix Analysis,  Cambridge University Press, 2012.


# FeatureBox: Feature Engineering on GPUs
for Massive-Scale Ads Systems

Weijie Zhao, Xuewu Jiao, Xinsheng Luo, Jingxue Li, Belhal Karimi, Ping Li
Cognitive Computing Lab, Baidu Research
Baidu Search Ads (Phoenix Nest), Baidu Inc.
10900 NE 8th St. Bellevue, Washington 98004, USA
No. 10 Xibeiwang East Road, Beijing 100193, China

{weijiezhao, jiaoxuewu, luoxinsheng, lijingxue01, belhalkarimi,
<EMAIL_ADDRESS>
###### Abstract

Deep learning has been widely deployed for online ads systems to predict
Click-Through Rate (CTR). Machine learning researchers and practitioners
frequently retrain CTR models to test their new extracted features. However,
the CTR model training often relies on a large number of raw input data logs.
Hence, the feature extraction can take a significant proportion of the
training time for an industrial-level CTR model. In this paper, we propose
FeatureBox, a novel end-to-end training framework that pipelines the feature
extraction and the training on GPU servers to save the intermediate I/O of the
feature extraction. We rewrite computation-intensive feature extraction
operators as GPU operators and leave the memory-intensive operator on CPUs. We
introduce a layer-wise operator scheduling algorithm to schedule these
heterogeneous operators. We present a light-weight GPU memory management
algorithm that supports dynamic GPU memory allocation with minimal overhead.
We experimentally evaluate FeatureBox and compare it with the previous in-
production feature extraction framework on two real-world ads applications.
The results confirm the effectiveness of our proposed method.

###### Index Terms:

CTR Prediction; GPU; Large-Scale Machine Learning Framework

##  I Introduction

Deep learning has been widely employed in many real-world applications, e.g.,
computer vision [10, 13, 28, 5], data mining [8, 22, 14, 17, 21, 33], and
recommendation systems [3, 2, 29, 35, 20, 15, 32]. In recent years, sponsored
online advertising also adopts deep learning techniques to predict the Click-
Through Rate (CTR) [7, 38, 19, 40, 9, 24, 27, 36, 34, 31]. Unlike common
machine learning applications, the accuracy of the CTR prediction is critical
to the revenue. In the context of a many-billion-dollar online ads industry,
even a $0.1\%$ accuracy increase will result in a noticeable revenue gain
[37]. In this work, we identify two major paths to improve the model accuracy.
The first area is to propose different and enhanced model architectures. Every
improvement in this direction is considered a fundamental milestone in the
deep learning community—and does not happen often in the CTR prediction
industry. The other (more practical) is feature engineering, i.e., to propose
and extract new features from the raw training data. The benefit of feature
engineering is usually neglected in common deep learning applications because
of the general belief that deep neural networks inherently extract the
features through their hidden layers. However, recall that CTR prediction
applications are accuracy-critical, hence, the gain from an improved feature
engineering strategy remains attractive for in-production CTR prediction
models. Therefore, in order to achieve a better prediction performance, CTR
deep learning models in real-world ads applications tend to utilize larger
models and more features extracted from raw data logs.

Testing on the historical and online data is the rule-of-the-thumb way to
determine whether a new feature is beneficial. Every new feature with positive
accuracy improvement (e.g., $0.1\%$) is included into the CTR model. Machine
learning researchers and practitioners keep this feature engineering trial-
and-error on top of the current in-production CTR model. As a result, the in-
production CTR model becomes larger and larger with more and more features. To
support the trial-and-error research for new features, it requires us to
efficiently train massive-scale models with massive-scale raw training data in
a timely manner. Previous studies [37] propose hierarchical GPU parameter
server that trains the out-of-memory model with GPU servers to accelerate the
training with GPUs and SSDs. With a small number of GPU servers, e.g., 4, can
obtain the same training efficiency as a CPU-only cluster with hundreds of
nodes. The training framework focuses on the training stage and assumes the
training data are well-prepared—the training data are accessed from a
distributed file system.

However, preparing the training data is not trivial for industrial level CTR
prediction models—with $\sim$ $10^{12}$ features. The feature extraction from
raw data logs can take a significant proportion of the training time. In
addition to the frequent retraining for new feature engineering trials, online
ads systems have to digest a colossal amount of newly incoming data to keep
the model up-to-date with the optimal performance. For the rapid training
demands, optimizing the feature extraction stage becomes one of the most
desirable goals of online ads systems. This latter point is the scope of our
contribution.

Figure 1: A visual illustration for the original feature extraction and
training workflow (upper); and our proposed FeatureBox (lower).

Training workflow. The upper part of Figure 1 depicts a visual illustration of
the feature extraction. Due to the large amount of raw data, the original
feature extraction task is constructed as MapReduce [4] jobs that compute
feature combinations, extract keywords with language models, etc. Those
MapReduce jobs frequently read and write intermediate files with the
distributed file system (i.e., HDFS [1]). The intermediate I/O can be as large
as 200 TB. Once the features are extracted, we also need to materialize them
to the $\sim$$15$ TB extracted features to the HDFS so that the following
distributed training framework can read them from the distributed file system.
This training workflow incurs rapid communication with HDFS that generates
heavy I/O overhead.

One straightforward question can be raised: Can we perform the feature
extraction within GPU servers to eliminate the communication overhead? In the
lower part of Figure 1, we depict an example for the proposed training
framework that combines the feature extraction and the training computation
within GPU servers. The intermediate I/O is eliminated by integrating the
feature extraction and the training computation into a pipeline: for each
batch of extracted features, we feed the batch to the model training without
writing them as intermediate files into HDFS.

Challenges & Approaches. However, moving the feature extraction to GPU servers
is non-trivial. Note that the number of GPU nodes is much fewer compared with
the CPU-only cluster. We acknowledge two main challenges in embedding the
feature extraction phase into GPU servers:

  1. 1.

Network I/O bandwidth. The network I/O bandwidth of GPU servers is by orders
of magnitude smaller than the bandwidth of CPU clusters because we have fewer
nodes—the total number of network adapters is lower. We materialize
frequently-used features as basic features so that we can reuse them without
extra I/O and computations. In addition, we use column-store that reads only
the required columns in the logs to reduce I/O.

  2. 2.

Computing Resources. With a smaller number of nodes, the CPU computing
capability on GPU servers is also orders of magnitudes less powerful than the
CPU cluster. We have to move the CPU computations to GPU operations to bridge
the computing power gap.

  3. 3.

Memory Usage. The feature extraction process contains many memory-intensive
operations, such as dictionary table lookup, sort, reduce, etc. It is desired
to have an efficient memory management system to efficiently perform dynamic
memory allocations on GPU servers with limited memory.

We summarize our contributions as follows:

  * •

We propose FeatureBox, a novel end-to-end training framework that pipelines
the feature extraction and the training on GPU servers.

  * •

We present a layer-wise operator scheduling algorithm that arranges the
operators to CPUs and GPU.

  * •

We introduce a light-weight GPU memory management algorithm that supports
dynamic GPU memory allocation with minimal overhead.

  * •

We experimentally evaluate FeatureBox and compare it with the previous in-
production feature extraction framework on two real-world ads applications.
The results confirm the effectiveness of our proposed methods.

##  II Preliminary

In this section, we present a brief introduction of CTR prediction models and
the hierarchical GPU parameter server. Both concepts are the foundations of
FeatureBox.

###  II-A CTR Prediction Models

About a decade ago, CTR prediction strategies with large-scale logistic
regression model on carefully engineered features are proposed in [6, 11].
With the rapid development of deep learning, deep neural networks (DNN)
attract a lot of attention in the CTR research community: The DNN model, with
wide embedding layers, obtains significant improvements over classical models.
The model takes a sparse high-dimensional vector as input and converts those
sparse features into dense vectors through sequential embedding layers. The
output dense vector is considered a low-dimensional representation of the
input and is then fed into the following layers in order to compute the CTR.
Most proposed CTR models share the same embedding layer architecture and only
focus on the following neural network layers, see for e.g., Deep Crossing
[26], Product-based Neural Network (PNN) [25], Wide&Deep Learning [2], YouTube
Recommendation CTR model [3], DeepFM [12], xDeepFM [18] and Deep Interest
Network (DIN) [39]. They introduce special neural layers for specific
applications that capture latent feature interactions.

Figure 2: An example for the CTR prediction network architecture.

We summarize those architectures in Figure 2. The input features are fed to
the neural network as a sparse high-dimensional vector. The dimension of the
vector can be $\sim$$10^{12}$ or more. The input features for CTR models are
usually from various resources with categorical values, e.g., query words, ad
keywords, and user portrait. The categorical values are commonly represented
as a one-hot or multi-hot encoding. Therefore, with categorical values with
many sources, the number of dimensions is high ($\sim$$10^{12}$) for industry
CTR prediction models. Note that, as demonstrated in [37], feature compression
or hashing strategies [30, 16] that reduce the number of dimensions are not
fully applicable to the CTR prediction model because those solutions
inevitably trade off the prediction accuracy for better computational
time—recall that even a small accuracy loss leads to a noticeable online
advertising revenue decrease, which is unacceptable. We embed the high-
dimensional features through an embedding layer to obtain a low-dimensional
($\sim$$10^{3}$) representation. The number of parameters in the embedding
layer can be 10 TB or more due to the high input dimension. After the low-
dimensional embedding is obtained, we fed this dense vector to the neural
network components to compute the CTR.

###  II-B Hierarchical GPU Parameter Server

Due to the extremely high dimension of the embedding layer, the model contains
more than 10 TB parameters which do not fit on most computing servers.
Conventionally, the huge model is trained on an MPI cluster. We partition the
model parameters across multiple computing nodes (e.g., 150 nodes) in the MPI
cluster. Every computing node is assigned a batch of training data streamed
directly from the HDFS. For each node, it retrieves the required parameters
from other nodes and computes the gradients for its current working mini-
batch. The gradients are then updated to the nodes that maintain the
corresponding parameters through MPI communications. Recently, hierarchical
GPU parameter servers [37] are proposed to train the massive-scale model on a
limited number of GPU servers. The key observation of the hierarchical GPU
parameter server is that the number of referenced parameters in a mini-batch
fits the GPU memory because the input vector is sparse. It maintains three
levels of hierarchical parameter servers on GPU, CPU main memory, and SSD. The
working parameters are stored in GPUs, the frequently used parameters are kept
in CPU main memory, and other parameters are materialized as files on SSDs.
The upper-level module acts as a high-speed cache of the lower-level module.
With 4 GPU nodes, the hierarchical GPU parameter server is able to be 2X
faster than 150 CPU-only nodes in an MPI cluster. Our proposed FeatureBox
follows the design of the training framework in the hierarchical GPU parameter
server and absorbs the feature engineering workload into GPUs to eliminate
excessive intermediate I/O.

##  III FeatureBox Overview

In this section, we present an overview of FeatureBox. We aim at allowing the
training framework to support pipeline processing with mini-batches so that we
can eliminate the excessive intermediate resulting I/O in conventional stage-
after-stage methods. Figure 3 depicts the detailed workflow of the FeatureBox
pipeline.

Figure 3: FeatureBox pipeline. Figure 4: An example for the heterogeneous
operator scheduling.

The workflow in Figure 3 has two major tracks: – extract features from input
views and – reading basic features. A view is a collection of raw data logs
from one source, e.g., user purchase history. CTR prediction models collect
features from multiple sources to obtain the best performance. The views are
read from the network file system HDFS. We need to clean the views by filling
null values and filtering out unrelated instances. Afterwards, the views are
joined with particular keys such as user id, ads id, etc. We extract features
from the joined views to obtain the desired features from the input views.
Then, these features are merged with the basic features, read in a parallel
path. We provide a detailed illustration for these operations as follows:

Read views and basic features. The views and basic features are streamed from
the distributed file system. The features are organized in a column-wise
manner so that we only need to read the required features.

Clean views. Views contain null values and semi-structured data, e.g., JSON
format [23]. At the view cleaning stage, we fill the null values and extract
required fields from the semi-structured data. Following the cleaning, all
columns have non-empty and simple type (as integer, float, or string) fields.
Note that the resulting views contain all the logged instances. For an
application, it may not need to include all instances, e.g., an application
for young people. A custom filter can be applied to filter out unrelated
instances of the current application.

Join views. We now have one structured table for each view. Data from
different views are concatenated by joining their keys, e.g., user id, ad id,
etc. We recall that the join step combines multiple views into a single
structured table.

Extract features. Every time CTR model engineers propose a new feature, an
operator that computes the new feature extraction on the structured table is
created. A collection of those operators are executed in the feature
extraction stage. The FeatureBox framework figures out the dependencies of
operators and schedules the execution of the operators.

Merge features. The extracted features are further merged with the basic
features read from HDFS. The merging is also realized by a join operation on
the instance id, which is a unique value generated when an instance is logged.
Subsequent to the merging, a mini-batch of training data is generated and is
fed to the neural network for the training.

##  IV Heterogeneous Operator Scheduling

The stages discussed above are represented as operators in the FeatureBox
pipeline. Note that those operators are heterogeneous: Some operators are
network I/O intensive, e.g., read views and read basic features; some
operators are computation-intensive, e.g., clean views and extract features;
and the remaining operators with joining, e.g., join views and merge features,
rely on heavy memory consumption for large table joins (which corresponds to a
large dictionary lookup). Therefore, we introduce a heterogeneous operator
scheduler that manages the operator execution on both CPUs and GPUs.

Scheduling. Figure 4 shows an example for the heterogeneous operator
scheduling algorithm. We first present a function call graph for operators in
Figure 4(a). Three operators and three major functions are displayed in the
example. Op1 calls Func3; Op2 calls Func1 and Func3; and Op3 calls Func2 and
Func3, where Func1 and Func2 are pre-processing calls, and Func3 is a post-
processing call. We make a fine granularity pipeline so that the initialing
overhead of the pipeline is minimized. The fine-granularity is obtained by
viewing each function call as a separate operator. Then, we obtain 5 more
operators: Op4 is a call for Func1; Op5 is a call for Func2; Op6, Op7, and Op8
are the Func3 calls from Op1, Op2, and Op3, respectively. Their dependency
graph is illustrated in Figure 4(b). Now we have a directed acyclic graph
(DAG) for the operators. As shown in Figure 4(c), we perform a topological
sort on the dependency graph, assign the operators with no dependencies (root
operators) to the first layer, and put the remaining operators to the
corresponding layer according to their depth from the root operators. With
this layer-wise partition, we observe that the operators in the same layer do
not have any execution dependency. We issue the operators in the same layer
together and perform a synchronization at the end of each layer to ensure the
execution dependency. We prefer to execute operators on GPUs unless an
operator requires a significant memory footprint that does not fit in the GPU
memory. For instance, Op5 (Func2) in Figure 4 is a word embedding table look
up operation that requires a considerable amount of memory. We assign this
operation to CPU workers and move its results from the CPU main memory to GPUs
as a host-to-device (H2D) CUDA call.

Inner-GPU operator launching. After the layer-wise DAG operator scheduling, we
have determined the execution device for each operator and the synchronization
barriers. However, CUDA kernel launching is has a noticeable overhead. We
report the CUDA kernel launch overhead in Table I.

TABLE I: The kernel launching overhead with an empty kernel on Nvidia Tesla V100-SXM2-32GB. #Launches | 1 | 10 | 100 | 1,000 | 10,000
---|---|---|---|---|---
Time (us) | 4 | 35 | 360 | 3,619 | 34,515

The test is performed on an Nvidia Tesla V100-SXM2-32GB GPU for an empty
kernel with 5 pointer-type arguments. The CUDA driver version is 10.2. The
average launching time for a kernel is around 3.5 us. Since we have fine-
granularity operators, we have to rapidly launch CUDA kernels to execute the
large number of operators. In order to eliminate the launching overhead, we
rewrite the operator kernel as a CUDA device function for each operator in the
same layer and create a meta-kernel that sequentially executes the operator
device functions in a runtime-compilation manner. The overhead of the meta-
kernel generation is disregarded—we only need to create this meta-kernel for
each layer once as a pre-processing of the training since we determine the
operator execution order before the actual training phase and keep the
scheduling fixed. With the generated meta-kernels, we only need to launch one
kernel for each layer.

##  V GPU Memory Management

Feature extraction operators usually need to cope with strings of varying
length, e.g., query keywords and ads titles. The execution of the operator
commonly dynamically allocates memory to process the strings. For example,
splitting a string with a delimiter needs to allocate an array to store the
result of the splitting operation. We propose a light-weight block-level GPU
memory pool to accelerate this dynamic allocation.

Figure 5: A visual illustration for the GPU memory pool architecture.

Figure 5 presents a visual illustration for our proposed block-level GPU
memory pool. The Thread Offsets denotes an array that stores the pointers to
the dynamically allocated memory in the GPU memory pool. The memory in the GPU
memory pool is pre-allocated in the GPU global memory. For each block, the
allocated memory is aligned in 128 bytes for a cache-friendly execution.

Dynamic GPU memory allocation. Algorithm 1 describes the workflow of the in-
kernel dynamic memory allocation. We maintain a global variable
idle_memory_head that stores the pointer of the head address of our pre-
allocated GPU memory pool. We assume each GPU thread in a block has computed
their required allocation size $\textit{size}_{i}$. We first compute an in-
block parallel prefix sum on $\textit{size}_{1..N}$ to obtain the prefix sum
$\textit{prefix}_{1..N}$, where $N$ is the number of threads in a block. The
prefix sum is used to compute the total size of the requested memory. In
addition, we can easily compute the thread offsets by adding the prefix sum to
the head of the allocated memory address. After that, we let one thread in the
block, e.g., thread 1, to apply the memory for the entire block—the total size
is $\textit{prefix}_{N}$. The memory allocation is implemented by an
atomic_add operation. Line 2 calls the CUDA atomic add that adds
$\textit{prefix}_{N}$ to idle_memory_head and returns the old value of
idle_memory_head to address in an atomic fashion—no data race within this
operation. Once the requested memory is allocated for the block, we increment
the idle_memory_head pointer in the memory pool. We finalize the allocation by
letting all threads in the block compute their corresponding offsets by adding
the prefix sum to the allocated address. The memory allocation is called
inside the meta-kernel that we generated in the operator scheduling. The
entire allocation process has very little overhead costs—it does not require
any inter-block synchronization or any kernel launches.

Algorithm 1 In-Kernel Dynamic Memory Allocation

Input: allocation memory size for the $i^{\textit{th}}$ thread,
$\textit{size}_{i}$; global memory pool head pointer, idle_memory_head;
Output: thread offsets, $\textit{offsets}_{i}$;

1:
$\textit{prefix}_{1..N}\leftarrow\textit{parallel\\_prefix\\_sum}(\textit{size}_{1..N})$

2:
$\textit{address}\leftarrow\textit{atomic\\_add}(\textit{idle\\_memory\\_head},\textit{prefix}_{N})$

3: for each thread $i$ in the current block concurrently do

4:
$\textit{offsets}_{i}\leftarrow\textit{address}+\textit{prefix}_{i}-\textit{prefix}_{1}$

5: end for

TABLE II: End-to-end training of MapReduce feature extraction with hierarchical GPU parameter server and FeatureBox. | Application A | Application B
---|---|---
#Instances | $\sim 1\times 10^{9}$ | $\sim 2\times 10^{9}$
Log Size |  $\sim$15 TB |  $\sim$25 TB
Framework | MapReduce + GPU | FeatureBox | MapReduce + GPU | FeatureBox
#Machines | 20 CPU + 1 GPU | 1 GPU | 30 CPU + 2 GPU | 2 GPU
Execution Time | 18h | 3.5h | 27h | 2.65h
Speedup | - | 5.14X | - | 10.19X
Intermediate I/O Saving | - |  $\sim$50 TB | - |  $\sim$100 TB

Reset GPU memory pool. Our light-weight memory allocation strategy only
maintains a pointer on a pre-allocated continuous global memory. However, the
single-pointer design does not support memory freeing. We have to maintain an
additional collection of freed memory and allocate the requested memory chunks
from this collection—the maintenance of this additional data structure leads
to significant memory allocation overhead. We observe that our operators are
in fine-granularity and are scheduled layer by layer. Therefore, we can assume
that the total required memory for dynamic allocations fits the GPU memory. We
perform the memory release in a batch fashion: the memory pool is reset after
each meta-kernel. The reset can be done in a constant time—we only need to set
idle_memory_head to the original allocated memory address for the memory pool
so that the allocation request in the meta-kernel for the following layer gets
the allocation from the beginning of the memory pool.

##  VI Experimental Evaluation

In this section, we investigate the effectiveness of our proposed framework
FeatureBox through a set of numerical experiments. Specifically, the
experiments are targeted to address the following questions:

  * •

How is the end-to-end training time of FeatureBox compared with the previous
MapReduce solution?

  * •

How much intermediate I/O is saved by the pipelining architecture?

  * •

What is the performance of FeatureBox in the feature extraction task?

Systems. The MapReduce feature extraction baseline is our previous in-
production solution to extract features for the training tasks. It runs in an
MPI cluster with CPU-only nodes in a data center. Commonly, a feature
extraction job requires 20 to 30 nodes. Each node is equipped with server-
grade CPUs ($\sim$100 threads). The training part is executed on GPU nodes.
Each GPU node has 8 cutting-edge 32 GB HBM GPUs, $\sim$1 TB main memory,
$\sim$20 TB RAID-0 NVMe SSDs, and a 100 Gb RDMA network adaptor. The training
framework is the hierarchical GPU parameter server. All nodes are inter-
connected through a high-speed Ethernet switch.

Models. We use CTR prediction models on two real-world online advertising
applications. The neural network backbones of both models follow the design in
Figure 2. The major difference between the two models is the number of input
features. Both models have more than $\sim$10 TB parameters. We collect real
user click history logs as the training dataset.

###  VI-A End-to-End Training

We report Table II specifications about the training data and the end-to-end
training comparison between our proposed FeatureBox and the MapReduce feature
extraction with hierarchical GPU parameter server training as a baseline. Both
training datasets contain billions of instances. The size of the logs is
$\sim$15 TB for application A, and $\sim$25 TB for application B. The end-to-
end training time includes the features extraction from the log time and the
model training time. FeatureBox uses 1 GPU server for application A and 2 GPU
servers for application B. In addition to the GPU servers, the baseline
solution also employs 20/30 CPU-only servers to perform feature extraction.
The baseline solution first extracts features using MapReduce, saves the
features as training data in HDFS, and streams the generated training data to
the GPU servers to train the model. On the other hand, FeatureBox processes
the data in a pipeline fashion: features are extracted on GPU servers and then
are immediately fed to the training framework on the same GPU server. For
application A, FeatureBox only takes 3.5 hours to finish the feature
extraction and the training while the baseline solution requires 18 hours—with
fewer number of machines, FeatureBox has a _5.14X_ speedup compared to the
baseline. Meanwhile, Application B presents a bigger volume of log instances.
Hence, we use two GPU servers to perform the training. We can observe a larger
gap between FeatureBox and the baseline when the data size scales up:
FeatureBox outperforms the baseline with a _10.19X_ speedup. One of the main
reasons of the speedup is that FeatureBox eliminates the huge intermediate I/O
from the MapReduce framework. We save $\sim$50-100 TB intermediate I/O while
using FeatureBox.

###  VI-B Feature Extraction

Although the improvement of FeatureBox in the end-to-end training time mainly
benefits from the pipeline architecture, we also investigate the feature
extraction performance to confirm that our proposed GPU feature extraction
framework is a better alternative to the baseline MapReduce solution.

Figure 6: Feature extraction time of MapReduce and FeatureBox.

We report, in Figure 6, the time to extract features from $10,000$ log
instances of Application B. MapReduce runs on 30 CPU-only servers and
FeatureBox runs on 2 GPU servers. The pre-processing time includes the stages
to prepare the data for the feature extraction, such as read, clean, and join
views. The pre-processing time of both methods are comparable because the
executed operations are mostly memory and network I/O. Regarding the time to
extract features, FeatureBox is more than 3 times faster than MapReduce.
FeatureBox only takes around half of the time to extract the features than the
baseline.

###  VI-C Discussion

Based on these results, we can answer the questions that drive the
experiments: The end-to-end training time of FeatureBox is 5-10 times faster
than the baseline. Due to the pipeline design, FeatureBox saves us 50-100 TB
intermediate I/O. For feature extraction only tasks, FeatureBox on 2 GPU
servers is 2X faster than MapReduce on 30 CPU-only servers.

##  VII Conclusions

In this paper, we introduce FeatureBox, a novel end-to-end training framework
that pipelines the feature extraction and the training on GPU servers to save
the intermediate I/O of the feature extraction. We rewrite computation-
intensive feature extraction operators as GPU operators and leave the memory-
intensive operator on CPUs. We introduce a layer-wise operator scheduling
algorithm to schedule these heterogeneous operators. We present a light-weight
GPU memory management algorithm that supports dynamic GPU memory allocation
with minimal overhead. We experimentally evaluate FeatureBox and compare it
with the previous in-production MapReduce feature extraction framework on two
real-world ads applications. The results show that FeatureBox is 5-10X faster
than the baseline.

## References

  * [1] Dhruba Borthakur.  HDFS architecture guide.  Hadoop apache project, 53(1-13):2, 2008.
  * [2] Heng-Tze Cheng, Levent Koc, Jeremiah Harmsen, Tal Shaked, Tushar Chandra, Hrishi Aradhye, Glen Anderson, Greg Corrado, Wei Chai, Mustafa Ispir, Rohan Anil, Zakaria Haque, Lichan Hong, Vihan Jain, Xiaobing Liu, and Hemal Shah.  Wide & deep learning for recommender systems.  In Proceedings of the 1st Workshop on Deep Learning for Recommender Systems, (DLRS@RecSys), pages 7–10, Boston, MA, 2016.
  * [3] Paul Covington, Jay Adams, and Emre Sargin.  Deep neural networks for youtube recommendations.  In Proceedings of the 10th ACM Conference on Recommender Systems (RecSys), pages 191–198, Boston, MA, 2016.
  * [4] Jeffrey Dean and Sanjay Ghemawat.  Mapreduce: simplified data processing on large clusters.  Communications of the ACM, 51(1):107–113, 2008.
  * [5] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In 9th International Conference on Learning Representations (ICLR), Virtual Event, Austria, 2021.
  * [6] Benjamin Edelman, Michael Ostrovsky, and Michael Schwarz.  Internet advertising and the generalized second-price auction: Selling billions of dollars worth of keywords.  American economic review, 97(1):242–259, 2007.
  * [7] Miao Fan, Jiacheng Guo, Shuai Zhu, Shuo Miao, Mingming Sun, and Ping Li.  MOBIUS: towards the next generation of query-ad matching in baidu’s sponsored search.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD), pages 2509–2517, Anchorage, AK, 2019.
  * [8] Hongliang Fei, Shulong Tan, and Ping Li.  Hierarchical multi-task word embedding learning for synonym prediction.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD), pages 834–842, Anchorage, AK, 2019.
  * [9] Hongliang Fei, Jingyuan Zhang, Xingxuan Zhou, Junhao Zhao, Xinyang Qi, and Ping Li.  GemNN: Gating-enhanced multi-task neural networks with feature interaction learning for CTR prediction.  In Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR), pages 2166–2171, Virtual Event, Canada, 2021.
  * [10] Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron C. Courville, and Yoshua Bengio.  Generative adversarial nets.  In Advances in Neural Information Processing Systems (NIPS), pages 2672–2680, Montreal, Canada, 2014.
  * [11] Thore Graepel, Joaquin Quiñonero Candela, Thomas Borchert, and Ralf Herbrich.  Web-scale bayesian click-through rate prediction for sponsored search advertising in microsoft’s bing search engine.  In Proceedings of the 27th International Conference on Machine Learning (ICML), pages 13–20, Haifa, Israel, 2010.
  * [12] Huifeng Guo, Ruiming Tang, Yunming Ye, Zhenguo Li, and Xiuqiang He.  DeepFM: A factorization-machine based neural network for CTR prediction.  In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence (IJCAI), pages 1725–1731, Melbourne, Australia, 2017.
  * [13] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, Las Vegas, NV, 2016.
  * [14] Tianyu Kang, Ping Chen, John Quackenbush, and Wei Ding.  A novel deep learning model by stacking conditional restricted boltzmann machine and deep neural network.  In Proceedings of the 26th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD), pages 1316–1324, Virtual Event, CA, 2020.
  * [15] Dingcheng Li, Xu Li, Jun Wang, and Ping Li.  Video recommendation with multi-gate mixture of experts soft actor critic.  In Proceedings of the 43rd International ACM SIGIR conference on research and development in Information Retrieval (SIGIR), pages 1553–1556, Virtual Event, China, 2020.
  * [16] Ping Li, Anshumali Shrivastava, Joshua L. Moore, and Arnd Christian König.  Hashing algorithms for large-scale learning.  In Advances in Neural Information Processing Systems (NIPS), pages 2672–2680, Granada, Spain, 2011.
  * [17] Xijun Li, Weilin Luo, Mingxuan Yuan, Jun Wang, Jiawen Lu, Jie Wang, Jinhu Lü, and Jia Zeng.  Learning to optimize industry-scale dynamic pickup and delivery problems.  In Proceedings of the 37th IEEE International Conference on Data Engineering (ICDE), pages 2511–2522, Chania, Greece, 2021.
  * [18] Jianxun Lian, Xiaohuan Zhou, Fuzheng Zhang, Zhongxia Chen, Xing Xie, and Guangzhong Sun.  xdeepfm: Combining explicit and implicit feature interactions for recommender systems.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD), pages 1754–1763, London, UK, 2018\.
  * [19] Bin Liu, Niannan Xue, Huifeng Guo, Ruiming Tang, Stefanos Zafeiriou, Xiuqiang He, and Zhenguo Li.  Autogroup: Automatic feature grouping for modelling explicit high-order feature interactions in CTR prediction.  In Proceedings of the 43rd International ACM SIGIR conference on research and development in Information Retrieval (SIGIR), pages 199–208, Virtual Event, China, 2020.
  * [20] Yifei Ma, Balakrishnan (Murali) Narayanaswamy, Haibin Lin, and Hao Ding.  Temporal-contextual recommendation in real-time.  In Proceedings of the 26th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD), pages 2291–2299, Virtual Event, CA, 2020.
  * [21] Charbel Merhej, Ryan J. Beal, Tim Matthews, and Sarvapali D. Ramchurn.  What happened next? using deep learning to value defensive actions in football event-data.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD), pages 3394–3403, Virtual Event, Singapore, 2021.
  * [22] Giang Nguyen, Stefan Dlugolinsky, Martin Bobák, Viet Tran, Álvaro López García, Ignacio Heredia, Peter Malík, and Ladislav Hluchỳ.  Machine learning and deep learning frameworks and libraries for large-scale data mining: a survey.  Artificial Intelligence Review, 52(1):77–124, 2019.
  * [23] Felipe Pezoa, Juan L. Reutter, Fernando Suárez, Martín Ugarte, and Domagoj Vrgoc.  Foundations of JSON schema.  In Proceedings of the 25th International Conference on World Wide Web (WWW), pages 263–273, Montreal, Canada, 2016.
  * [24] Tao Qi, Fangzhao Wu, Chuhan Wu, and Yongfeng Huang.  PP-Rec: News recommendation with personalized user interest and time-aware news popularity.  In Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (ACL), pages 5457–5467, Virtual Event, 2021.
  * [25] Yanru Qu, Han Cai, Kan Ren, Weinan Zhang, Yong Yu, Ying Wen, and Jun Wang.  Product-based neural networks for user response prediction.  In Proceedings of the 2016 IEEE 16th International Conference on Data Mining (ICDM), pages 1149–1154, Barcelona, Spain, 2016.
  * [26] Ying Shan, T. Ryan Hoens, Jian Jiao, Haijing Wang, Dong Yu, and J. C. Mao.  Deep crossing: Web-scale modeling without manually crafted combinatorial features.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), pages 255–262, San Francisco, CA, 2016.
  * [27] Xiang-Rong Sheng, Liqin Zhao, Guorui Zhou, Xinyao Ding, Binding Dai, Qiang Luo, Siran Yang, Jingshan Lv, Chi Zhang, Hongbo Deng, and Xiaoqiang Zhu.  One model to serve all: Star topology adaptive recommender for multi-domain CTR prediction.  In Proceedings of the 30th ACM International Conference on Information and Knowledge Management (CIKM), pages 4104–4113, Virtual Event, Queensland, Australia, 2021.
  * [28] Athanasios Voulodimos, Nikolaos Doulamis, Anastasios D. Doulamis, and Eftychios Protopapadakis.  Deep learning for computer vision: A brief review.  Comput. Intell. Neurosci., 2018:7068349:1–7068349:13, 2018.
  * [29] Jian Wei, Jianhua He, Kai Chen, Yi Zhou, and Zuoyin Tang.  Collaborative filtering and deep learning based recommendation system for cold start items.  Expert Systems with Applications, 69:29–39, 2017.
  * [30] Kilian Q. Weinberger, Anirban Dasgupta, John Langford, Alexander J. Smola, and Josh Attenberg.  Feature hashing for large scale multitask learning.  In Proceedings of the 26th Annual International Conference on Machine Learning (ICML), pages 1113–1120, Montreal, Canada, 2009.
  * [31] Tan Yu, Zhipeng Jin, Jie Liu, Yi Yang, Hongliang Fei, and Ping Li.  Boost ctr prediction for new advertisements via modeling visual content.  arXiv preprint arXiv:2209.11727, 2022.
  * [32] Siamak Zamany, Dingcheng Li, Hongliang Fei, and Ping Li.  Towards deeper understanding of variational auto-encoders for binary collaborative filtering.  In Proceedings of the 2022 ACM SIGIR International Conference on the Theory of Information Retrieval (ICTIR), pages 254–263, Madrid, Spain, 2022.
  * [33] George Zerveas, Srideepika Jayaraman, Dhaval Patel, Anuradha Bhamidipaty, and Carsten Eickhoff.  A transformer-based framework for multivariate time series representation learning.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (kDD), pages 2114–2124, Virtual Event, Singapore, 2021.
  * [34] Kai Zhang, Hao Qian, Qing Cui, Qi Liu, Longfei Li, Jun Zhou, Jianhui Ma, and Enhong Chen.  Multi-interactive attention network for fine-grained feature learning in CTR prediction.  In Proceedings of the Fourteenth ACM International Conference on Web Search and Data Mining (WSDM), pages 984–992, Virtual Event, Israel, 2021\.
  * [35] Shuai Zhang, Lina Yao, Aixin Sun, and Yi Tay.  Deep learning based recommender system: A survey and new perspectives.  ACM Computing Surveys (CSUR), 52(1):1–38, 2019.
  * [36] Pu Zhao, Chuan Luo, Cheng Zhou, Bo Qiao, Jiale He, Liangjie Zhang, and Qingwei Lin.  RLNF: reinforcement learning based noise filtering for click-through rate prediction.  In Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR), pages 2268–2272, Virtual Event, Canada, 2021.
  * [37] Weijie Zhao, Deping Xie, Ronglai Jia, Yulei Qian, Ruiquan Ding, Mingming Sun, and Ping Li.  Distributed hierarchical GPU parameter server for massive scale deep learning ads systems.  In Proceedings of Machine Learning and Systems 2020 (MLSys), Austin, TX, 2020.
  * [38] Weijie Zhao, Jingyuan Zhang, Deping Xie, Yulei Qian, Ronglai Jia, and Ping Li.  AIBox: CTR prediction model training on a single node.  In Proceedings of the 28th ACM International Conference on Information and Knowledge Management (CIKM), pages 319–328, Beijing, China, 2019\.
  * [39] Guorui Zhou, Xiaoqiang Zhu, Chengru Song, Ying Fan, Han Zhu, Xiao Ma, Yanghui Yan, Junqi Jin, Han Li, and Kun Gai.  Deep interest network for click-through rate prediction.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD), pages 1059–1068, London, UK, 2018\.
  * [40] Jieming Zhu, Jinyang Liu, Weiqi Li, Jincai Lai, Xiuqiang He, Liang Chen, and Zibin Zheng.  Ensembled CTR prediction via knowledge distillation.  In Proceedings of the 29th ACM International Conference on Information and Knowledge Management (CIKM), pages 2941–2958, Virtual Event, Ireland, 2020.