where the last equality is due to the property that

$(z_{i}-1+dy_{i})\cdot dy=0,$ (53)

and the rest are all algebraic manipulations.

Putting together,

$\displaystyle\widehat{G}g(z_{i},\bm{u}^{1:L})$
$\displaystyle=\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(g(z_{i}+a_{i}-dv_{i})-g(z_{i})\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}$
$\displaystyle\mspace{15.0mu}+\lambda_{i}r\left(g(z_{i}-1+dy_{i})-g(z_{i})\right)$
$\displaystyle=\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{{\bm{k}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(2a_{i}\cdot
z_{i}+a_{i}^{2}-(dv_{i})^{2}-2dv_{i}\cdot\eta_{\max}\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}$
$\displaystyle\mspace{15.0mu}+\lambda_{i}r\cdot\left(-2z_{i}+1-(dy_{i})^{2}\right).$

After recombining the terms, we get

$\displaystyle\mspace{23.0mu}\eta_{\max}\cdot\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}dv_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}\right]$
$\displaystyle=\mathbb{E}\left[\left(\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}a_{i}\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}-\lambda_{i}r\right)\cdot
Z_{i}\right]$
$\displaystyle\mspace{23.0mu}+\frac{1}{2}\mathbb{E}\left[\sum_{\ell=1}^{L}\sum_{({\bm{k}},\bm{a}^{\prime})}\gamma_{\widehat{{\bm{K}}}^{\ell},({\bm{k}}^{\prime},\bm{a})}\left(a_{i}^{2}-(dv_{i})^{2}\right)\mathbbm{1}_{\\{\bm{\eta}^{\ell}=\bm{0}\\}}+\lambda_{i}r\cdot(1-(dy_{i})^{2})\right]$

This finishes the calculation.

###  B.4 Proof of Conversion Theorem without Assuming Irreducibility

In this subsection, we prove Theorem 3 without assuming
${\bm{k}}^{0}$-irreducibility of the subroutine $\overline{\sigma}$.
Specifically, suppose that we have a Markovian single-server policy
$\overline{\sigma}$ and an initial distribution $p_{j}$ over its recurrent
classes $S_{j}$ for $j=1,2\dots J$, we will construct an infinite-server
policy $\sigma$ such that (5) and (6) still hold. The basic idea is to
decompose a general Markovian single-server policy $\overline{\sigma}$ into
multiple ${\bm{k}}^{0}$-irreducible Markovian policies, each induces one
recurrent class and preserves stationary distribution and the throughput of
$\overline{\sigma}$ on that recurrent class, as stated in Lemma 5 below.

######  Lemma 5 (Decomposing The Reducible Policy).

Let $\overline{\sigma}$ be a general single-server Markovian policy with
recurrent classes $S_{j}$ for $j=1,2,\dots,J$. Then for each $j$ exists a
Markovian policy $\overline{\sigma}^{j}$ such that

  * •

The induced Markov chain is ${\bm{k}}^{0}$-irreducible with the unique
recurrent class being $S_{j}$;

  * •

The stationary distribution is the same as the stationary distribution under
$\overline{\sigma}$ starting from a configuration in $S_{j}$.

###### Proof.

For each $j=1,2,\dots,J$, we define the policy $\overline{\sigma}^{j}$ as
follows: when the system has configuration ${\bm{k}}\in S_{j}$, the policy
$\overline{\sigma}^{j}$ makes the same decisions as $\overline{\sigma}$; when
${\bm{k}}\notin S_{j}$ and ${\bm{k}}=\bm{0}$, the policy starts a timer whose
duration follows an exponential distribution with rate $1$ and immediately
adds ${\bm{k}}^{0}$ many jobs of each type when the timer ticks for some
arbitrary ${\bm{k}}^{0}\in S_{j}$; when ${\bm{k}}\notin S_{j}$ and
${\bm{k}}\neq\bm{0}$, the policy does not request any jobs.

We show that under the new policy $\overline{\sigma}^{j}$, $S_{j}$ is also a
recurrent class of the induced Markov chain. This is because if the system
starts from a configuration ${\bm{k}}\in S_{j}$, then it will stay in $S_{j}$
since it makes the same decisions and have the same transitions as under the
policy $\overline{\sigma}$. Because $S_{j}$ is a recurrent class under
$\overline{\sigma}$, it is still a recurrent class under
$\overline{\sigma}^{j}$.

To show that the Markov chain induced by $\overline{\sigma}^{j}$ is
${\bm{k}}^{0}$-irreducible, observe that starting from any ${\bm{k}}\notin
S_{j}$, the system state will return to $S_{j}$. Specifically,

  * •

If the system starts from a configuration ${\bm{k}}$ such that ${\bm{k}}\notin
S_{j}$ and ${\bm{k}}\neq\bm{0}$, then no new jobs will be requested until
either ${\bm{k}}\in S_{j}$ or ${\bm{k}}=\bm{0}$. In the latter case, by the
construction of the policy, the system jumps to a configuration in $S_{j}$
after the next transition.

  * •

If the system starts from a configuration ${\bm{k}}$ such that ${\bm{k}}\notin
S_{j}$ and ${\bm{k}}=\bm{0}$, the system jumps to a configuration in $S_{j}$
after the next transition.

The claim that the stationary distribution under $\overline{\sigma}^{j}$ is
the same as the stationary distribution under $\overline{\sigma}$ starting
from a configuration in $S_{j}$ is trivial to show, because when the system
initializes from any configuration in $S_{j}$, it stays in $S_{j}$ and the
transitions are exactly the same under the two policies. ∎

The JRS policy with a general Markovian single-server policy as its subroutine
$\overline{\sigma}$ is constructed using the ${\bm{k}}^{0}$-irreducible
policies $\overline{\sigma}^{j}$’s obtained from the decomposition of
$\overline{\sigma}$ and the probabilities $p^{j}$’s.

  1. 1.

We divide all the servers into $J$ server pools, each with infinite many
servers. Let the $j$-th server pool run the JRS policy with subroutine
$\overline{\sigma}^{j}$ (defined in Section 4.2) for each $j=1,\dots J$.

  2. 2.

Whenever we see an arrival of type $i$, we route the job to the $j$-th
infinite-server system with probability
$\frac{p^{j}\overline{\lambda}_{i}^{j}}{\sum_{j}p^{j}\overline{\lambda}_{i^{\prime}}^{j}}$
for each $j=1,\dots J$.

To analyze the policy $\sigma$, let $\pi^{j}$ and
$\overline{\lambda}_{i}^{j}$’s be the stationary distribution and throughput
of the policy $\overline{\sigma}^{j}$ for $j=1,\dots,J$. By Lemma 5, we have
the following relationships:

$\displaystyle\pi({\bm{k}})$
$\displaystyle=\textstyle\sum_{j=1}^{J}p^{j}\pi^{j}({\bm{k}}),\quad\forall{\bm{k}}\in\mathcal{K},$
(54) $\displaystyle\overline{\lambda}_{i}$
$\displaystyle=\textstyle\sum_{j=1}^{J}p^{j}\overline{\lambda}_{i}^{j},\quad\forall
i\in\mathcal{I}.$ (55)

Based on the above relationships, we can prove the general version of
Conversion Theorem that does not require ${\bm{k}}^{0}$-irreducibility.

###### Proof of Theorem 3..

For each $j=1,\dots,J$, Lemma 5 implies that the policy
$\overline{\sigma}^{j}$ and stationary distribution $\pi^{j}$ form a feasible
solution to the single-server optimization problem
$\overline{\mathcal{P}}\big{(}\big{(}p^{j}\overline{\lambda}_{i}^{j}\overline{N}\big{)}_{i\in\mathcal{I}},\epsilon\big{)}$,
and the corresponding objective value is $p^{j}\overline{N}$. Consider the
infinite-server system with arrival rates
$\big{(}p^{j}\overline{\lambda}_{i}^{j}\overline{N}\big{)}_{i\in\mathcal{I}}$
and budget $\epsilon$. As we have proved Theorem 3 for the JRS policy with
${\bm{k}}^{0}$-irreducible subroutine $\overline{\sigma}^{j}$, it follows that

$\displaystyle\left\lvert\textstyle\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}^{j}]-\left\lceil
p^{j}\overline{N}\right\rceil\cdot(1-\pi^{j}(\bm{0}))\right\rvert$
$\displaystyle=O\left(r^{0.5}\right),$ (56)
$\displaystyle\left\lvert\textstyle\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}^{j}]-\left\lceil
p^{j}\overline{N}\right\rceil\cdot\textstyle\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\pi^{j}({\bm{k}})\right\rvert$
$\displaystyle=O\left(r^{0.5}\right),$ (57)

where $X_{{\bm{k}}}^{j}$ is the random variable representing the steady-state
number of servers in configuration ${\bm{k}}$ in the infinite-server system
under the JRS policy with subroutine $\overline{\sigma}^{j}.$

By the construction of $\sigma$, the arrival rate to the $j$-th server pool is
equal to
$\frac{p^{j}\overline{\lambda}_{i}^{j}}{\sum_{j}p^{j}\overline{\lambda}_{i^{\prime}}^{j}}\cdot\lambda_{i}r=\frac{p^{j}\overline{\lambda}_{i}^{j}}{\overline{\lambda}_{i}}\cdot\lambda_{i}r=p^{j}\overline{\lambda}_{i}^{j}\overline{N}$,
where the first equality is due to (55), and the second equality is due to the
condition that $\lambda_{i}r=\overline{\lambda}_{i}^{j}\cdot L$. Therefore,
(56) and (57) hold, and we have

$\displaystyle\left\lvert\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}]-\overline{N}\cdot(1-\pi(\bm{0}))\right\rvert$
$\displaystyle=\left\lvert\sum_{j=1}^{J}\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}^{j}]-\sum_{j=1}^{J}p^{j}\overline{N}\cdot(1-\pi^{j}(\bm{0}))\right\rvert$
$\displaystyle\leq\sum_{j=1}^{J}\left\lvert\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}^{j}]-\left\lceil
p^{j}\overline{N}\right\rceil\cdot(1-\pi^{j}(\bm{0}))\right\rvert+O(1)$
$\displaystyle=O(r^{0.5}).$

Here we use (56) and the relationship between $\pi({\bm{k}})$ and
$\pi^{j}({\bm{k}})$. Similarly,

$\displaystyle\left\lvert\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\mathbb{E}[X_{{\bm{k}}}]-\overline{N}\cdot\sum_{j=1}^{J}\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\pi({\bm{k}})\right\rvert$
$\displaystyle=\left\lvert\sum_{j=1}^{J}\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}^{j}]-\sum_{j=1}^{J}p^{j}\overline{N}\cdot\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\pi^{j}({\bm{k}})\right\rvert$
$\displaystyle=\sum_{j=1}^{J}\left\lvert\sum_{{\bm{k}}\neq\bm{0}}\mathbb{E}[X_{{\bm{k}}}^{j}]-\left\lceil
p^{j}\overline{N}\right\rceil\cdot\sum_{{\bm{k}}\neq\bm{0}}h({\bm{k}})\pi^{j}({\bm{k}})\right\rvert+O(1)$
$\displaystyle=O(r^{0.5}).$

After re-indexing the servers, we get (5) and (6). The bounds on $N(\sigma)$
and $C(\sigma)$, (7) and (8), follow from (5) and (6). They can be verified in
the same way as that in the proof for the irreducible case, so we omit the
argument here. ∎

##  Appendix C Analysis of Single-Server System

In this section, we solve the single-server optimization problem in (3)
through the linear program (LP) in (9). We first derive (9) as a linear
program relaxation of the single-server optimization problem, and then
construct a single-server policy that achieves the optimal value of the LP,
which implies the optimality of the policy. Moreover, the constructed single-
server policy is a Markovian policy, so it can be converted to an
asymptotically optimal policy for the infinite-server optimization problem
through JRS.

###  C.1 Lower Bound via LP Relaxation

In this subsection we derive an LP relaxation of the optimization problem (3),
restated below:

$\displaystyle\underset{\overline{N},\mspace{10.0mu}\overline{\sigma},\mspace{10.0mu}\pi}{\text{minimize}}$
$\displaystyle\overline{N}$ (3 revisited) subject to
$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]\leq\epsilon,$
$\displaystyle\overline{N}\cdot\overline{\lambda}_{i}=\lambda_{i}r,\quad\forall
i\in\mathcal{I}.$

Observe that both $\overline{{\bm{K}}}(\infty)$ and $\overline{\lambda}_{i}$
depend on the stationary distribution $\pi$, but the constraints in terms of
$\pi$ are implicit. To derive an LP relaxation, we give an explicit
characterization of the constraints that must be satisfied by the stationary
distribution $\pi$ induced by any feasible policy $\overline{\sigma}$ (not
necessarily Markovian).

To do this, we derive a version of the stationary equation in terms of a
quantity called transition frequency. The transition frequency of type $i$
jobs is a function $u_{i}\colon\mathcal{K}\to\mathbb{R}$ describing the
steady-state frequency of requesting a type $i$ job when the system has
configuration $\bm{k}$. To rigorously define transition frequency, we first
introduce a concept called nominal transition.

######  Definition 1 (Nominal Transition).

Consider a single-server system under any policy. When the configuration
$\overline{{\bm{K}}}(t)$ transitions from ${\bm{k}}$ to
${\bm{k}}^{\prime}+\bm{a}$ for some ${\bm{k}},{\bm{k}}^{\prime}\in\mathcal{K}$
with $\bm{a}=(a_{i})_{i\in\mathcal{I}}$ new jobs added into service, we
decompose the transition by adding intermediate configurations as illustrated
below, where ${\bm{k}}$ first goes to ${\bm{k}}^{\prime}$ if
${\bm{k}}^{\prime}\neq{\bm{k}}$, then add jobs of each type one by one.

${\bm{k}}\to{\bm{k}}^{\prime}\to({\bm{k}}^{\prime}+\bm{e}_{i_{1}})\to\dots\to({\bm{k}}^{\prime}+a_{i_{1}}\bm{e}_{i_{1}})\to\dots\to({\bm{k}}^{\prime}+a_{i_{1}}\bm{e}_{i_{1}}+\dots+a_{i_{\left\lvert\mathcal{I}\right\rvert}}\bm{e}_{i_{\left\lvert\mathcal{I}\right\rvert}}),$

where $(i_{1},i_{2},\dots,i_{\left\lvert\mathcal{I}\right\rvert})$ is a fixed
ordering of the set of phases $\mathcal{I}$. We call each short transition in
the diagram a nominal transition.

For ${\bm{k}}^{1},{\bm{k}}^{2}\in\mathcal{K}$, we denote
$F({\bm{k}}^{1},{\bm{k}}^{2},t)$ as the cumulative number of nominal
transitions from ${\bm{k}}^{1}$ to ${\bm{k}}^{2}$ during the time interval
$[0,t]$, which is a random variable with a distribution depending on the
single-server policy and initial distribution of configurations.

Note that for any $i\in\mathcal{I}$ and ${\bm{k}}\in\mathcal{K}$ s.t.
$k_{i}\geq 1$, $F({\bm{k}},{\bm{k}}-\bm{e}_{i},t)$ counts the number of times
that a type $i$ job departs when being in configuration ${\bm{k}}$ . As a
result,

$F({\bm{k}},{\bm{k}}-\bm{e}_{i},t)=\mathcal{N}\left(\int_{0}^{t}k_{i}\mu_{i\perp}\mathbbm{1}_{\\{K(s)={\bm{k}}\\}}ds\right),$

where $\mathcal{N}(t)$ denotes a unit rate Poisson process. If we take
expectation, divide both sides by $t$, and let $t\to\infty$, we have

$\lim_{t\to\infty}\frac{1}{t}\mathbb{E}[F({\bm{k}},{\bm{k}}-\bm{e}_{i},t)]=k_{i}\mu_{i\perp}\cdot\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}\mathbb{P}({\bm{K}}(s)={\bm{k}})ds=k_{i}\mu_{i\perp}\pi({\bm{k}}).$
(58)

Similarly, $F({\bm{k}},{\bm{k}}-\bm{e}_{i}+\bm{e}_{i^{\prime}},t)$ counts the
number of times a job in phase $i$ transitions to phase $i^{\prime}$ when
being in configuration ${\bm{k}}$ for any $i,i^{\prime}\in\mathcal{I}$,
${\bm{k}}\in\mathcal{K}$ s.t. $i^{\prime}\neq i$ and $k_{i}\geq 1$, so

$\lim_{t\to\infty}\frac{1}{t}\mathbb{E}[F({\bm{k}},{\bm{k}}-\bm{e}_{i}+\bm{e}_{i^{\prime}},t)]=k_{i}\mu_{ii^{\prime}}\pi({\bm{k}}).$
(59)

We define transition frequency as follows.

######  Definition 2 (Transition Frequency).

Transition frequency of type $i$ jobs at state ${\bm{k}}$ is the long-run
average number of nominal transitions from configuration ${\bm{k}}$ to
${\bm{k}}+\bm{e}_{i}$ per unit time,

$u_{i}({\bm{k}})\triangleq\lim_{t\to\infty}\frac{1}{t}\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)].$
(60)

The transition frequencies allow us to derive the following version of
stationary equation.

######  Lemma 6 (Stationary Equation).

Under any policy, the stationary distribution $\pi$ and the transition
frequency $u_{i}$ satisfy the following equation:

$\displaystyle\mspace{23.0mu}\sum_{i}u_{i}({\bm{k}}-\bm{e}_{i})\mathbbm{1}_{\\{k_{i}\geq
1\\}}+\sum_{i}(k_{i}+1)\mu_{i\perp}\pi({\bm{k}}+\bm{e}_{i})$ (61)
$\displaystyle\mspace{23.0mu}+\sum_{i,i^{\prime}\colon i\neq
i^{\prime}}(k_{i}+1)\mu_{ii^{\prime}}\pi({\bm{k}}+\bm{e}_{i}-\bm{e}_{i^{\prime}})\mathbbm{1}_{\\{k_{i^{\prime}}\geq
1\\}}$
$\displaystyle=\sum_{i}u_{i}({\bm{k}})+\Big{(}\sum_{i}k_{i}\mu_{i\perp}+\sum_{i,i^{\prime}\colon
i\neq i^{\prime}}k_{i}\mu_{ii^{\prime}}\Big{)}\pi({\bm{k}})$

for any state ${\bm{k}}\in\mathcal{K}$, and $\sum_{i},\sum_{i,i^{\prime}}$ are
shorthand for $\sum_{i\in\mathcal{I}},\sum_{i,i^{\prime}\in\mathcal{I}}$.

###### Proof.

For each configuration ${\bm{k}}\in\mathcal{K}$, if we look at difference
between the number of nominal transitions into configuration ${\bm{k}}$ and
that out of configuration ${\bm{k}}$ by time $t$, we have the following
equation,

$\displaystyle\mspace{23.0mu}\mathbbm{1}_{\\{\overline{{\bm{K}}}(t)={\bm{k}}\\}}-\mathbbm{1}_{\\{\overline{{\bm{K}}}(0)={\bm{k}}\\}}$
(62)
$\displaystyle=\sum_{i}F({\bm{k}}-\bm{e}_{i},{\bm{k}},t)\mathbbm{1}_{\\{k_{i}\geq
1\\}}+\sum_{i}F({\bm{k}}+\bm{e}_{i},{\bm{k}},t)+\sum_{i,i^{\prime}\colon i\neq
i^{\prime}}F({\bm{k}}+\bm{e}_{i}-\bm{e}_{i^{\prime}},{\bm{k}},t)\mathbbm{1}_{\\{k_{i^{\prime}}\geq
1\\}}$
$\displaystyle\mspace{23.0mu}-\sum_{i}F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)-\sum_{i}F({\bm{k}},{\bm{k}}-\bm{e}_{i},t)\mathbbm{1}_{\\{k_{i}\geq
1\\}}-\sum_{i,i^{\prime}\colon i\neq
i^{\prime}}F({\bm{k}},{\bm{k}}-\bm{e}_{i}+\bm{e}_{i^{\prime}},t)\mathbbm{1}_{\\{k_{i}\geq
1\\}}$

By Definition 2 and (58)-(59), if we divide both sides of (62) by $t$ and let
$t\to\infty$, we get the stationary equation in (61). ∎

Since the stationary equation in (61) is linear in $u_{i}({\bm{k}})$ and
$\pi({\bm{k}})$, we can write it in matrix form:

$A\bm{\pi}+\textstyle\sum_{i\in\mathcal{I}}B_{i}\bm{u}_{i}=0,$ (63)

where $\bm{\pi}$ and $\bm{u}_{i}$ are column vectors representing
$\pi(\cdot)$, $u_{i}(\cdot)$, and $A,B_{i}$ are matrices that make (63)
equivalent to (61). Therefore, the following three conditions are necessary
for any tuple $(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ to be a possible
pair of stationary distribution and transition frequencies for a Markovian
policy.

$\displaystyle A\bm{\pi}+\textstyle\sum_{i\in\mathcal{I}}B_{i}\bm{u}_{i}$
$\displaystyle=0$ (64) $\displaystyle\textstyle\sum_{{\bm{k}}}\pi({\bm{k}})$
$\displaystyle=1$ $\displaystyle\bm{\pi},\bm{u}_{i}$ $\displaystyle\geq
0,\quad\forall i\in\mathcal{I}$

Now we can convert the single-server control problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ in (3) to
the linear program
$\overline{\mathcal{LP}}((\lambda_{i})_{i\in\mathcal{I}},\epsilon)$ defined in
(9).

Note that (63) is the necessary condition for any policy, thus any feasible
solution to (3) is also feasible to (9). Let $\overline{N}^{*}$ be the optimal
value of (3) and $\Phi^{*}$ be the optimal value of (9). We have

$\overline{N}^{*}\geq\frac{r}{\Phi^{*}}.$ (65)

###  C.2 Policy Construction

In this subsection, we describe a procedure that allows us to construct a
policy that achieves the lower bound given by the LP relaxation in (9).
Specifically, given a feasible solution
$(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ to (9), we define an LP-based
policy that requests jobs as follows:

  * •

Case 1: When the system enters a configuration ${\bm{k}}$ with
$\pi({\bm{k}})\neq 0$, for each $i\in\mathcal{I}$, the policy starts a timer
whose duration follows an exponential distribution with rate
$\frac{u_{i}({\bm{k}})}{\pi({\bm{k}})}$. The policy requests a type $i$ job
when the the $i$-th timer ticks. When the configuration changes, all timers
are canceled.

  * •

Case 2: When the system enters a configuration ${\bm{k}}$ with
$\pi({\bm{k}})=0$ and $\sum_{i^{\prime}}u_{i^{\prime}}({\bm{k}})\neq 0$, the
policy immediately requests a type $i$ job with probability
$\frac{u_{i}({\bm{k}})}{\sum_{i^{\prime}}u_{i^{\prime}}({\bm{k}})}$.

  * •

Case 3: When the system enters a configuration ${\bm{k}}$ with
$\pi({\bm{k}})=0$ and $\sum_{i^{\prime}}u_{i^{\prime}}({\bm{k}})=0$, the
policy does not request any jobs.

We denote the LP-based policy based on the solution
$(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ as
$\overline{\sigma}(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$.

######  Remark 2.

Note that the definition of the LP-based policy here is stated from a view
different from the view in Section 4.2: here each request only adds one job to
the server, and one request can happen immediately after another; while in
Section 4.2 each request can add multiple jobs to the server, and there is
only one request happening at a time. We refer to the view here as the
_impulsive view_ , because multiple requests happening at the same time can be
thought of as having an infinite request rate. In contrast, we call the view
in Section 4.2 the _non-impulsive view_.

The LP-based policy can be alternatively described using the non-impulsive
view if we see multiple requests happening at the same time as one request
that adds multiple jobs to the server. More specifically, each reactive
request of the LP-based policy is initiated by an internal transition or
departure and consists of one or multiple requests of Case 2; each proactive
request of the LP-based policy consists of one request in Case 1 and possibly
several requests of Case 2.

The following lemma characterizes the steady-state behavior of a single-server
system under the LP-based policy.

######  Lemma 7 (Properties of LP-based Policies).

Consider a single-server system under the LP-based policy
$\overline{\sigma}(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$, where
$(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ is a feasible solution to (9). We
have that $\bm{\pi}$ is a stationary distribution under policy
$\overline{\sigma}$, and $(\bm{u}_{i})_{i\in\mathcal{I}}$ are the transition
frequencies corresponding to $\overline{\sigma}$.

The proof of Lemma 7 is based on (62), following the same argument as the
proof of Lemma 6, as well as an induction argument.

###### Proof.

Let $(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ be a feasible solution to the
LP in (9). To show that $\bm{\pi}$ and $(\bm{u}_{i})_{i\in\mathcal{I}}$ are
also the actual stationary distribution and transition frequencies, it
suffices to show that if the initial distribution follows $\bm{\pi}$, i.e.

$\mathbb{P}\big{(}\overline{{\bm{K}}}(0)={\bm{k}}\big{)}=\pi({\bm{k}})$

then under the policy
$\overline{\sigma}(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$, we have

$\displaystyle\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathbb{P}\big{(}\overline{{\bm{K}}}(t)={\bm{k}}\big{)}dt=\pi({\bm{k}}),$
(66)
$\displaystyle\lim_{T\to\infty}\frac{1}{T}\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},T)]=u_{i}({\bm{k}}),$
(67)

where $F$ is the cumulative number of nominal transitions under the policy
$\overline{\sigma}(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$ and the initial
distribution $\bm{\pi}$.

Our proof is based on the following equation:

$\displaystyle\mspace{23.0mu}\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})\Big{|}_{t=0}$
(68)
$\displaystyle=\sum_{i}\frac{d}{dt}\mathbb{E}\left[F({\bm{k}}-\bm{e}_{i},{\bm{k}},t)\right]\Big{|}_{t=0}-\sum_{i}\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}$
$\displaystyle\mspace{23.0mu}+\sum_{i}(k_{i}+1)\mu_{i\perp}\pi({\bm{k}}+\bm{e}_{i})+\sum_{i,i^{\prime}}(k_{i}+1)\mu_{ii^{\prime}}\pi({\bm{k}}+\bm{e}_{i}-\bm{e}_{i^{\prime}})\mathbbm{1}_{\\{k_{i^{\prime}}\geq
1\\}}$
$\displaystyle\mspace{23.0mu}-\sum_{i}k_{i}\mu_{i\perp}\pi({\bm{k}})-\sum_{i,i^{\prime}}k_{i}\mu_{ii^{\prime}}\pi({\bm{k}}),$

The equation is a straightforward consequence of (62), following the same
argument as the proof of Lemma 6.

We prove the following two equations by induction on
$\sum_{i\in\mathcal{I}}k_{i}$.

$\displaystyle\lim_{t\to
0}\frac{1}{t}\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)]$
$\displaystyle=u_{i}({\bm{k}}),$ (69)
$\displaystyle\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})|_{t=0}$
$\displaystyle=0.$ (70)

We first consider the base case when $\sum_{i\in\mathcal{I}}k_{i}=0$. In this
case, ${\bm{k}}=\bm{0}$ and we have

$\frac{d}{dt}\mathbb{E}\left[F({\bm{k}}-\bm{e}_{i},{\bm{k}},t)\right]\Big{|}_{t=0}=u_{i}({\bm{k}}-\bm{e}_{i})\mathbbm{1}_{\\{k_{i}\geq
1\\}}=0,$ (71)

for all $i$. This reduces (68) to

$\displaystyle\mspace{23.0mu}\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})\Big{|}_{t=0}$
(72) $\displaystyle=\sum_{i}u({\bm{k}}-\bm{e}_{i})\mathbbm{1}_{\\{k_{i}\geq
1\\}}-\sum_{i}\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}$
$\displaystyle\mspace{23.0mu}+\sum_{i}(k_{i}+1)\mu_{i\perp}\pi({\bm{k}}+\bm{e}_{i})+\sum_{i,i^{\prime}}(k_{i}+1)\mu_{i,i^{\prime}}\pi({\bm{k}}+\bm{e}_{i}-\bm{e}_{i^{\prime}})\mathbbm{1}_{\\{k_{i^{\prime}}\geq
1\\}}$
$\displaystyle\mspace{23.0mu}-\sum_{i}k_{i}\mu_{i\perp}\pi({\bm{k}})-\sum_{i,i^{\prime}}k_{i}\mu_{i,i^{\prime}}\pi({\bm{k}}),$

Now we discuss based on whether $\pi({\bm{k}})=0$. If $\pi({\bm{k}})\neq 0$,
by the definition of our policy, for all $i$,

$\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}=\frac{u_{i}({\bm{k}})}{\pi({\bm{k}})}\cdot\mathbb{P}(\overline{{\bm{K}}}(0)={\bm{k}})=u_{i}({\bm{k}}),$
(73)

which is (69). Combining the above equation and the stationary equation (61)
satisfied by $(\bm{\pi},(\bm{u}_{i})_{i\in\mathcal{I}})$, we conclude that the
RHS of (72) is zero, i.e.

$\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})|_{t=0}=0,$

which is (70). For the case when $\pi({\bm{k}})=0$ and
$\sum_{i}u_{i}({\bm{k}})\neq 0$, because the system immediately leave the
configuration ${\bm{k}}$ after reaching it through a nominal transition,

$\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})|_{t=0}=0,$

i.e., the LHS of (72) is $0$. Again we compare (72) against the stationary
equation (61) and get

$\sum_{i}\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}=\sum_{i}u_{i}({\bm{k}}).$

By the definition of our policy, we have

$\displaystyle\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}$
$\displaystyle=\frac{u_{i}({\bm{k}})}{\sum_{i^{\prime}}u_{i^{\prime}}({\bm{k}})}\cdot\sum_{i^{\prime}}\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i^{\prime}},t)\right]\Big{|}_{t=0}=u_{i}({\bm{k}}).$
(74)

which is (69). For the case when $\pi({\bm{k}})=0$ and
$\sum_{i}u_{i}({\bm{k}})=0$, (61) implies that $u_{i}({\bm{k}}-\bm{e}_{i})=0$,
$\pi({\bm{k}}+\bm{e}_{i})=0$, $\pi({\bm{k}}+\bm{e}_{i}-\bm{e}_{i^{\prime}})=0$
for any $i$. Then (72) is further reduced to

$\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})\Big{|}_{t=0}=-\sum_{i}\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}.$

Because $\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})\geq 0$ and
$\mathbb{P}(\overline{{\bm{K}}}(0)={\bm{k}})=0$, the LHS of the above
expression is non-negative. However, the RHS of the above expression is non-
positive. Therefore, both sides equal to zero, thus we have
$\frac{d}{dt}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})\Big{|}_{t=0}=0$ and
$\frac{d}{dt}\mathbb{E}\left[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)\right]\Big{|}_{t=0}=u_{i}({\bm{k}})$.

Having proved the base case, we do the induction step. Suppose we have proved
(69) and (70) for all ${\bm{k}}$ such that $\sum_{i\in\mathcal{I}}k_{i}\leq
m-1$ for some integer $m\geq 1$. We consider ${\bm{k}}$ with
$\sum_{i\in\mathcal{I}}k_{i}=m$. By the induction hypothesis,

$\frac{d}{dt}\mathbb{E}\left[F({\bm{k}}-\bm{e}_{i},{\bm{k}},t)\right]\Big{|}_{t=0}=u_{i}({\bm{k}}-\bm{e}_{i})\mathbbm{1}_{\\{k_{i}\geq
1\\}}.$ (75)

Then we repeat the arguments after (71) of the base case verbatim. By
induction, we have proved (69) and (70).

Therefore, given the policy and initial distribution, the distribution of
$\overline{{\bm{K}}}(t)$ is stationary, i.e., we always have
$\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})=\pi({\bm{k}})$ for all
${\bm{k}}\in\mathcal{K}$. As a result, an analogue of (69) holds for all
$t\geq 0$: $\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)]$ is differentiable
with respect to $t$ and

$\displaystyle\frac{d}{dt}\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)]$
$\displaystyle=u_{i}({\bm{k}}),$

for all ${\bm{k}}\in\mathcal{K}$ and all $i\in\mathcal{I}$. Therefore,

$\displaystyle\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathbb{P}(\overline{{\bm{K}}}(t)={\bm{k}})dt$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\cdot\pi({\bm{k}})\cdot
T=\pi({\bm{k}}),$
$\displaystyle\lim_{T\to\infty}\frac{1}{T}\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},T)]$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\frac{d}{dt}\mathbb{E}[F({\bm{k}},{\bm{k}}+\bm{e}_{i},t)]dt$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\cdot u_{i}({\bm{k}})\cdot
T=u_{i}({\bm{k}}).$

This completes the proof. ∎

###  C.3 Proof of Theorem 4

###### Proof.

By Lemma 7, under the policy
$\overline{\sigma}(\bm{\pi}^{*},(\bm{u}^{*}_{i})_{i\in\mathcal{I}})$,
$\bm{\pi}^{*}$ is a stationary distribution, and
$(\bm{u}^{*}_{i})_{i\in\mathcal{I}})$ are the corresponding transition
frequencies. Recall the single-server optimization
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$ problem

$\displaystyle\underset{\overline{N},\mspace{10.0mu}\overline{\sigma},\mspace{10.0mu}\pi}{\text{minimize}}$
$\displaystyle\overline{N}$ (3 revisited) subject to
$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]\leq\epsilon,$
$\displaystyle\overline{N}\cdot\overline{\lambda}_{i}=\lambda_{i}r,\quad\forall
i\in\mathcal{I}.$

Observe that under the policy
$\overline{\sigma}(\bm{\pi}^{*},(\bm{u}^{*}_{i})_{i\in\mathcal{I}})$, the cost
of resource contention is $\bm{h}^{T}\bm{\pi}\leq\epsilon(1-\pi(\bm{0}))$, the
request rate of type $i$ jobs is
$\overline{\lambda}_{i}=\bm{1}_{o}^{T}\bm{u}^{*}_{i}=\Phi^{*}\cdot\lambda_{i}$,
so

$\displaystyle\mathbb{E}\left[h\big{(}\overline{{\bm{K}}}(\infty)\big{)}\middle|\overline{{\bm{K}}}(\infty)\neq\bm{0}\right]=\frac{\bm{h}^{T}\bm{\pi}}{1-\pi(\bm{0})}\leq\epsilon,$
$\displaystyle\overline{\lambda}_{i}=\Phi^{*}\cdot\lambda_{i},\quad\forall
i\in\mathcal{I}.$

where we have used the fact that $h(\bm{0})=0$ in the first equality.
Therefore,
$(\Phi^{*}/r,\overline{\sigma}(\bm{\pi}^{*},(\bm{u}^{*}_{i})_{i\in\mathcal{I}}),\bm{\pi})$
is a feasible solution to the single-server problem
$\overline{\mathcal{P}}((\lambda_{i}r)_{i\in\mathcal{I}},\epsilon)$, achieving
the objective value of $r/\Phi^{*}$, which is the optimal value because
$r/\Phi^{*}\leq\overline{N}^{*}$. ∎


# Many-to-One Distribution Learning and $K$-Nearest Neighbor Smoothing for
Thoracic Disease Identification

Yi Zhou1, Lei Huang2, Tianfei Zhou3, Ling Shao4
Corresponding author: Yi Zhou

###### Abstract

Chest X-rays are an important and accessible clinical imaging tool for the
detection of many thoracic diseases. Over the past decade, deep learning, with
a focus on the convolutional neural network (CNN), has become the most
powerful computer-aided diagnosis technology for improving disease
identification performance. However, training an effective and robust deep CNN
usually requires a large amount of data with high annotation quality. For
chest X-ray imaging, annotating large-scale data requires professional domain
knowledge and is time-consuming. Thus, existing public chest X-ray datasets
usually adopt language pattern based methods to automatically mine labels from
reports. However, this results in label uncertainty and inconsistency. In this
paper, we propose many-to-one distribution learning (MODL) and $K$-nearest
neighbor smoothing (KNNS) methods from two perspectives to improve a single
model’s disease identification performance, rather than focusing on an
ensemble of models. MODL integrates multiple models to obtain a soft label
distribution for optimizing the single target model, which can reduce the
effects of original label uncertainty. Moreover, KNNS aims to enhance the
robustness of the target model to provide consistent predictions on images
with similar medical findings. Extensive experiments on the public NIH Chest
X-ray and CheXpert datasets show that our model achieves consistent
improvements over the state-of-the-art methods.

## Introduction

Chest X-rays are one of the most common radiology exams. Hundreds of millions
of such images are acquired in hospitals and clinics all over the world, for
identifying a wide range of diseases. Advances in deep learning present a
powerful means of developing automated systems to help radiologists interpret
chest X-ray images. Recent efforts have shown promise in improving the
identification of different chest lesions (Rajpurkar et al. 2017; Allaouzi and
Ahmed 2019), localization of their corresponding positions and sizes (Li et
al. 2018b; Liu et al. 2019), and generation of human-readable diagnostic
reports (Jing, Xie, and Xing 2018; Li et al. 2018a; Zhang et al. 2020).
However, there exist many drawbacks to chest X-rays, such as the lack of 3D
information, inconsistent visual characteristics of various lesions, and label
uncertainty, which still make automated chest X-ray diagnosis a very
challenging task.

Figure 1: (a) An input image. (b) The ground-truth binary labels automatically
extracted from reports. (c) Random smoothed labels proposed by (Pham et al.
2019). (d) Our many-to-one soft label distribution. (e) The bounding boxes of
three different colors are annotated by three radiologists, separately, in the
left image. The right image is the heat map, showing abnormal regions,
obtained by our method.

For chest X-ray images, large anonymized public datasets are available to
researchers, and have facilitated the development of deep learning models.
However, providing accurate clinical labels for the very large image sets
required for training deep models is difficult and time-consuming. Most
datasets (Wang et al. 2017; Irvin et al. 2019) have applied rule-based natural
language processing (NLP) to radiology reports to automatically mine the
labels, but this often leads to inconsistencies and noise for training.
Although previous efforts have tried to address this problem, for example, by
converting all the uncertain labels into positive or negative (Irvin et al.
2019), or applying label smoothing regularization (Pham et al. 2019), these
methods are not quite effective and only bring a slight improvement. Moreover,
most efforts focus on model ensembles when pursuing a further increase of
accuracy, rather than enhancing a single model. Combining multiple models,
however, will largely increase the computational costs and complexities.
Therefore, in this paper, we aim to address these two challenges by exploring
label distributions to reduce the effects of uncertainty and noise, and
enhancing a single model’s ability rather than using a model ensemble. Figure
1 illustrates our motivation and solutions.

To train deep learning models for identifying thoracic diseases, multi-label
supervision indicating different lesions is usually adopted. The optimization
of such models is treated as an image classification problem that builds a
mapping from the instances to the labels via a Softmax cross-entropy or binary
cross-entropy loss for each label. However, binary labels (positive or
negative) cannot usually describe in detail the severity of a certain disease
in an image. As shown in Figure 1 (e), even radiologists are unable to always
provide consistent annotations. With the label uncertainty and error further
introduced by automatic label extraction, we believe using soft labels of
probabilities will contribute to more reasonable training. Thus, in addition
to multi-label classification, we propose to optimize the Kullback-Leibler
divergence between the ground-truth soft label distribution and the predicted
one. Moreover, to obtain a single model with competitive performance to model
ensembles, we compute $K$ nearest neighbors using the soft label distributions
integrated from multiple models, and then constrain a single target model to
provide similar predictions to these scans, which can enhance its robustness.
Our main contributions are highlighted as follows:

1. A many-to-one distribution learning (MODL) method is proposed to exploit multiple reference models and integrate their predictions as a soft label distribution to optimize the target network, which can reduce the effects of label uncertainty and ambiguity.

2. To enhance the robustness of a target model for providing consistent predictions on images with similar medical findings, we introduce a $K$-nearest neighbor smoothing (KNNS) module as an auxiliary loss for optimization. The local similarity score is adopted to enhance the learning of the neighbor smoothing.

3. Extensive experiments are conducted to demonstrate the effectiveness of the proposed MODL and KNNS. Consistent improvements by our method are achieved on both the NIH Chest X-ray and CheXpert datasets for thoracic disease identification, compared to state-of-the-art models. Another big advantage of our method is that our single target model does not require additional computational resources during the testing phase.

## Related Work

### Thoracic Disease Diagnosis

Recently, deep convolutional neural networks (CNNs) have been extensively
investigated for thoracic disease diagnosis (Litjens et al. 2017; Shin et al.
2016; Zhou et al. 2021; Fan et al. 2020). The common tasks can be categorized
into multi-disease identification (Kumar, Grewal, and Srivastava 2018;
Rajpurkar et al. 2017; Zhou, Li, and Wang 2018; Baltruschat et al. 2019),
localization (Pesce et al. 2019; Yao et al. 2018; Liu et al. 2019), and
automatic report generation (Wang et al. 2018; Jing, Xie, and Xing 2018; Li et
al. 2018a). To facilitate research into thoracic disease diagnosis, various
public chest X-ray datasets (Wang et al. 2017; Irvin et al. 2019; Johnson et
al. 2019) have been released. Several single- or multi-CNN models (Rajpurkar
et al. 2017) have been developed for disease classification, while weakly
supervised mechanisms using limited bounding-box annotations are typically
used for localization. Li et al. (Li et al. 2018b) proposed a patch slicing
layer to resize CNN features by max-pooling or bilinear interpolation and
applied a fully convolutional recognition network on these regions, to improve
identification and localization, simultaneously. Moreover, a contrast induced
attention network was introduced by (Liu et al. 2019) to exploit the highly
structured property of chest X-ray images for localizing diseases. In
addition, the automatic generation of medical image reports has also been
investigated; for example, (Jing, Xie, and Xing 2018) exploited a co-attention
mechanism to identify abnormal regions and used a hierarchical LSTM to
generate sentences. However, while several previous works have mentioned
issues with label uncertainty and inconsistency of annotation, but none have
provided an effective solution.

### Reducing Label Inconsistency and Ambiguity

Label inconsistency, uncertainty, and ambiguity are common problems that
affect the training of machine learning models in numerous areas. Several
strategies, such as confident learning (CL) (Northcutt, Jiang, and Chuang
2019), label distribution learning (LDL) (Geng 2016; Gao et al. 2017; Chen et
al. 2020), label enhancement (LE) (Xu, Lv, and Geng 2019; Xu, Tao, and Geng
2018), and active learning (AL) (Wu et al. 2018), have been proposed in recent
years to mitigate these problems. CL is based on the principles of pruning
noisy data, counting to estimate noise, and ranking examples to train with
confidence. It can directly estimate the joint distribution of noisy and true
labels, and find the label errors. LDL treats the labels representing the
degree to which each label describes the instance. It is a more general
learning framework of which single-label and multi-label classification are
special cases. LDL has been adopted in many applications, including facial
landmark detection (Su and Geng 2019) and age estimation (Hou et al. 2017). LE
methods aim to solve the unavailability of label distributions. For example, a
partial label learning method was presented by (Xu, Lv, and Geng 2019), in
which label distributions are recovered by leveraging the topological
information of the feature space. Moreover, AL adopts a small set of clean
data, and decreases the noise of the training set by re-labeling other
uncertain samples.

Figure 2: Overview of the proposed method. Network-T is the single target
network, which we aim to optimize. $N$ reference networks are pre-trained to
compute soft label distributions for helping to optimize Network-T with a
distribution learning loss, which can reduce the uncertainty and ambiguity of
the original binary labels. Moreover, $K$ nearest neighbors are computed based
on the soft label distribution, and then used to constrain the Network-T to
provide consistent predictions using an auxiliary neighbor smoothing loss.

## Proposed Methods

### Problem Formulation

First of all, the main notations used in this paper are defined as follows.
Given the $i$-th input X-ray image $x_{i}$ and its corresponding binary ground
truth labels $y^{j}_{i}$, where $j$ denotes the $j$-th disease label, we aim
to optimize a single target network $f_{T}$ for thoracic disease
identification. Due to the uncertainty of the original binary labels, we train
multiple reference models $f_{n},n\in\\{1,2,...,N\\}$, where $N$ is the number
of networks, and integrate their pseudo probability predictions as the input
image’s soft label distribution
$l_{i}=(l_{x_{i}}^{y_{1}},l_{x_{i}}^{y_{2}},...,l_{x_{i}}^{y_{C}})$, where $C$
is the number of label classes. Thus, in addition to the standard multi-label
classification loss $\mathcal{L}_{cls}$, we propose a label distribution
learning loss $\mathcal{L}_{distri}$ to better optimize the target $f_{T}$ for
reducing the effects of the label uncertainty. Moreover, we indicate the input
image as an anchor, and compute its $K$ nearest neighbors
$x_{i}^{k},k\in\\{1,2,...,K\\}$ based on the averaged soft label predictions.
An auxiliary neighbor smoothing loss $\mathcal{L}_{neigh}$ is introduced to
enhance the robustness of the target $f_{T}$. Therefore, the target model can
be optimized by minimizing the overall loss function as:

$\displaystyle\mathcal{L}=\sum_{i,j}[\mathcal{L}_{cls}(y_{i}^{j},f_{T}(x_{i}|\theta))+\lambda\mathcal{L}_{distri}(l_{i},f_{T}(x_{i}|\theta))$
(1)
$\displaystyle+\gamma\sum_{k}\mathcal{L}_{neigh}(f_{T}(x_{i}|\theta),f_{T}(x_{i}^{k}|\theta))]\text{,}$

where $\theta$ is the parameter vector, and $\lambda$ and $\gamma$ balance the
weights of different losses. An overview of our proposed method is illustrated
in Figure 2.

### Basic Multi-Label Classification and Label Smoothing Regularization

The basic thoracic disease identification framework used in most previous
works (Wang et al. 2017; Rajpurkar et al. 2017), usually adopts a classic deep
CNN (i.e. backbone) and optimizes it with a binary cross-entropy (BCE) loss
for multi-label classification. DenseNet-121 (Huang et al. 2017) is the most
commonly used backbone for this task, since it improves the flow of
information and gradients through the network, making the optimization of deep
networks tractable. Therefore, we select DenseNet-121 as our target network
$f_{T}$ for final evaluation, and use the BCE loss as $\mathcal{L}_{cls}$ for
training.

To deal with the label uncertainty caused by the ground-truth labels not being
annotated by human experts for images provided in the training dataset,
several preliminary studies have been conducted. For example, for the large-
scale public CheXpert (Irvin et al. 2019) dataset, uncertain labels are often
extracted due to the unavoidable ambiguities in radiology reports. These
labels can all be $ignored$ (U-Ignore), converted to $positive$ labels
(U-Ones), or converted to $negative$ labels (U-Zeros), proposed in (Irvin et
al. 2019). Moreover, in (Pham et al. 2019), label smoothing regularization
(LSR) was applied to better handle these uncertain labels, with the main goal
of preventing the model from training with overconfident labels that might
contain misclassified data. Specifically, in the U-Ones setting, the
uncertainty labels are mapped to a random number close to one, subject to a
uniform distribution $\mathcal{U}(a,b)$, where $a$ and $b$ are the hyper-
parameters. In contrast, in the U-Zeros setting, a similar process is used to
convert the uncertainty labels to a random number close to zero. Although
these solutions provide some improvements, they are based on too much on
brute-force to address the mislabeled data reasonably.

### Many-to-One Distribution Learning

Incorporating the opinions of multiple radiologists is usually more reliable
than only adopting a single diagnosis. Similarly, although the target network
$f_{T}$ based on the single backbone DenseNet-121 can achieve satisfactory
results, model ensemble of several network architectures usually further makes
a big improvement on the mean accuracy. According to many observations from
extensive experiments, different network architectures, and even different
epoch checkpoints for one network architecture, provide inconsistent
predictions and have their own strengths on different diseases. Thus,
integrating probability predictions from multiple models can reduce the
effects of the uncertainty and ambiguity from the original binary labels. In
this paper, we propose a many-to-one distribution learning method to optimize
a single target network by taking the essence and discarding the dregs of
multiple networks, which can also avoid the extra computational resources
introduced by a model ensemble in the testing phase.

To select a strong set of reference models, we follow (Pham et al. 2019) and
adopt six state-of-the-art CNNs ($N=6$), including DenseNet-121, DenseNet-169,
DenseNet-201 (Huang et al. 2017), Inception-ResNet-v2 (Szegedy et al. 2017),
Xception (Chollet 2017), and NASNetLarge (Zoph et al. 2018). Once the training
is separately done for each network, the probability predictions of all the
trained networks are simply averaged, and can then be viewed as the soft label
distribution. By getting an average, those evident lesions with high
probabilities predicted by most of the reference networks still get positive-
toward values, and healthy scans still get negative-toward values. However,
those lesions getting inconsistent predictions which have ambiguities will be
smoothed to reduce the effects of uncertainty. Our label distribution
smoothing tends to be a major voting manner, which is more reasonable than the
random LSR method. Therefore, this soft label distribution for each image
should be closer to the correct labels to indicate the ground truths. We use
the Kullback-Leibler (KL) divergence for the distribution learning loss
$\mathcal{L}_{distri}$ and minimize it to optimize the target network $f_{T}$:

$\displaystyle\mathcal{L}_{distri}=KL\\{l_{i}||f_{T}(x_{i}|\theta)\\}=\sum_{i}l_{i}\log(\frac{l_{i}}{f_{T}(x_{i}|\theta)})\text{.}$
(2)

Please note that all the reference networks $f_{n}$ are pre-trained in
advance. Their network parameters are fixed when training the target network
$f_{T}$. Thus, the soft label distribution $l_{i}$ is fixed to compute the KL
divergence as the relative entropy of $l_{i}$ with respect to
$f_{T}(x_{i}|\theta)$.

###  $K$-Nearest Neighbor Smoothing

Our second important concern is how to make the target network as robust and
competitive as the model ensemble. According to the smoothness assumption
(Chapelle, Scholkopf, and Zien 2009), we assume that if two samples are
similar in the averaged label distribution space, they should also be similar
in the label space predicted by the target model. Therefore, before training
the target network, the $K$ nearest neighbors of an image (marked as the
anchor image) are pre-computed based on the fixed soft label distribution.

To make the disease predictions of the anchor image $f_{T}(x_{i}|\theta)$ as
close as the predictions of its corresponding $k$-th neighbor
$f_{T}(x_{i}^{k}|\theta)$, we generate a pool that contains $K$ neighbor
images of each anchor during the training, for enhancing local similarity. Let
$k$ denote a sample in the pool. Inspired by (Xu, Tao, and Geng 2018), we
specify the local similarity $s_{i}^{k}$ which is defined as:

$\displaystyle
s_{i}^{k}=\begin{cases}\mathrm{exp}(-\frac{||l_{i}-l_{i}^{k}||^{2}}{2\sigma^{2}})&\mbox{if
$k\in\\{1,...,K\\}$}\\\ 0&\mbox{otherwise}\text{,}\end{cases}$ (3)

where $\sigma$ denotes the width parameter for similarity calculation, which
is fixed as 1 in our settings. Thus, the more similar $l_{i}$ and $l_{i}^{k}$
are, the higher $s_{i}^{k}$ will be, which means $f_{T}(x_{i}|\theta)$ and
$f_{T}(x_{i}^{k}|\theta)$ should be closer. We propose a neighbor smoothing
loss $\mathcal{L}_{neigh}$ as follows:

$\displaystyle\mathcal{L}_{neigh}=\sum_{i,k}s_{i}^{k}f_{T}(x_{i}^{k}|\theta)\log\frac{f_{T}(x_{i}^{k}|\theta)}{f_{T}(x_{i}|\theta)}\text{.}$
(4)

By optimizing $\mathcal{L}_{neigh}$, the target model $f_{T}$ is constrained
to give similar predictions of neighboring images which have similar abnormal
findings. This enhances the robustness of $f_{T}$ for more consistent
identification performance.

Methods | Atel. | Card. | Effu. | Infi. | Mass | Nod. | Pne.1 | Pne.2 | Cos. | Ede. | Emp. | Fibr. | P.T. | Hern | Mean
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
(Wang et al. 2017) | 70.03 | 81.00 | 75.85 | 66.14 | 69.33 | 66.87 | 65.80 | 79.93 | 70.32 | 80.52 | 83.30 | 78.59 | 68.35 | 87.17 | 74.51
(Yao et al. 2018) | 73.30 | 85.80 | 80.60 | 67.50 | 72.70 | 77.80 | 69.00 | 80.50 | 71.70 | 80.60 | 84.20 | 75.70 | 72.40 | 82.40 | 76.73
(Yao et al. 2017) | 77.20 | 90.40 | 85.90 | 69.50 | 79.20 | 71.70 | 71.30 | 84.10 | 78.80 | 88.20 | 82.90 | 76.70 | 76.50 | 91.40 | 80.30
(Li et al. 2018b) | 80.00 | 81.00 | 87.00 | 70.00 | 83.00 | 75.00 | 67.00 | 87.00 | 80.00 | 88.00 | 91.00 | 78.00 | 79.00 | 77.00 | 80.21
(Rajpurkar et al. 2017) | 80.94 | 92.48 | 86.38 | 73.45 | 86.76 | 78.02 | 76.80 | 88.87 | 79.01 | 88.78 | 93.71 | 80.47 | 80.62 | 91.64 | 84.13
Baseline (B) | 82.42 | 89.80 | 89.68 | 68.46 | 86.84 | 80.34 | 75.31 | 85.34 | 80.37 | 91.62 | 90.82 | 84.60 | 80.64 | 93.97 | 84.30
B+MODL | 83.91 | 91.66 | 91.58 | 72.84 | 88.68 | 81.65 | 76.73 | 89.32 | 82.19 | 93.72 | 93.83 | 86.51 | 82.92 | 95.69 | 86.52
B+KNNS | 84.44 | 91.73 | 91.62 | 72.73 | 88.99 | 82.12 | 76.97 | 89.52 | 82.11 | 93.75 | 93.81 | 86.78 | 83.74 | 95.77 | 86.72
B+MODL+KNNS | 85.20 | 92.24 | 92.09 | 73.01 | 89.74 | 82.94 | 78.30 | 90.46 | 83.04 | 94.42 | 94.32 | 87.64 | 84.26 | 95.28 | 87.35

Table 1: Identification results using AUC (%) metric on the NIH Chest X-ray
dataset.

Methods | Atelectasis | Cardiomegaly | Consolidation | Edema | Pleural Effusion | Mean
---|---|---|---|---|---|---
(Allaouzi and Ahmed 2019) LP | 72.0 | 87.0 | 77.0 | 87.0 | 90.0 | 82.6
(Allaouzi and Ahmed 2019) BR | 72.0 | 88.0 | 77.0 | 87.0 | 90.0 | 82.8
(Allaouzi and Ahmed 2019) CC | 70.0 | 87.0 | 74.0 | 86.0 | 90.0 | 81.4
(Irvin et al. 2019) U-Zeros * | 81.1 | 84.0 | 93.2 | 92.9 | 93.1 | 88.8
(Irvin et al. 2019) U-Ones * | 85.8 | 83.2 | 89.9 | 94.1 | 93.4 | 89.3
(Pham et al. 2019) U-Zeros+CT+LSR | 80.6 | 83.3 | 92.9 | 93.3 | 92.1 | 88.4
(Pham et al. 2019) U-Ones+CT+LSR | 82.5 | 85.5 | 93.7 | 93.0 | 92.3 | 89.4
Baseline (B) | 87.29 | 84.47 | 92.09 | 91.66 | 91.34 | 89.37
B+MODL | 87.64 | 87.84 | 93.34 | 94.58 | 93.35 | 91.35
B+KNNS | 89.15 | 88.24 | 92.40 | 94.92 | 92.67 | 91.48
B+MODL+KNNS | 89.46 | 89.67 | 93.80 | 95.87 | 94.35 | 92.63

Table 2: Identification results using AUC (%) metric on the CheXpert dataset.
* indicates that the result is obtained by the ensemble of 30 checkpoints.

### Implementation Details

The training procedure of our proposed method consists of two steps. In the
first step, all the reference models based on different backbones are pre-
trained using the basic multi-label classification loss $\mathcal{L}_{cls}$.
All the images are resized into a standard size $256\times 256$ and randomly
cropped to $224\times 224$ patches as inputs. The Adam optimizer is adopted
with an initial learning rate of 0.001 and default parameters
$\beta_{1}=0.9,\beta_{2}=0.999$. The network parameters are initialized with a
model pre-trained on ImageNet (Deng et al. 2009). The mini-batch size is set
to 32 for training over 10 epochs and the learning rate is reduced by a factor
of 3 after every two epochs. Once all the reference models are well trained,
the soft label distribution of each image is computed and used for training in
the second step.

In the second step of the training phase, in addition to $\mathcal{L}_{cls}$,
we add the proposed MODL loss $\mathcal{L}_{distri}$ and KNNS loss
$\mathcal{L}_{neigh}$ for optimizing the single target model. For each image
in a training batch (we marked it as an anchor), its $K$ neighbor images and
the corresponding local similarity scores $s_{i}^{k}$ are stored in a pool for
optimizing $\mathcal{L}_{neigh}$. The different loss weights $\lambda$ and
$\gamma$ in Eq. 1 are both set to 0.1 in the experiments. All other training
settings in the second step are the same as Step 1. Once the training is
completed, only the target model $f_{T}$ is deployed in the testing phase,
without any additional computation consumption introduced by model ensemble.

## Experiments and Results

### Datasets

#### NIH Chest X-ray (Wang et al. 2017)

is comprised of 112,120 X-ray images with disease labels from 30,805 unique
patients, 51,708 of which contain one or more pathologies, while the remaining
images have no findings. The 14 different labels are extracted using NLP to
text-mine disease classifications from the associated radiology reports, but
more than 10% of the annotations are incorrect. The dataset split setting of
experiments for our model follows (Rajpurkar et al. 2017).

#### CheXpert (Irvin et al. 2019)

consists of 224,316 chest radiographs of 65,240 patients. The disease labels
of the training set are labeled for the presence of 14 observations as
positive, negative, uncertain, or unmentioned, which also contain
inconsistency and ambiguity. However, CheXpert releases a separate validation
set of 200 studies, which is annotated by three board-certified radiologists
for five diseases: atelectasis, cardiomegaly, consolidation, edema, and
pleural effusion. The evaluation on this set is more convincing.

Figure 3: ROC curves of different methods for the five pathologies on the
CheXpert dataset.

### Comparison with State-of-the-Arts

To evaluate the multi-disease identification task, the area-under-the-curve
(AUC) of the receiver operating characteristic (ROC) is usually adopted. Table
1 provides the identification results on the NIH Chest X-ray dataset. Model
(Wang et al. 2017), using a ResNet-50 backbone and weighted binary cross-
entropy loss, obtains a preliminary mean AUC of 74.51% over 14 diseases.
CheXnet (Rajpurkar et al. 2017) adopts DenseNet-121 and obtains a large
improvement, since DenseNets improve the flow of information and gradients
through the network, making the optimization of deep networks tractable. In
(Yao et al. 2017), a long short-term memory (LSTM) architecture is used to
model dependencies among labels, but achieves poor results. For most diseases,
our method with many-to-one distribution learning and $K$-nearest neighbor
smoothing obtains consistent increases in identification performance, compared
to the best-performing method (Rajpurkar et al. 2017). The only exception is
on cardiomegaly and infiltration. The mean AUC is increased by 3.22%.

The results on CheXpert are compared in Table 2. Most methods exploit
ensembles of multiple models to pursue higher scores. For example, an ensemble
of 30 generated checkpoints of DenseNet-121 (Irvin et al. 2019) achieves a
mean AUC of 89.3%. In this subsection, we mainly focus on comparing the
results achieved by single models. Binary relevance (BR), label powerset (LP),
and classifier chain (CC) are explored in terms of label dependencies in
(Allaouzi and Ahmed 2019), but no satisfactory performance is achieved.
Conditional training (CT) and label smoothing regularization (LSR) are used to
enhance the performance of a single DenseNet with a small gain achieved. In
our implementation, the U-Ones setting is adopted to deal with those
uncertainty labels when training the Baseline (B) (based on the backbone
DenseNet-121), as well as when training the reference models in the first
training stage. Our method B+MODL+KNNS without any ensemble operation
increases the mean AUC by 3.26%.

### Ablation Studies

#### Effectiveness of MODL

To evaluate the effectiveness of the proposed many-to-one distribution
learning, the soft label distribution obtained is only used to optimize the
target network through the auxiliary loss $\mathcal{L}_{distri}$. Clear
improvements in mean AUC 2.22% and 1.98% are achieved on NIH Chest X-ray and
CheXpert, respectively, compared to the baselines. This illustrates that the
soft label distribution can represent the ground truths of a scan better than
the original automatically extracted binary labels which have uncertainty and
ambiguity.

Figure 4: Comparison of the mean AUC increase based on different $K$ in KNNS.

#### Effectiveness of KNNS

The $K$-nearest neighbor smoothing is also separately validated. An anchor
image and its $K$ neighbor images are passed forward to the target network,
which is updated by the neighbor smoothing loss $\mathcal{L}_{neigh}$. The
results show that constraining the model to provide consistent predictions for
images with similar medical findings also largely improves the performance.
The mean AUC is increased by 2.42% and 2.11% on NIH Chest X-ray and CheXpert,
respectively. Moreover, combining $\mathcal{L}_{distri}$ and
$\mathcal{L}_{neigh}$ can further enhance the target model without requiring
any extra computational resources during the testing phase. Figure 3 provides
the ROC curves of different methods for the five pathologies on the CheXpert
validation set.

#### Investigation of the Number $K$ of KNNS

Another concern is to determine the best-performing $K$ of KNNS. We conduct an
experiment on the CheXpert dataset, where $K$ is increased from 3 to 11 by a
step of 2. The mean AUC improvement in terms of $K$ is shown in Figure 4. We
observe that configuring three neighbor images to constrain the model
prediction consistency can already obtain a significant increase in mean AUC
of 1.14%. Improvements of 1.68% and 1.98% are achieved when we increase $K$ to
5 and 7, respectively. However, with the further growth of $K$, the increase
rate becomes smaller and the method requires much higher computational memory
during the training phase. Therefore, considering the trade-off between
performance and computational costs, we set $K$ as 9 in KNNS.

Methods | Mean AUC | Params.
---|---|---
DenseNet-121 (B) | 89.66 | 7.98 M
B+MODL+KNNS | 92.63 | 7.98 M
Model Ensemble 1 | 92.19 | 42.13 M
Model Ensemble 2 | 93.96 | 208.33 M

Table 3: Comparison with model ensembles on the CheXpert dataset. M denotes
million. Figure 5: Heatmap visualization on the CheXpert dataset. Each row
shows two cases of one disease. For each case, the left is the input image
with ground truths. The middle and right are the heatmaps obtained by the
Baseline (B) and B+MODL+KNNS, respectively. The value on the top-left of each
heatmap is the predicted probability of the corresponding disease.

#### Comparison with Model Ensembles

In addition to the mitigation of label uncertainty and enhancement of model
robustness, another big advantage of our method is the zero-increase in
network complexity during the testing phase, compared to model ensemble
methods. To validate this point, the mean AUC and number of network parameters
of some ensemble models are compared in Table 3. Model Ensemble 1 is a
combination of DenseNet-121, DenseNet-169, and DenseNet-201. Model Ensemble 2
is an ensemble of the six reference models used for learning the soft label
distribution in our method. As shown from the results, our method B+MODL+KNNS
outperforms Model Ensemble 1 by 0.44% of the mean AUC, with much fewer
parameters. Moreover, although Model Ensemble 2 obtains better identification
performance than ours, it has 26 times more parameters than our target model,
making it too heavy to deploy in practical applications.

### Qualitative Results

To better visualize the model performance, we adopt the probabilistic class
activation map (PCAM) (Ye et al. 2020) pooling to predict the lesion heatmaps
with only image-level supervision. As shown in Figure 5, we visualize the
heatmaps of the five pathologies in the CheXpert dataset. All the examples do
have the corresponding pathology and are labeled positive. For each case, we
compare the heatmaps obtained by the Baseline (B) and B+MODL+KNNS. In addition
to the improvements in identification results demonstrated before, we observe
that a more satisfactory localization performance can also be achieved. For
all the pathologies, the lesion regions predicted by B+MODL+KNNS are more
precise in location and size than those predicted by the Baseline (B). For
example, the cardiomegaly heatmaps predicted by the Baseline (B) usually have
an abnormal region around the location of the mediastinum, while B+MODL+KNNS
can better highlight the enlarged heart shape region. To localize the position
of pleural effusion, B+MODL+KNNS can better identify the left heart border,
costophrenic angle and hemidiaphragm are obscured, and slight blunting of the
right costophrenic angle. Moreover, as shown from the predicted probabilities
in the top-left of each heatmap, B+MODL+KNNS identifies the positive diseases
as true with higher confidence scores, compared to the Baseline (B).

## Conclusion

In this paper, we proposed a many-to-one distribution learning and a
$K$-nearest neighbor smoothing method to improve the performance of thoracic
disease diagnosis. The methods can reduce the effects of label uncertainty and
ambiguity, and also constrain the target model to provide consistent
predictions on images with similar medical findings. Extensive experiments
demonstrated that our method is effective and achieves significant
improvements without requiring any additional computational costs in the
testing phase. In our future work, we will explore visual attention methods to
learn better representations of thoracic diseases for improvements.

## References

  * Allaouzi and Ahmed (2019) Allaouzi, I.; and Ahmed, M. B. 2019.  A novel approach for multi-label chest X-ray classification of common thorax diseases.  _IEEE Access_ 7: 64279–64288.
  * Baltruschat et al. (2019) Baltruschat, I. M.; Nickisch, H.; Grass, M.; Knopp, T.; and Saalbach, A. 2019.  Comparison of deep learning approaches for multi-label chest X-ray classification.  _Scientific reports_ 9(1): 1–10.
  * Chapelle, Scholkopf, and Zien (2009) Chapelle, O.; Scholkopf, B.; and Zien, A. 2009.  Semi-supervised learning (chapelle, o. et al., eds.; 2006)[book reviews].  _IEEE Transactions on Neural Networks_ 20(3): 542–542.
  * Chen et al. (2020) Chen, S.; Wang, J.; Chen, Y.; Shi, Z.; Geng, X.; and Rui, Y. 2020.  Label Distribution Learning on Auxiliary Label Space Graphs for Facial Expression Recognition.  In _CVPR_ , 13984–13993. IEEE.
  * Chollet (2017) Chollet, F. 2017.  Xception: Deep learning with depthwise separable convolutions.  In _CVPR_ , 1251–1258. IEEE.
  * Deng et al. (2009) Deng, J.; Dong, W.; Socher, R.; Li, L.-J.; Li, K.; and Fei-Fei, L. 2009.  Imagenet: A large-scale hierarchical image database.  In _CVPR_ , 248–255. IEEE.
  * Fan et al. (2020) Fan, D.-P.; Zhou, T.; Ji, G.-P.; Zhou, Y.; Chen, G.; Fu, H.; Shen, J.; and Shao, L. 2020.  Inf-net: Automatic covid-19 lung infection segmentation from ct images.  _IEEE Transactions on Medical Imaging_ 39(8): 2626–2637.
  * Gao et al. (2017) Gao, B.-B.; Xing, C.; Xie, C.-W.; Wu, J.; and Geng, X. 2017.  Deep label distribution learning with label ambiguity.  _IEEE Transactions on Image Processing_ 26(6): 2825–2838.
  * Geng (2016) Geng, X. 2016.  Label distribution learning.  _IEEE Transactions on Knowledge and Data Engineering_ 28(7): 1734–1748.
  * Hou et al. (2017) Hou, P.; Geng, X.; Huo, Z.-W.; and Lv, J.-Q. 2017.  Semi-supervised adaptive label distribution learning for facial age estimation.  In _AAAI_.
  * Huang et al. (2017) Huang, G.; Liu, Z.; Van Der Maaten, L.; and Weinberger, K. Q. 2017.  Densely connected convolutional networks.  In _CVPR_ , 4700–4708. IEEE.
  * Irvin et al. (2019) Irvin, J.; Rajpurkar, P.; Ko, M.; Yu, Y.; Ciurea-Ilcus, S.; Chute, C.; Marklund, H.; Haghgoo, B.; Ball, R.; Shpanskaya, K.; et al. 2019.  Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison.  In _AAAI_ , 590–597.
  * Jing, Xie, and Xing (2018) Jing, B.; Xie, P.; and Xing, E. 2018.  On the automatic generation of medical imaging reports.  In _ACL_ , 2577–2586.
  * Johnson et al. (2019) Johnson, A. E.; Pollard, T. J.; Berkowitz, S.; Greenbaum, N. R.; Lungren, M. P.; Deng, C.-y.; Mark, R. G.; and Horng, S. 2019.  MIMIC-CXR: A large publicly available database of labeled chest radiographs.  _arXiv preprint arXiv:1901.07042_ .
  * Kumar, Grewal, and Srivastava (2018) Kumar, P.; Grewal, M.; and Srivastava, M. M. 2018.  Boosted cascaded convnets for multilabel classification of thoracic diseases in chest radiographs.  In _International Conference Image Analysis and Recognition_ , 546–552. Springer.
  * Li et al. (2018a) Li, Y.; Liang, X.; Hu, Z.; and Xing, E. P. 2018a.  Hybrid retrieval-generation reinforced agent for medical image report generation.  In _NeurIPS_ , 1530–1540.
  * Li et al. (2018b) Li, Z.; Wang, C.; Han, M.; Xue, Y.; Wei, W.; Li, L.-J.; and Fei-Fei, L. 2018b.  Thoracic disease identification and localization with limited supervision.  In _CVPR_ , 8290–8299. IEEE.
  * Litjens et al. (2017) Litjens, G.; Kooi, T.; Bejnordi, B. E.; Setio, A. A. A.; Ciompi, F.; Ghafoorian, M.; Van Der Laak, J. A.; Van Ginneken, B.; and Sánchez, C. I. 2017\.  A survey on deep learning in medical image analysis.  _Medical Image Analysis_ 42: 60–88.
  * Liu et al. (2019) Liu, J.; Zhao, G.; Fei, Y.; Zhang, M.; Wang, Y.; and Yu, Y. 2019.  Align, Attend and Locate: Chest X-ray Diagnosis via Contrast Induced Attention Network with Limited Supervision.  In _ICCV_ , 10632–10641. IEEE.
  * Northcutt, Jiang, and Chuang (2019) Northcutt, C. G.; Jiang, L.; and Chuang, I. L. 2019.  Confident learning: Estimating uncertainty in dataset labels.  _arXiv preprint arXiv:1911.00068_ .
  * Pesce et al. (2019) Pesce, E.; Withey, S. J.; Ypsilantis, P.-P.; Bakewell, R.; Goh, V.; and Montana, G. 2019.  Learning to detect chest radiographs containing pulmonary lesions using visual attention networks.  _Medical Image Analysis_ 53: 26–38.
  * Pham et al. (2019) Pham, H. H.; Le, T. T.; Tran, D. Q.; Ngo, D. T.; and Nguyen, H. Q. 2019.  Interpreting chest X-rays via CNNs that exploit disease dependencies and uncertainty labels.  _arXiv preprint arXiv:1911.06475_ .
  * Rajpurkar et al. (2017) Rajpurkar, P.; Irvin, J.; Zhu, K.; Yang, B.; Mehta, H.; Duan, T.; Ding, D.; Bagul, A.; Langlotz, C.; Shpanskaya, K.; et al. 2017.  Chexnet: Radiologist-level pneumonia detection on chest x-rays with deep learning.  _arXiv preprint arXiv:1711.05225_ .
  * Shin et al. (2016) Shin, H.-C.; Roth, H. R.; Gao, M.; Lu, L.; Xu, Z.; Nogues, I.; Yao, J.; Mollura, D.; and Summers, R. M. 2016.  Deep convolutional neural networks for computer-aided detection: CNN architectures, dataset characteristics and transfer learning.  _IEEE Transactions on Medical Imaging_ 35(5): 1285–1298.
  * Su and Geng (2019) Su, K.; and Geng, X. 2019.  Soft facial landmark detection by label distribution learning.  In _AAAI_ , volume 33, 5008–5015.
  * Szegedy et al. (2017) Szegedy, C.; Ioffe, S.; Vanhoucke, V.; and Alemi, A. A. 2017.  Inception-v4, inception-resnet and the impact of residual connections on learning.  In _AAAI_.
  * Wang et al. (2017) Wang, X.; Peng, Y.; Lu, L.; Lu, Z.; Bagheri, M.; and Summers, R. M. 2017.  Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases.  In _CVPR_ , 2097–2106. IEEE.
  * Wang et al. (2018) Wang, X.; Peng, Y.; Lu, L.; Lu, Z.; and Summers, R. M. 2018.  Tienet: Text-image embedding network for common thorax disease classification and reporting in chest x-rays.  In _CVPR_ , 9049–9058. IEEE.
  * Wu et al. (2018) Wu, J.; Ruan, S.; Lian, C.; Mutic, S.; Anastasio, M. A.; and Li, H. 2018.  Active learning with noise modeling for medical image annotation.  In _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_ , 298–301. IEEE.
  * Xu, Lv, and Geng (2019) Xu, N.; Lv, J.; and Geng, X. 2019.  Partial label learning via label enhancement.  In _AAAI_ , volume 33, 5557–5564.
  * Xu, Tao, and Geng (2018) Xu, N.; Tao, A.; and Geng, X. 2018.  Label Enhancement for Label Distribution Learning.  In _IJCAI_ , 2926–2932.
  * Yao et al. (2017) Yao, L.; Poblenz, E.; Dagunts, D.; Covington, B.; Bernard, D.; and Lyman, K. 2017\.  Learning to diagnose from scratch by exploiting dependencies among labels.  _arXiv preprint arXiv:1710.10501_ .
  * Yao et al. (2018) Yao, L.; Prosky, J.; Poblenz, E.; Covington, B.; and Lyman, K. 2018.  Weakly supervised medical diagnosis and localization from multiple resolutions.  _arXiv preprint arXiv:1803.07703_ .
  * Ye et al. (2020) Ye, W.; Yao, J.; Xue, H.; and Li, Y. 2020.  Weakly Supervised Lesion Localization With Probabilistic-CAM Pooling.  _arXiv preprint arXiv:2005.14480_ .
  * Zhang et al. (2020) Zhang, Y.; Wang, X.; Xu, Z.; Yu, Q.; Yuille, A.; and Xu, D. 2020.  When Radiology Report Generation Meets Knowledge Graph.  In _AAAI_ , 12910–12917.
  * Zhou, Li, and Wang (2018) Zhou, B.; Li, Y.; and Wang, J. 2018.  A weakly supervised adaptive densenet for classifying thoracic diseases and identifying abnormalities.  _arXiv preprint arXiv:1807.01257_ .
  * Zhou et al. (2021) Zhou, Y.; Zhou, T.; Zhou, T.; Fu, H.; Liu, J.; and Shao, L. 2021.  Contrast-Attentive Thoracic Disease Recognition with Dual-Weighting Graph Reasoning.  _IEEE Transactions on Medical Imaging_ 1–1.  doi:10.1109/TMI.2021.3049498.
  * Zoph et al. (2018) Zoph, B.; Vasudevan, V.; Shlens, J.; and Le, Q. V. 2018.  Learning transferable architectures for scalable image recognition.  In _CVPR_ , 8697–8710. IEEE.