\section{Introduction}
\label{sec:intro}

The dynamics of hadrons with different flavor degrees of freedom provides different windows for
our understanding of the underlying theory of strong interaction, quantum chromodynamics (QCD).
With regard to hadrons with charm, so far spectroscopy has been the most visible
and definitely by far the most interesting branch of research. Indeed, the
large number of structures observed in experiments at energies above the open charm
production threshold provides a challenge for our standard (but obviously naive) picture
that mesons are composed out of quark-antiquark pairs and baryons out
of three quarks. See Refs.~\cite{Guo:2017,Lebed:2016hpi,Chen:2016,Esposito:2016}
for recent overviews and discussions of these structures, commonly
referred to as X, Y and Z states.

Some proposals for experiments at sites such as J-PARC \cite{JPARC}
and FAIR \cite{CBM,Wiedner:2011} aim at exploring also other aspects of charm physics.
Specifically, so-called charm factories would allow
the targeted production of charmed hadrons like the $D$-meson or the ${\Lambda}_c$ and
${\Sigma}_c$ hyperons and a study of their interaction with ordinary hadrons.
The expectation of possible experiments in the not so far future has triggered a variety
of theoretical investigations. In particular, it has led to a renewed interest in the
interaction of the ${\Lambda}_c$ with nucleons and with nuclei over the last few years
\cite{Liu:2011xc,Huang:2013,Gal:2014,Garcilazo:2015qha,Maeda:2015hxa,Shyam:2016uxa,Ohtani:2017wdc}.
Those studies join the ranks of a long history of speculations about bound nuclear systems
involving the ${\Lambda}_c$, the lightest charmed
baryon~\cite{Dover:1977jw,Bando:1981ti,Bando:1983yt,Gibson:1983zw,Bunyatov:1992in,Tsushima:2002ua,Tsushima:2003dd,Kopeliovich:2007kd} {\textemdash} see also the recent reviews~\cite{Hosaka:2016ypm,Krein:2017usp}.
Indeed, in most of the investigations so far, the ${Y}_c N$ interaction
(${Y}_c = \Lambda_c, \Sigma_c$), derived within the meson-exchange
framework~\cite{Liu:2011xc,Dover:1977jw,Bando:1981ti,Bando:1983yt,Gibson:1983zw} or in
the constituent quark model~\cite{Maeda:2015hxa,Froemel:2004ea}, turns out to be
strongly attractive.

Interestingly, quite the opposite picture emerged from recent (2+1)-flavor lattice QCD (LQCD)
simulations by the HAL QCD Collaboration~\cite{Miyamoto:2016,Miyamoto:2016A,Miyamoto:2017}.
Pertinent calculations, performed for unphysical quark masses corresponding to pion masses\footnote{The
Gell-Mann-Oakes-Renner relation states that the squared pion mass is proportional to the average
light quark mass. Therefore, the notions ``quark mass'' and ``pion mass'' are used synonymously.}
of $m_\pi = 410 - 700$ MeV,
suggest that the ${\Lambda}_c N$ and ${\Sigma}_c N$ interactions could be much less attractive
than predicted by the phenomenological potentials mentioned above.
While initial preliminary studies \cite{Miyamoto:2016,Miyamoto:2016A} indicated
an extremely weak ${\Lambda}_c N$ interaction, the recently published final LQCD results
\cite{Miyamoto:2017} rectify that conjecture and imply a somewhat stronger,
though still only moderately attractive ${\Lambda}_c N$ interaction\footnote{We note that the
method employed in these works and also other methods are
presently under discussion in the LQCD
community~\cite{Aoki:2017,Yamazaki:2017,Davoudi:2017,Aoki:2017A}.}.

In the present work we provide predictions for the ${\Lambda}_c N$ interaction at the physical point
based on the LQCD simulations by the HAL QCD Collaboration \cite{Miyamoto:2017}.
The extrapolation of the LQCD results, available for $m_\pi = 410 - 700$ MeV,
to $m_\pi = 138$ MeV is performed within the framework of chiral effective field
theory (EFT)~\cite{Epelbaum:2005,Epelbaum:2008ga,Machleidt:2011zz}.
Thereby we follow a strategy that has been already employed by us in the past in the
analysis of other LQCD results on baryon-baryon interactions in the strangeness sector,
notably the $\Lambda\Lambda$~\cite{Haidenbauer:2011a}, $\Xi\Xi$ \cite{Haidenbauer:2011b},
and $\Omega\Omega$ \cite{Haidenbauer:2017} systems:
At first chiral EFT is utilized to establish a ${\Lambda}_c N$ potential for baryon and meson
masses that correspond to those in the lattice simulation. In particular, open parameters
are determined by a fit to pertinent LQCD results (phase shifts, scattering lengths).
Then the potential is extrapolated to the physical point. Thereby the pion-mass
dependence of the ingredients (pseudoscalar-meson exchange, four-baryon
contact terms) is taken into account explicitly, and in line with chiral EFT.

The paper is structured as follows: In Sect.~2 we provide an outline of the
employed formalism.
Results for ${\Lambda}_c N$ phase shifts (for the $^1\!S_0$ and $^3\!S_1$ partial waves)
are reported in Sect.~3, for pion masses corresponding to those in the
LQCD simulation and for the physical value, $m_\pi = 138$ MeV.
The corresponding scattering lengths are evaluated too and turn out to be in
the order of $a\approx -1$~fm. Finally, consequences of our results for the
existence of bound ${\Lambda}_c$ hypernuclei are discussed.
The paper closes with a short summary.

\section{Formalism}
\label{sec:formalism}

The $Y_c N$ interaction is constructed by using chiral EFT as guideline.
Thereby we follow closely our application of this scheme to the $\Lambda N$ and
$\Sigma N$ systems in Refs.~\cite{Polinder,Haidenbauer:2013} where corresponding
potentials have been obtained up to NLO in the Weinberg power counting
\cite{Epelbaum:2008ga,Machleidt:2011zz}.
In this framework the potential is given in terms of pion exchanges
and a series of contact interactions with an increasing number of derivatives.
The latter represent the short-range part of the baryon-baryon force and are
parameterized by low-energy constants (LECs), that need to be fixed in a fit
to data.
Both classes of contributions depend on the quark mass (or, equivalently, the pion mass).
At LO the only quark-mass dependence of the ${Y}_c N$ potential is through the
pion mass that appears in the propagator of the pion-exchange potential, cf. below.
However, at NLO the contact terms as well as the pion coupling constants depend on the
pion mass. For details, we refer to
Refs.~\cite{Beane:2002vs,Beane03,Epe02,Epe02a,Baru:2015,Baru:2016},
where the quark mass dependence of the $NN$ interaction has been investigated;
see also Ref.~\cite{Petschauer:2013}.

Let us start by introducing the contact interaction that we employ. For the partial waves
considered in the present study ($^1\!S_0$, $^3\!S_1$-$^3\!D_1$) it is given by
\begin{eqnarray}
V(^1\!S_0) &=& {\tilde{C}}_{^1\!S_0} + \tilde{D}_{^1\!S_0} m^2_\pi \nonumber \\
&& + \, ({C}_{^1\!S_0}
+ {D}_{^1\!S_0} m^2_\pi)\, ({p}^2+{p}'^2) \ , \nonumber \\
V(^3\!S_1) &=& {\tilde{C}}_{^3\!S_1} + \tilde{D}_{^3\!S_1} m^2_\pi \nonumber \\
&& + \, ({C}_{^3\!S_1} \,
+ {D}_{^3\!S_1} m^2_\pi)\, ({p}^2+{p}'^2) \ , \nonumber \\
V(^3\!D_1 -\, ^3\!S_1) &=& {C}_{\varepsilon_1}\, {p'}^2 \ , \nonumber \\
V(^3\!S_1 -\, ^3\!D_1) &=& {C}_{\varepsilon_1}\, {p}^2 \ ,
\label{LEC}
\end{eqnarray}
with $p = |{\bf p}\,|$ and ${p}' = |{\bf p}\,'|$  being the initial and final
center-of-mass (c.m.) momenta in the ${\Lambda}_c N$ or ${\Sigma}_c N$ systems.
The quantities ${\tilde{C}}_{i}$, ${\tilde{D}}_{i}$, ${C}_{i}$, ${D}_{i}$ are the
aforementioned LECs that need to be fixed by a fit to lattice data (phase shifts).
The ansatz Eq.~(\ref{LEC}) is motivated by the corresponding expression in the standard Weinberg
counting up to NLO \cite{Epe02,Petschauer:2013} but differs from it by the terms proportional to
$m^2_\pi\, ({p}^2+{p}'^2)$ which are nominally of higher order. Nevertheless, we include these
terms because they allow us to obtain an optimal description of the LQCD results at $m_\pi=410$ MeV
as well as at $570$ MeV and, thereby, enable us a better constrained extrapolation to lower pion
masses. Contrary to our study of the $YN$ interaction in Refs.~\cite{Polinder,Haidenbauer:2013},
here we do not impose SU(3) (or SU(4)) flavor symmetry.
In any case, given that there is no explicit information on the ${\Sigma}_c N$ channel from LQCD in
Ref.~\cite{Miyamoto:2017} and the ${\Lambda}_c N$ interaction is determined as an effective
single-channel potential, we treat the ${\Lambda}_c N$ interaction likewise as an effective
single-channel problem.

The contribution of pion exchange to the ${Y}_c N$ potential is given by
\begin{equation}
V^{OPE}_{{B N\to B' N}} =-f_{{B B'\pi}}f_{{NN\pi}}
\frac{\left({\bf \sigma}_1\cdot {\bf q} \right)
\left({\bf \sigma}_2\cdot {\bf q} \right)}{ {\bf q}^{\,2}+m_{\pi}^2}
\mathcal{I}_{{B N\to B' N}}\ ,
\label{OPE}
\end{equation}
where $B,B'$ stand for ${\Lambda}_c$ and/or ${\Sigma}_c$,
${\bf q}$ is the transferred momentum, ${\bf q} = {\bf p'} - {\bf p}$,
and $\cal{I}$ is a pertinent isospin factor \cite{Polinder}.
As already mentioned, we do not assume the validity of SU(4) flavor symmetry in the
present study. This concerns also the coupling constants $f_{{B B'\pi}}$.
The ${\Lambda}_c{\Sigma}_c\pi$ coupling constant can be determined from the experimentally
known ${\Sigma}_c \to {\Lambda}_c\pi$ decay rate, see Refs.~\cite{Albertus:2005,Can:2016}.
With regard to the ${\Sigma}_c\Sc\pi$ coupling constant we resort to LQCD results~\cite{Alexandrou:2016}.
Besides their value at the physical point, we need also
the $m_\pi$ dependence of the ${\Lambda}_c{\Sigma}_c\pi$ and ${\Sigma}_c\Sc\pi$
coupling constants as well as the one for the $NN\pi$ vertex,
\begin{equation}
f_{BB'\pi} (m^2_\pi) = \frac{g_A^{BB'}(m^2_\pi)}{2\,F_\pi (m^2_\pi) } \ .
\label{fpi}
\end{equation}
LQCD results for the $m^2_\pi$ dependence of the pion decay constant $F_\pi$ are
readily available in the literature, e.g. in Ref.~\cite{Durr:2013}. From that
reference, one deduces the values
$F_\pi \approx 112$ MeV at $m_\pi = 410$ MeV,
$F_\pi \approx 129$ MeV at $m_\pi = 570$ MeV, and
$F_\pi \approx 141$ MeV at $m_\pi = 700$ MeV. To obtain the
latter value, a linear $m^2_\pi$ dependence of $F_\pi$ has been assumed,
as suggested by {Fig.~5} in Ref.~\cite{Durr:2013}.
LQCD results for the $m^2_\pi$ dependence of the axial-vector strengths $g^{BB}_A$
can be found in Ref.~\cite{Alexandrou:2016} though only up to
$m_\pi = 500$ MeV. For $g^{{\Sigma}_c\Sc}_A$ a rather moderate increase with $m_\pi$ is
suggested by LQCD~\cite{Alexandrou:2016}. $g^{NN}_A$ is found to be practically
independent of $m_\pi$, though, unfortunately, the lattice results do not match
well with the known value at the physical point. There is no information on
the variation of the ${\Lambda}_c{\Sigma}_c\pi$ vertex with $m_\pi$. Because of these reasons
we neglect the dependence of the $g_A$'s on $m_\pi$ in our calculation and
assume that $g_{A}^{BB'}\equiv g_{A}^{BB'}(m_\pi=138~{\rm MeV})$.
Specifically, we use $g_A^{NN}=1.27$~\cite{PDG},
$g_A^{{\Sigma}_c\Sc}=0.71$~\cite{Alexandrou:2016}
and $g_A^{{\Lambda}_c{\Sigma}_c}=0.74$~\cite{Albertus:2005,Can:2016}.
Note that the by far strongest dependence of $f_{BB'\pi}$ on $m^2_\pi$
comes via $F_\pi (m^2_\pi)$ \cite{Durr:2013} and this circumstance is
adequately taken into account in our calculation.

Since $f_{{\Lambda}_c\Lc\pi} \equiv 0$ under the assumption that isospin is
conserved, there is no one-pion exchange contribution to
the ${\Lambda}_c N \to {\Lambda}_c N$ potential. However, we include the coupling
of ${\Lambda}_c N$ to ${\Sigma}_c N$ via pion exchange, which is known to play an important
role in case of the $\Lambda N$ and $\Sigma N$ systems \cite{Haidenbauer:2013}.
The resulting effective two-pion exchange contribution to the
${\Lambda}_c N$ potential is generated by solving a coupled-channel
Lippmann-Schwinger (LS) equation, see below. In this context let us note that the pertinent
contribution arises anyway at NLO, even in a single-channel treatment \cite{Beane:2003}.
In principle, at NLO there are further contributions from two-pion exchange
\cite{Haidenbauer:2013}. However, in the present study we omit those for simplicity
reasons and assume that they are effectively absorbed into
the contact terms. Furthermore, contributions from $\eta$ and $D-$meson
exchanges that would arise under the assumption of SU(3) (or SU(4)) symmetry,
are likewise delegated to the contact interactions.

As already mentioned above, in line with the analysis of the LQCD data in
Ref.~\cite{Miyamoto:2017}, we treat the ${\Lambda}_c N$ interaction as an effective
single-channel problem. Possible effects from the interaction in the
${\Sigma}_c N \to {\Sigma}_c N$ channel are thereby absorbed into the LECs of the
${\Lambda}_c N$ potential.
A further issue that we have ignored is the role of heavy quark spin
symmetry~\cite{Liu:2011xc,Lu:2017}. Indeed the ${\Sigma}_c N$ and ${{\Sigma}_c^*} N$ thresholds
are just about $65$ MeV apart \cite{PDG} so that the coupling between those systems
could be important \cite{Liu:2011xc}. However, since already the ${{\Sigma}_c} N$ channel is
not explicitly considered, there is no point in including ${{\Sigma}_c^*} N$ in our analysis.
If there are any effects of it, these are likewise absorbed into the ${\Lambda}_c N$ LECs.
Anyway, since $M_{{\Sigma}_c^*}-M_{{\Lambda}_c} = {234}$ MeV,
there should be less influence on the ${\Lambda}_c N$ amplitude. Indeed, because of the larger
mass difference $M_{{\Sigma}_c}-M_{{\Lambda}_c} \approx 167$~MeV, as compared to $M_{\Sigma}-M_{\Lambda}
\approx 78$~MeV, one would expect that even the channel coupling to ${{\Sigma}_c} N$ plays a
less important role for the charm sector than in the strangeness sector.

The reaction amplitudes are obtained from the solution of a coupled-channel
LS equation for the interaction potentials, which is
given in partial-wave projected form by
\begin{eqnarray}
T^{\nu'\nu,J}_{\rho'\rho}(p',p;\sqrt{s})&=& V^{\nu'\nu,J}_{\rho'\rho}(p',p) \nonumber\\
&& + \sum_{{\rho''},\nu''}\int_0^\infty \frac{dp''p''^2}{(2\pi)^3} \,
V^{\nu'\nu'',J}_{\rho'\rho''}(p',p'') \nonumber \\
&&\times\, \frac{2\mu_{\rho''}}{p^2_{\rho}-p''^2+i\eta}
T^{\nu''\nu,J}_{\rho''\rho}(p'',p;\sqrt{s}).
\label{LS}
\end{eqnarray}
Here, the label $\rho$ indicates the particle channels and the label $\nu$ the partial wave.
$\mu_\rho$ is the pertinent reduced mass. The on-shell momentum in the intermediate state,
$p_{\rho}$, is defined by
$\sqrt{s}=\sqrt{m^2_{B_{1,\rho}}+p_{\rho}^2}+\sqrt{m^2_{B_{2,\rho}}+p_{\rho}^2}$.
Following the practice of the HAL QCD Collaboration~\cite{Miyamoto:2017},
phase shifts will be given as functions of the kinetic energy in the ${\Lambda}_c N$
c.m. frame, $E = \sqrt{s} - M_{{\Lambda}_c}- M_N$.

For baryon-baryon potentials constructed along the lines of chiral EFT
(cf. Eqs. (\ref{LEC}) and (\ref{OPE})) a regularization is required when solving the
LS equation (\ref{LS}) \cite{Epelbaum:2005,Machleidt:2011zz}.
Ideally, the resulting reaction amplitude $T$ should be completely independent of
the employed regularization scheme. In practise, there is still a dispute about how
regularization should be done in the application of chiral EFT to $NN$ scattering
(and accordingly in the $YN$ or $Y_cN$ case) to satisfy all formal aspects from a
field theory point of view and there is no generally accepted scheme,
cf. Refs.~\cite{Epelbaum:2005,Machleidt:2011zz} and references therein.
A commonly accepted procedure is the introduction of a cutoff into the
Lippmann-Schwinger equation or (equivalently) to the potential.
The controversial issue is, however, how one should then proceed in
detail in order to achieve the desired cutoff independence of the results
see, e.g.\ \cite{Lepage:1997,Epe05,Nogga:2005,Mach:2013,Phillips:2013,Epe15,Epelbaum:2017}.

In the present work we refrain from elaborating on this still open question.
Rather we follow strictly the procedure used in the application of chiral EFT to $NN$ scattering
by Epelbaum and collaborators \cite{Epe05} and in our $YN$ studies \cite{Haidenbauer:2013},
where the potentials in the LS equation are cut off in momentum space by multiplication
with a regulator function, $f(p',p) = \exp\left[-\left(p'^4+p^4\right)/\Lambda^4\right]$,
so that the high-momentum components of the baryon and pseudoscalar meson fields are
removed. The cutoff parameter $\Lambda$ in the regulator is typically in the order of
$\Lambda \approx 500$ MeV \cite{Epelbaum:2008ga,Machleidt:2011zz}. It is kept finite in the
calculation. Approximate cutoff independence is achieved by going to higher orders
in the perturbative expansion of the potential where the sucessively
arising contact terms allow one to absorb/remove the cutoff dependence
more and more efficiently \cite{Epe15}.
The actual values we employ for the cut-off are $\Lambda =500$-$600$ MeV,
in line with the range that yielded optimal and stable results in our NLO study
of the ${\Lambda} N$ and ${\Sigma} N$ interactions \cite{Haidenbauer:2013}.
The variations of the results with the cutoff, i.e. the residual cutoff dependence,
reflect uncertainties that will be indicated by bands in the plots we show in the
next section.

In the analysis of the LQCD simulations we follow closely the strategy
of our previous works in Refs.~\cite{Haidenbauer:2011a,Haidenbauer:2011b,Haidenbauer:2017}:
(i) The LECs, i.e. the only free parameters in the potential, are determined by
a fit to LQCD results (phase shifts) employing the inherent baryon and meson masses
of the lattice simulation;
(ii) Results at the physical point are obtained via a calculation in which
the pertinent physical masses of the mesons are substituted in the evaluation of
the potential and those of the baryons in the baryon-baryon propagators appearing
in the LS equation.
The baryon masses corresponding to the LQCD simulations at $m_\pi= 410$ and
$570$~MeV are taken from Ref.~\cite{Miyamoto:2017}.
For the calculation at the physical point we use the masses from the
PDG \cite{PDG}, i.e. $M_{{\Lambda}_c} = 2286.5$ MeV, $M_{{\Sigma}_c} = 2455$ MeV.

Alternative ways to implement the quark mass dependence of the nuclear forces can
be found in Refs.~\cite{Baru:2015,Berengut:2013}.

\begin{table*}
\renewcommand{\arraystretch}{1.4}
\centering
\caption{Low-energy constants employed in the present work.
$\tilde C_{^1\!S_0} \equiv {\tilde{C}}_{^1\!S_0} + \tilde{D}_{^1\!S_0} m^2_\pi$
at the specified pion mass, etc., see Eq.~(\ref{LEC}). The values
for $\tilde C$ are in 10$^4$~GeV$^{-2}$, those for $C$ in 10$^4$~GeV$^{-4}$.
}
\label{tab:lec}
\vskip 0.2cm
\begin{center}
\begin{tabular}{||l||c|c||c|c||c|c||}
\hline
\hline
& \multicolumn{2}{c||}{$570$ MeV} & \multicolumn{2}{c||}{$410$ MeV} & \multicolumn{2}{c||}{$138$ MeV} \\
\hline
\hline
{$\Lambda = 500$ MeV}
& $\tilde C$ & $C$ & $\tilde C$ & $C$ & $\tilde C$ & $C$ \\
\hline
$^1S_0$ &  -0.008576 & -0.006105  & -0.01758 & -0.01362  & -0.02615 & -0.02077 \\
$^3S_1$ &  0.03917   & 0.04569    & 0.1409   & 0.1270    & 0.2377   & 0.2043  \\
$^3D-{}^3S_1$ &      &-0.1136     &          & -0.1136   &          & -0.1136 \\
\hline
\hline
{$\Lambda = 600$ MeV}
& $\tilde C$ & $C$ & $\tilde C$ & $C$ & $\tilde C$ & $C$ \\
\hline
$^1S_0$ &   -0.008026 &  0.001564  &  -0.01393 & -0.002252  &  -0.01954 & -0.005879 \\
$^3S_1$ &    0.04304  &0.07553  & 0.1619  &0.1688   & 0.2749  & 0.2574 \\
$^3D-{}^3S_1$ &      &-0.1344     &          & -0.1344   &          & -0.1344 \\
\hline
\hline
\end{tabular}
\end{center}
\end{table*}
\renewcommand{\arraystretch}{1.0}

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

LQCD results for phase shifts are available for the ${\Lambda}_c N$ $^1\!S_0$ and $^3\!S_1$
partial waves for $m_\pi = 410,\, 570,\, 700$ MeV \cite{Miyamoto:2017}.
We determine the LECs of the contact interaction, cf. Eq.~(\ref{LEC}),
by a fit to the lattice data at the two lower pion masses. Specifically,
the pion-mass dependence exhibited by the LQCD simulation is exploited to determine
the LECs $\tilde D_i$ and $D_i$ that encode the pion-mass dependence of the contact
interaction.
The fits are done to the phase shifts, generated from the parameterized version
of the ${\Lambda}_c N$ potentials provided in Ref.~\cite{Miyamoto:2017}, for energies
up to $30$~MeV.
Alternative fits taking into account the HAL QCD results up
to $50$~MeV were performed too. In this context it should be noted that applications
of chiral EFT to $NN$ scattering and specifically to the $^1\!S_0$ $np$ state reveal
that NLO interactions are expected to provide quantitative results up to roughly
$50$ MeV \cite{Epelbaum:2008ga,Machleidt:2011zz}.
The actual values of the LECs at the fitted pion masses and
at the physical value are summarized in Table~\ref{tab:lec}.

\begin{figure*}[t]
\begin{center}
\includegraphics[height=80mm,angle=-90]{A0E.ps}
\includegraphics[height=80mm,angle=-90]{Ca0.ps}
\caption{Results for the ${\Lambda}_c N$ $^1\!S_0$ partial wave.
Left: Phase shifts of the HAL QCD Collaboration \cite{Miyamoto:2017}
at $m_\pi = 570$ MeV (blue open circles) and $410$ MeV (green filled circles),
as a function of the c.m. kinetic energy, together with our fits (lines/bands).
The broader (red) band is the prediction for $m_\pi = 138$ MeV.
Right: Dependence of the $^1\!S_0$ scattering length on the pion mass.
The bands represent the cutoff variation $\Lambda = 500-600$ MeV, see text.
}
\label{ph1s0E}
\end{center}
\end{figure*}

Results for the $^1\!S_0$ partial wave are presented in Fig.~\ref{ph1s0E}.
The phase shifts for pion masses $570$, $410$ and $138$ MeV are shown
on the left side while the dependence of the scattering length on the
pion mass is depicted on the right side.
The bands represent the dependence of the results on variations of the
cutoff $\Lambda$.
One can see that the lattice results at $m_\pi=410$~MeV are reproduced
quantitatively by our potential up to c.m. kinetic energies of around $40$~MeV,
as expected for an NLO interaction, while those at $570$~MeV are remarkably
well described over the whole energy range shown.
In both cases the cutoff dependence is negligibly small at low energies so that
the bands are hardly visible.
The phase shift obtained from our interaction when extrapolated to the physical
point (red bands) do exhibit a noticeable but still rather moderate cutoff
dependence. A maximum of around 20 degrees of the $^1\!S_0$ phase shift is
predicted in this case.
The pion-mass dependence of the $^1\!S_0$ scattering length,
shown on the right-hand side of Fig.~\ref{ph1s0E}, is fairly smooth and almost
linear in $m_\pi$. Only close to the physical point a somewhat
stronger $m_\pi$ dependence is visible. The value predicted at
$138$ MeV is $a= -0.85 {\,\cdot\cdot\cdot\,} -1.0$~fm.

\begin{figure*}[t]
\begin{center}
\includegraphics[height=80mm,angle=-90]{A1E.ps}
\includegraphics[height=80mm,angle=-90]{Ca1.ps}
\caption{Results for the ${\Lambda}_c N$ $^3\!S_1$ partial wave.
Left: Phase shifts of the HAL QCD Collaboration \cite{Miyamoto:2017}
at $m_\pi = 570$ MeV (blue open circles) and $410$ MeV (green filled circles),
as a function of the c.m. kinetic energy, together with our fits (lines/bands).
The broader (red) band is the prediction for $m_\pi = 138$ MeV.
Right: Dependence of the $^3\!S_1$ scattering length on the pion mass.
The bands represent the cutoff variation $\Lambda = 500-600$ MeV, see text.
}
\label{ph3s1E}
\end{center}
\end{figure*}

Results for the $^3\!S_1$ partial wave are presented in Fig.~\ref{ph3s1E}.
Again phase shifts as well as the pion-mass dependence of the scattering
length are shown. Those are very similar to the ones of the $^1\!S_0$ state,
even on a quantitative level. Indeed, it was already noted by the HAL
QCD Collaboration that the corresponding ($^1S_0$ and $^3S_1$) potentials
they extracted are almost identical -- at 410 as well as at 570 MeV \cite{Miyamoto:2017},
and we see that this characteristic feature persists even in our extrapolation to the
physical point. This is certainly a remarkable feature in view of the fact that in case
of the $^3\!S_1$ partial wave there is a coupling to the $^3\!D_1$ induced by the
tensor force, and this coupling is taken into account in our analysis.

\begin{table*}
\renewcommand{\arraystretch}{1.2}
\centering
\caption{${\Lambda}_c N$ scattering lengths in the $^1\!S_0$ ($a_s$) and $^3\!S_1$ ($a_t$)
partial waves (in fm) for different pion masses. Predictions for the corresponding
effective range parameters $r$ (in fm) at the physical point are likewise given.
}
\label{tab:a}
\vskip 0.2cm
\begin{center}
\begin{tabular}{ l|l||c|c|c|c }
\hline
\hline
\multicolumn{2}{c||}{$m_\pi$} & $700$ MeV & $570$ MeV & $410$ MeV & $138$ MeV \\
\hline
\hline
$a_s$ & HAL QCD \cite{Miyamoto:2017} & $-0.13\pm 0.11$ & $-0.24\pm 0.13$ & $-0.49\pm 0.18$ & $-$ \\
        &our results    & $-0.04{\,\cdot\cdot\cdot\,} -0.13$ & $-0.21$ & $ -0.45 $ & $-0.85{\,\cdot\cdot\cdot\,} -1.00$ \\
$r_s$ &                              &                 & &        & $2.88{\,\cdot\cdot\cdot\,} 2.61$ \\
\hline
$a_t$ & HAL QCD \cite{Miyamoto:2017} & $\quad -0.17\pm 0.10\quad$ & $\quad -0.29\pm 0.16\quad$
& $\quad -0.51\pm 0.20\quad$ & $-$ \\
 & our results & $-0.09{\,\cdot\cdot\cdot\,} -0.15$ & $-0.25{\,\cdot\cdot\cdot\,} -0.26$ & $ -0.46{\,\cdot\cdot\cdot\,} -0.47 $ & $-0.81{\,\cdot\cdot\cdot\,} -0.98$ \\
$r_t$ &                              &                 & &                   & $3.50{\,\cdot\cdot\cdot\,} 3.15$ \\
\hline
\hline
\end{tabular}
\end{center}
\end{table*}
\renewcommand{\arraystretch}{1.0}

For completeness we have summarized the results for the $^1\!S_0$ and $^3\!S_1$
scattering lengths, $a_s$ and $a_t$, in Table~\ref{tab:a}.
From those numbers one can see that the variation in the extrapolated values due
to the employed regularization scheme is in the order of $0.2$ fm.
Additional fits to the phase shifts where the considered energy range was extended
up to $50$ MeV led to changes in the scattering lengths at the physical point of
around $0.1$ fm, with a clear tendency to smaller values.
Further exploratory fits carried out by us indicate that variations in the
scattering length of $\pm 0.05$~fm at $m_\pi = 410$~MeV amount to differences
of about $\pm 0.10$~fm at $m_\pi = 138$ MeV. Combining these observations with
the uncertainty of $\pm 0.2$~fm given by the HAL QCD Collaboration for their
result at $m_\pi = 410$ MeV suggests that the scattering lengths
at the physical point could be about $0.3$~fm larger, i.e. as large as $-1.3$~fm.
Finally, for estimating possible effects from the ${\Sigma}_c N$ channel we performed
calculations with the pion-exchange contribution (\ref{OPE}) to the ${\Sigma}_c N$ potential
included. Adding its contribution within a coupled-channel framework (\ref{LS})
has a negligible influence on the results at $m_\pi = 410$ MeV and $570$ MeV, but
leads to noticeable variations in the predictions for $m_\pi = 138$ MeV that
amount to roughly $0.3$~fm in the scattering lengths.
Of course, a complete calculation has to include also a contact term in the ${\Sigma}_c N$
channel, which could facilitate a reduction of the effect.
Obviously, additional information on the ${\Sigma}_c N$ phase shifts, as had been
provided in the preliminary study of the HAL QCD Collaboration \cite{Miyamoto:2016A},
and promised in Ref.~\cite{Miyamoto:2017} for the future,
would be helpful for reducing the uncertainty in the extrapolation.
In any case, a more quantitative estimate of the overall uncertainty should
be attempted, once lattice data with better statistics are
available and/or results for pion masses closer to the physical point.

Clearly, and in line with the LQCD results for larger pion masses, we do not
get any ${\Lambda}_c N$ bound states.
However, what are the consequences of the results presented in the preceding
paragraphs for the possible existence of bound ${\Lambda}_c$ nuclei?
Let us look
at the strangeness sector and, specifically, at the lightest $\Lambda$ nucleus
that is experimentally observed, namely the hypertriton $^3\!_\Lambda \rm H$.
The experimental value for the $^3\!_\Lambda \rm H$ binding energy is
$(-2.354\pm 0.050)$ MeV,
which implies a separation energy for the $\Lambda$ of only $(0.13\pm 0.05)$ MeV.
Faddeev calculations of the coupled $\Lambda NN$--$\Sigma NN$ three-body systems
for realistic $YN$ potentials have been reported in Refs.~\cite{Nogga:2002,Haidenbauer:2007,Nogga:2013}.
Those suggest that, for $YN$ interactions which provide sufficient attraction
so that the hypertriton is bound, the $\Lambda N$ scattering lengths are in the order of
$-2.9$ to $-2.6$ fm for the $^1\!S_0$ state and $-1.5$ to $-1.7$ fm for $^3\!S_1$.

The scattering lengths for the ${\Lambda}_c N$ system at the physical point, deduced
from the LQCD simulations of the HAL QCD Collaboration, are $-0.8$ to $-1.0$ fm
for the $^1\!S_0$ as well as the $^3\!S_1$ states.
Thus, in the case of the $^1\!S_0$ partial wave the value is considerably smaller
than its counterpart in the strangeness sector while for the $^3\!S_1$ channel the
difference is less dramatic. Since the average $\Lambda N$ (or ${\Lambda}_c N$) potential
that is relevant for the hypertriton is dominated by the spin-singlet channel, i.e.
$V_{\Lambda N}\,(V_{{\Lambda}_c N}) \propto \frac{3}{4}V_s+\frac{1}{4}V_t$, see Ref.~\cite{Gibson:1994}
for details, it seems rather unlikely that an only moderately attractive singlet ${\Lambda}_c N$
interaction can support the existence of a $^{\ \, 3}_{{\Lambda}_c} \rm H$ bound state.
Even effects due to a reduction of the kinetic energy associated with the ${\Lambda}_c$
induced by the larger mass of the ${\Lambda}_c$ as compared to the $\Lambda$, emphasized
in several works in the past \cite{Garcilazo:2015qha,Dover:1977jw,Gibson:1983zw},
cannot compensate for this large difference in the strength.
Actually, in Ref.~\cite{Gibson:1983zw} an estimate for the binding energy of charmed
three- and four-body systems is provided based on an exact solution of
corresponding scattering equations. Interestingly, Model~4 considered in that work
yields scattering lengths very similar to the ones of our analysis, namely
$a_s = -1.075$ fm and $a_t = -0.828$ fm (cf. Table I in that work). For this interaction
no bound state is found for the $^{\ \, 3}_{{\Lambda}_c} {\rm H}$ nucleus. However, the 4-body
systems $^{\ \, 4}_{{\Lambda}_c} {\rm He}$ and $^{\ \, 4}_{{\Lambda}_c} {\rm Li}$ could be already
bound, though possibly only weakly, even when considering the additional
repulsive effect from the Coulomb force due to the charge of the
${\Lambda}_c^+$ as discussed in Ref.~\cite{Gibson:1983zw}.

With regard to recent few-body calculations, the ${Y}_c N$ interaction employed in
Ref.~\cite{Maeda:2015hxa} is strongly attractive and leads already to ${\Lambda}_c N$
bound states. The binding energies for ${\Lambda}_c NN$
are then in the order of $20$ MeV. Those results do not allow us to draw any
conclusions about what would happen for significantly less attractive ${Y}_c N$
interactions like the one we deduce from the analysis of the LQCD simulations
by the HAL QCD Collaboration.
In case of the three-body calculation presented in Ref.~\cite{Garcilazo:2015qha},
a charmed hypertriton with $J=3/2$ (and total isospin $I=0$) is predicted.
For that state the ${\Lambda}_c N$ spin-triplet interaction is dominant~\cite{Miyagawa:1995}.
However, since there is no information on the scattering length and effective
range parameters of the ${\Lambda}_c N$ potential employed in
Ref.~\cite{Garcilazo:2015qha}, a proper assessment of that result is difficult.

\section{Summary}

In the present work, we have used the framework of chiral effective field theory
to extrapolate lattice QCD results for the $\Lambda_c N$ interaction at $m_\pi=410
- 570$~MeV by the HAL QCD Collaboration~\cite{Miyamoto:2017} to the physical point.
Thereby, we have followed a strategy employed previously in the
strangeness sector~\cite{Haidenbauer:2011a,Haidenbauer:2011b,Haidenbauer:2017}.
However, contrary to the calculations performed in those works, no recourse to
SU(3) or SU(4) flavor symmetry has been made in the present study.
Furthermore, in the present work the pion-mass dependence of all components that
constitute the $\Lambda_c N$ potential up to next-to-leading order
(pion-exchange diagrams and four-baryon contact terms)
are taken into account. Information from lattice QCD simulations is utilized to
implement these features.

Our analysis of the HAL QCD results points to a moderately attractive $\Lambda_c N$
interaction at the physical point with scattering lengths of $a \approx -1$~fm for
the $^1\!S_0$ as well as for the $^3\!S_1$ partial waves.
Such an interaction leads to the possibility of bound four- and/or five-baryons
systems with a $\Lambda_c$ baryon and presumably of heavier ${\Lambda}_c$ hypernuclei.
On the other hand, two-body (${\Lambda}_c N$) bound states as advocated in some recent
investigations based on phenomenological potentials~\cite{Liu:2011xc,Maeda:2015hxa}
can be definitely excluded, even if one considers uncertainties in the extrapolation
of the lattice results.
Also, the existence of a hypertriton-like $\Lambda_c NN$ three-body bound state
(with $J=1/2$) seems rather unlikely if one takes past investigations
\cite{Gibson:1983zw} as benchmark.

\vskip 0.2cm \noindent
{\bf Acknowledgments} \\
Work partially supported by Conselho Nacional de Desenvolvimento Cient\'{i}fico
e Tecnol\'{o}gico (CNPq), Grant No. 305894/2009-9 (G.K.) and Grant No. 464898/2014-5 (G.K.)
(INCT F\'{\i}sica  Nuclear e Applica\c{c}\~oes), Funda{\c c}\~{a}o de Amparo \`{a} Pesquisa
do Estado de S\~{a}o Paulo (FAPESP), Grant No. 2013/01907-0 (G.K.)


\section{Introduction}
Automatic image analysis is becoming an increasingly important technique to support clinical diagnosis and treatment planning.
It is usually formulated as a classification problem where medical imaging abnormalities
are identified as different clinical conditions~\cite{shi2017multimodal,chen2015glaucoma, shin2016learning, wang2017detecting,zhang2017mdnet}.
In clinical practice, visual evidence that supports the classification result, such as spatial localization~\cite{akselrod2016region} or segmentation~\cite{zhao2016multiscale, zilly2017glaucoma} of sites of abnormalities is an indispensable part of clinical diagnosis which provides interpretation and insights. Therefore, it is of vital importance that the image analysis method is able to provide both classification results and the associated visual evidence with high accuracy.

 \begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/overview.pdf}
    \caption{Overview of our chest X-ray image analysis network for thoracic disease diagnosis. The network reads chest X-ray images and produces prediction scores and localization for the diseases}
    \label{fig:overview}
    \vspace{-1.5em}
 \end{figure}

Figure~\ref{fig:overview} is an overview of our approach. We focus on chest X-ray image analysis. Our goal is to both classify the clinical conditions and identify the abnormality locations. A chest X-ray image might contain multiple sites of abnormalities with monotonous and homogeneous image features. This often leads to the inaccurate classification of clinical conditions. It is also difficult to identify the sites of abnormalities because of their variances in the size and location.
For example, as shown in Figure \ref{fig:chestimage},
the presentation of ``Atelectasis'' (alveoli are deflated down) is usually limited to local regions of a lung \cite{gylys2017medical} but possible to appear anywhere on both sides of lungs; while ``Cardiomegaly'' (enlarged heart) always covers half of the chest and is always around the heart.

The lack of large-scale datasets also stalls the advancement of automatic chest X-ray diagnosis.
Wang \etal provides one of the largest publicly available chest x-ray datasets with disease labels\footnote{While abnormalities, findings, clinical conditions, and diseases have distinct meanings in the medical domain, here we simply refer them as diseases and disease labels for the focused discussion in computer vision.}
along with a small subset with region-level annotations (bounding boxes) for evaluation~\cite{wang2017chestx}\footnote{The method proposed in \cite{wang2017chestx} did not use the bounding box information for localization training.}.
As we know, the localization annotation is much more informative than just a single disease label to improve the model performance as demonstrated in~\cite{liu2017attention}. However, getting detailed disease localization annotation can be difficult and expensive. Thus, designing models that can work well with only a small amount of localization annotation is a crucial step for the success of
clinical applications.

In this paper, we present a unified approach that simultaneously improves disease identification and localization with only a small amount of X-ray images containing disease location information. Figure~\ref{fig:overview} demonstrates an example of the output of our model. Unlike the standard object detection task in computer vision, we do not strictly predict bounding boxes.
Instead, we produce regions that indicate the diseases, which aligns with the purpose of visualizing and interpreting the disease better. Firstly, we apply a CNN to the input image so that the model learns the information of the entire image and implicitly encodes both the class and location information for the disease~\cite{redmon2016you}. We then slice the image into a patch grid to capture the local information of the disease.
For an image with bounding box annotation, the learning task becomes a fully supervised problem since the disease label for each patch can be determined by the overlap between the patch and the bounding box. For an image with only a disease label, the task is formulated as a multiple instance learning (MIL) problem~\cite{babenko2008multiple}---at least one patch in the image belongs to that disease. If there is no disease in the image, all patches have to be disease-free. In this way, we have unified the disease identification and localization into the same underlying prediction model but with two different loss functions.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{figures/chestimage.pdf}
  \caption{Examples of chest X-ray images with the disease bounding box. The disease regions are annotated in the yellow bounding boxes by radiologists. }
  \label{fig:chestimage}
  \vspace{-1.5em}
\end{figure}

We evaluate the model on the aforementioned chest X-ray image dataset provided in~\cite{wang2017chestx}.
Our quantitative results show that the proposed model achieves significant accuracy improvement over the published state-of-the-art on both disease identification and localization, despite the limited number of bounding box annotations of a very small subset of the data. In addition, our qualitative results reveal a strong correspondence between the radiologist's annotations and detected disease regions, which might produce further interpretation and insights of the diseases.

\section{Related Work}
\textbf{Object detection.}
Following the R-CNN work \cite{girshick2014rich}, recent progresses has focused on processing all regions with only one shared CNN \cite{he2014spatial,girshick2015fast}, and on eliminating explicit
region proposal methods by directly predicting the bounding boxes.
In \cite{ren2015faster}, Ren \etal developed a region proposal network (RPN) that regresses from anchors to regions of interest (ROIs).
However, these approaches could not be easily used for images without enough annotated bounding boxes.
To make the network process images much faster, Redmon \etal proposed a grid-based object detection network, YOLO, where an image is partitioned into $S\times S$ grid cells, each of which is responsible to predict the coordinates and confidence scores of $B$ bounding boxes \cite{redmon2016you}. The classification and bounding box prediction are formulated into one loss function to learn jointly. A step forward, Liu \etal partitioned the image into multiple grids with different sizes proposing a multi-box detector overcoming the weakness in YOLO and achieved better performance \cite{liu2016ssd}.
Similarly, these approaches are not applicable for the images without bounding boxes annotation.
Even so, we still adopt the idea of handling an image as a group of grid cells and treat each patch as a classification target.

\textbf{Medical disease diagnosis.}
Zhang \etal proposed a dual-attention model using images and optional texts to make accurate prediction \cite{zhang2017tandemnet}.
In \cite{zhang2017mdnet}, Zhang \etal proposed an image-to-text model to establish a direct mapping from medical images to diagnostic reports.
Both models were evaluated on a dataset of bladder cancer images and corresponding diagnostic reports.
Wang \etal took advantage of a large-scale chest X-ray dataset to formulate the disease diagnosis problem as multi-label classification, using class-specific image feature transformation \cite{wang2017chestx}. They also applied a thresholding method to the feature map visualization \cite{zeiler2014visualizing} for each class and derived the bounding box for each disease.
Their qualitative results showed that the model usually generated much larger bounding box than the ground-truth.
Hwang \etal \cite{hwang2016self} proposed a self-transfer learning framework to learn localization from the globally pooled class-specific feature maps supervised by image labels.
These works have the same essence with class activation mapping \cite{zhou2016learning} which handles natural images.
The location annotation information was not directly formulated into the loss function in the none of these works.
Feature map pooling based localization did not effectively capture the precise disease regions.

\textbf{Multiple instance learning.}
In multiple instance learning (MIL), an input is a labeled bag (e.g., an image) with many instances (e.g., image patches)~\cite{babenko2008multiple}. The label is assigned at the bag level.
Wu \etal assumed each image as a dual-instance example, including its object proposals and possible text annotations \cite{wu2015deep}.  The framework achieved convincing performance in vision tasks including classification and image annotation. In medical imaging domain, Yan \etal utilized a deep MIL framework for body part recognition \cite{yan2016multi}.
Hou \etal first trained a CNN on image patches and then an image-level decision fusion model by patch-level prediction histograms to generate the image-level labels \cite{hou2016patch}.
By ranking the patches and defining three types of losses for different schemes, Zhu \etal proposed an end-to-end deep multi-instance network to achieve mass classification for whole mammogram \cite{zhu2017deep}.
We are building an end-to-end unified model to make great use of both image level labels and bounding box annotations effectively.

\begin{figure*}
    \centering
    \vspace{-1em}
    \includegraphics[width=1.9\columnwidth]{figures/architecture.pdf}
    \caption{Model overview. (a) The input image is firstly processed by a CNN. (b) The patch slicing layer resizes the convolutional features from the CNN using max-pooling or bilinear interpolation. (c) These regions are then passed to a fully-convolutional recognition network. (d) During training, we use multi-instance learning assumption to formulate two types of images; during testing, the model predicts both labels and class-specific localizations. The red frame represents the ground truth bounding box. The green cells represent patches with positive labels, and brown is negative. Please note during training, for unannotated images, we assume there is at least one positive patch and the green cells shown in the figure are not deterministic.}
    \label{fig:architecture}
    \vspace{-1.5em}
 \end{figure*}

\section{Model}

Given images with disease labels and limited bounding box information, we aim to design a unified model that simultaneously produces disease identification and localization. We have formulated two tasks into the same underlying prediction model so that 1) it can be jointly trained end-to-end and 2) two tasks can be mutually beneficial. The proposed architecture is summarized in Figure \ref{fig:architecture}.

\subsection{Image model}\label{subsec:imagemodel}
\textbf{Convolutional neural network.}
As shown in Figure~\ref{fig:architecture}(a), we use the residual neural network (ResNet) architecture \cite{he2016deep} given its dominant performance in ILSVRC competitions \cite{russakovsky2015imagenet}. Our framework can be easily extended to any other advanced CNN models. The recent version of pre-act-ResNet \cite{he2016identity} is used (we call it ResNet-v2 interchangeably in this paper).
After removing the final classification layer and global pooling layer, an input image with shape $h\times w\times c$ produces a feature tensor with shape $h'\times w'\times c'$ where $h$, $w$, and $c$ are the height, width, and number of channels of the input image respectively while $h'=\frac{h}{32}$, $w'=\frac{w}{32}$, $c'=2048$.
The output of this network encodes the images into a set of abstracted feature maps.

\textbf{Patch slicing.}
Our model divides the input image into $P \times P$ patch grid,
and for each patch, we predict $K$ binary class probabilities, where $K$ is the number of possible disease types.
As the CNN gives $c'$ input feature maps with size of $h' \times w'$, we down/up sample the input feature maps to $P \times P$ through a patch slicing layer shown in Figure~\ref{fig:architecture}(b).
Please note that $P$ is an adjustable hyperparameter.
In this way, a node in the same spatial location across all the feature maps corresponds to one patch of the input image.
We upsample the feature maps If their sizes are smaller than expected patch grid size.
Otherwise, we downsample them.

\textit{Upsampling.}
We use a simple bilinear interpolation to upsample the feature maps to the desired patch grid size.
As interpolation is, in essence, a fractionally stridden convolution, it can be performed in-network for end-to-end learning and is fast and effective \cite{long2015fully}.
A deconvolution layer \cite{zeiler2014visualizing} is not necessary to cope with this simple task.

\textit{Downsampling.}
The bilinear interpolation makes sense for downsampling only if the scaling factor is close to $1$.
We use max-pooling to down sample the feature maps.
In general cases, the spatial size of the output volume is a function of the input width/height ($w$), the filter (receptive field) size ($f$), the stride ($s$), and the amount of zero padding used ($p$) on the border.
The output width/height ($o$) can be obtained by $\frac{w-f+2p}{s}+1$.
To simplify the architecture, we set $p=0$ and $s=1$, so that $f=w-o+1$.

\textbf{Fully convolutional recognition network.}\label{subsubsec:recognition_net}
We follow \cite{long2015fully} to use fully convolution layers as the recognition network.
Its structure is shown in Figure~\ref{fig:architecture}(c).
The $c'$ resized feature maps are firstly convolved by $3\times3$ filters into a smaller set of feature maps with $c^\ast$ channels, followed by batch normalization \cite{ioffe2015batch} and rectified linear units (ReLU) \cite{glorot2011deep}. Note that the batch normalization also regularizes the model.
We set $c^\ast=512$ to represent patch features.
The abstracted feature maps are then passed through a $1\times 1$ convolution layer to generate a set of $P \times P$ final predictions with $K$ channels.
Each channel gives prediction scores for one class among all the patches,
and the prediction for each class is normalized by a logistic function (sigmoid function) to $[0,1]$.
The final output of our network is the $P \times P \times K$ tensor of predictions.
The image-level label prediction for each class in $K$ is calculated across $P \times P$ scores, which is described in Section~ \ref{subsec:formulation}.

\subsection{Loss function}\label{subsec:formulation}
\textbf{Multi-label classification.}
Multiple disease types can be often identified in one chest X-ray image and these disease types are not mutually exclusive. Therefore, we define a binary classifier for each class/disease type in our model. The binary classifier outputs the class probability. Note that the binary classifier is not applied to the entire image, but to all small patches. We will show how this can translate to image-level labeling below.

\textbf{Joint formulation of localization and classification.}
Since we intend to build $K$ binary classifiers, we will exemplify just one of them, for example, class $k$. Note that $K$ binary classifiers will use the same features and only differ in their last logistic regression layers. The $i$th image $x_i$ is partitioned into a set $\mathcal{M}$ of patches equally, $x_i = [x_{i1}, x_{i2}, ..., x_{im}]$, where $m=|\mathcal{M}|=P\times P$

\textit{Images with annotated bounding boxes.}
As shown in Figure~\ref{fig:architecture}(d), suppose an image is annotated with class $k$ and a bounding box. We denote $n$ be the number of patches covered by the bounding box, where $n<m$. Let this set be $\mathcal{N}$.
Each patch in the set $\mathcal{N}$ as positive for class $k$ and each patch outside the bounding box as negative.
Note that if a patch is covered partially by the bounding box of class $k$, we still consider it a positive patch for class $k$. The bounding box information is not lost.
For the $j$th patch in $i$th image, let $p_{ij}^k$ be the foreground probability for class $k$. Since all patches have their labels,
the probability of an image being positive for class $k$ is defined as,
\begin{align}\label{eqn:bbox_prod}
\small
 \textstyle  p(y_k|x_i, {\rm bbox}^k_i) = \prod_{j\in \mathcal{N}} p_{ij}^k\cdot \prod_{j \in \mathcal{M} \setminus \mathcal{N}} (1-p_{ij}^k),
\end{align}
where $y_k$ is the $k$th network output denoting whether an image is a positive example of class $k$.
For example, for a class other than $k$, this image is treated as the negative sample without a bounding box. We define a patch as positive to class $k$ when it is overlapped with a ground-truth box, and negative otherwise.

\textit{Images without annotated bounding boxes.} If the $i$th image is labeled as class $k$ without any bounding box,
we know that there must be at least one patch classified as $k$ to make this image a positive example of class $k$.
Therefore, the probability of this image being positive for class $k$ is defined as the image-level score
\footnote{Later on, we notice an similar definition \cite{liao2017evaluate} for this multi-instance problem. We argue that our formulation is in a different context of solving classification and localization in a unified way for images with limited bounding box annotation.
Yet, this related work can be viewed as a successful validation of our multi-instance learning based formulation.},
\begin{align}\label{eqn:nobbox_prod}
\small
 \textstyle   p(y_k|x_i) = 1-\prod_{j\in \mathcal{M}}(1-p_{ij}^k).
\end{align}
At test time, we calculate $p(y_k|x_i)$ by Eq.~\ref{eqn:nobbox_prod} as the prediction probability for class $k$.

\textit{Combined loss function.} Note that $p(y_k|x_i, {\rm bbox}^k_i)$ and $p(y_k|x_i)$ are the image-level probabilities. The loss function for class $k$ can be expressed as minimizing the negative log likelihood of all observations as follows,
\begin{small}
\begin{align}
    \mathcal{L}_k =  & \textstyle - \lambda_{{\rm bbox}}\sum_i \eta_i p(y^\ast_k|x_i, {\rm bbox}^k_i)\log(p(y_k|x_i,{\rm bbox}^k_i)) \nonumber \\
    &- \textstyle  \lambda_{{\rm bbox}}\sum_i \eta_i (1 - p(y^\ast_k|x_i, {\rm bbox}^k_i))\log(1-p(y_k|x_i,{\rm bbox}^k_i)) \nonumber\\
    & \textstyle - \sum_i (1-\eta_i) p(y^\ast_k|x_i)\log(p(y_k|x_i)) \nonumber \\
    & \textstyle - \sum_i (1-\eta_i)(1 - p(y^\ast_k|x_i))\log(1-p(y_k|x_i)), \label{eqn:loss}
\end{align}
\end{small}where $i$ is the index of a data sample, $\eta_i$ is $1$ when the $i$th sample is annotated with bounding boxes, otherwise $0$. $\lambda_{{\rm bbox}}$ is the factor balancing the contributions from annotated and unannotated samples.
$p(y^\ast_k|x_i) \in \{0,1\}$ and $p(y^\ast_k|x_i, {\rm bbox}^k_i)\in \{0,1\}$ are the observed probabilities for class $k$. Obviously, $p(y^\ast_k|x_i, {\rm bbox}^k_i)\equiv1$, thus equation \ref{eqn:loss} can be re-written as follows,
\begin{small}
\begin{align}
\mathcal{L}_k = &-\textstyle \lambda_{{\rm bbox}}\sum_i \eta_i \log(p(y_k|x_i,{\rm bbox}^k_i)) \nonumber\\
    & - \textstyle \sum_i (1-\eta_i) p(y^\ast_k|x_i)\log(p(y_k|x_i)) \nonumber \\
    & - \textstyle\sum_i (1-\eta_i)(1 - p(y^\ast_k|x_i))\log(1-p(y_k|x_i)). \label{eqn:loss2}
\end{align}
\end{small}
In this way, the training is strongly supervised (per patch) by the given bounding box; it is also supervised by the image-level labels if the bounding boxes are not available.

To enable end-to-end training across all classes, we sum up the class-wise loss to define the total loss as,
\begin{align*}
\small
 \textstyle   \mathcal{L} = \sum_k\mathcal{L}_k.
\end{align*}

\subsection{Localization generation}
The full model predicts a probability score for each patch in the input image.
We define a score threshold $T_s$ to distinguish the activated patches against the non-activated ones.
If the probability score $p_{ij}^k$ is larger than $T_s$, we consider the $j$th patch in the $i$th image belongs to the localization for class $k$.
We set $T_s=0.5$ in this work.
Please note that we do not predict strict bounding boxes for the regions of disease---the combined patches representing the localization information can be a non-rectangular shape.

\subsection{Training}\label{subsec:training}
We use ResNet-v2-50 as the image model and select the patch slicing size from $\{12,16,20\}$.
The model is pre-trained on the ImageNet 1000-class dataset \cite{deng2009imagenet} with Inception \cite{szegedy2015going} preprocessing method where the image is normalized to $[-1,1]$ and resized to $299\times299$.
We initialize the CNN with the weights from the pre-trained model, which helps the model converge faster than training from scratch.
During training, we also fine-tune the image model, as we believe the feature distribution of medical images differs from that of natural images.
We set the batch size as $5$ to load the entire model to the GPU, train the model with $500$k iterations of minibatch, and decay the learning rate by $0.1$ from $0.001$ every $10$ epochs of training data.
We add L2 regularization to the loss function to prevent overfitting.
We optimize the model by Adam \cite{kingma2014adam} method with asynchronous training on 5 Nvidia P100 GPUs. The model is implemented in TensorFlow \cite{tensorflow2015-whitepaper}.

\begin{table*}
    \centering
    \vspace{-1.1em}
	\resizebox{1.7\columnwidth}{!}{
		\begin{tabular}{cccccccc}
		\hline\hline
		\multicolumn{1}{c|}{Disease}       & Atelectasis               & Cardiomegaly              & Consolidation             & Edema                     & Effusion                   & Emphysema                 & Fibrosis                   \\ \hline
\multicolumn{1}{c|}{baseline} & 0.70                      & 0.81                      & 0.70                      & 0.81                      & 0.76                       & 0.83                      & \textbf{0.79}                       \\ \hline
\multicolumn{1}{c|}{ours}          & $\textbf{0.80} \pm 0.00 $ & $\textbf{0.87} \pm 0.01 $ & $ \textbf{0.80}\pm 0.01 $ & $ \textbf{0.88}\pm 0.01 $ & $ \textbf{0.87}\pm 0.00 $  & $ \textbf{0.91}\pm 0.01 $ & $ 0.78\pm 0.02 $  \\ \hline
\multicolumn{1}{c|}{Disease}       & Hernia                    & Infiltration              & Mass                      & Nodule                    & Pleural Thickening         & Pneumonia                 & Pneumothorax               \\ \hline
\multicolumn{1}{c|}{baseline} & \textbf{0.87}                            & 0.66                      & 0.69                      & 0.67                      & 0.68                       & 0.66                      & 0.80                       \\ \hline
\multicolumn{1}{c|}{ours}          & $ 0.77 \pm 0.03 $         & $ \textbf{0.70}\pm 0.01 $ & $ \textbf{0.83}\pm 0.01 $ & $ \textbf{0.75}\pm 0.01 $ & $ \textbf{0.79} \pm 0.01 $ & $ \textbf{0.67}\pm 0.01 $ & $ \textbf{0.87} \pm 0.01 $ \\ \hline
		\end{tabular}
	}
\caption{AUC scores comparison with the reference baseline model. Results are rounded to two decimal digits for table readability. Bold values denote better results. The results for the reference baseline are obtained from the latest update of \cite{wang2017chestx}.}
\label{tbl:classification_comparison}
\vspace{-1.5em}
\end{table*}

\textbf{Smoothing the image-level scores.}
In Eq.~\ref{eqn:bbox_prod} and \ref{eqn:nobbox_prod}, the notation $\prod$ denotes the product of a sequence of probability terms ($[0,1]$), which often leads to the a product value of $0$ due to the computational underflow if $m=|\mathcal{M}|$ is large. The log loss in Eq.~\ref{eqn:loss} mitigates this for Eq.~\ref{eqn:bbox_prod}, but does not help Eq.~\ref{eqn:nobbox_prod}, since the log function can not directly affect its product term.
To mitigate this effect, we normalize the patch scores $p_{ij}^k$ and $1 - p_{ij}^k$ from $[0,1]$ to $[0.98,1]$ to make sure the image-level scores $p(y_k|x_i, {\rm bbox}^k_i)$ and $p(y_k|x_i)$ smoothly varies within the range of $[0,1]$.
Since we are thresholding the image-level scores in the experiments, we found this normalization works quite well. See supplementary material for a more detailed discussion on this.

\textbf{More weights on images with bounding boxes.}
In Eq.~\ref{eqn:loss2}, the parameter $\lambda_{\rm bbox}$ weighs the contribution from the images with annotated bounding boxes. Since the amount of such images is limited, and if we treat them equally with the images without bounding boxes, it often leads to worse performance. We thus increase the weight for images with bounding boxes
to $\lambda_{{\rm bbox}}=5$ by cross validation.

\section{Experiments}
\textbf{Dataset and preprocessing.}
NIH Chest X-ray dataset \cite{wang2017chestx} consists of $112,120$ frontal-view X-ray images with 14 disease labels (each image can have multi-labels). These labels are obtained by analyzing the associated radiology reports.
The disease labels are expected to have accuracy of above $90\%$~\cite{wang2017chestx}. We take the provided labels as ground-truth for training and evaluation in this work.
Meanwhile, the dataset also contains $984$ labelled bounding boxes for $880$ images by board-certified radiologists. Note that the provided bounding boxes correspond to only 8 types of disease instances. We separate the images with provided bounding boxes from the entire dataset. Hence we have two sets of images called ``annotated'' ($880$ images) and ``unannotated'' ($111,240$ images).

We resize the original 3-channel images from $1024\times1024$ to $512\times512$ pixels for fast processing. The pixel values in each channel are normalized to $[-1,1]$. We do not apply any data augmentation techniques

\begin{figure*}
  \centering
  \vspace{-1em}
  \includegraphics[width=1.95\columnwidth]{figures/auc_scores_revised.pdf}
  \caption{AUC scores for models trained using different data combinations. Training set: annotated samples,\{left: $80\%$ ($704$ images), right: $0\%$ (baseline, $0$ images)\} for each disease type; unannotated samples, \{80\% ($88,892$), 60\% ($66,744$), 40\% ($44,496$), 20\% ($22,248$), 0\%($0$)\} from left to right for each disease type. The evaluation set is $20\%$ annotated and unannotated samples which are not included in the training set. No result for 0\% annotated and 0\% unannotated images. Using $80\%$ annotated images and certain amount of unannotated images improves the AUC score compared to using the same amount of unannotated images (same colored bars in two groups for the same disease), as the joint model benefits from the strong supervision of the tiny set of bounding box annotations.}
  \label{fig:classification}
  \vspace{-1.5em}
\end{figure*}
\subsection{Disease identification}\label{subsec:classification}
We conduct a 5-fold cross-validation.
For each fold, we have done two experiments.
In the first one, we train the model using $70\%$ of bounding-box annotated and $70\%$ unannotated images to compare the results with the reference model~\cite{wang2017chestx} (Table~\ref{tbl:classification_comparison}). To our knowledge, the reference model has the published state-of-the-art performance of disease identification on this dataset.
In the second experiment, we explore two data ratio factors of annotated and unannotated images to demonstrate the effectiveness of the supervision provided by the bounding boxes (Figure~\ref{fig:classification}).
We decrease the amount of images without bounding boxes from $80\%$ to $0\%$ by a step of $20\%$. And then for each of those settings, we train our model by adding $80\%$ or none of bounding-box annotated images.
For both experiments, the model is always evaluated on the fixed $20\%$ annotated and unannotated images for this fold.

\textbf{Evaluation metrics.} We use AUC scores (the area under the ROC\footnote{Here ROC is the Receiver Operating Characteristic, which measures the true positive rate (TPR) against the false positive rate (FPR) at various threshold settings ($200$ thresholds in this work).} curve) to measure the performance of our model~\cite{fawcett2006introduction}. A higher AUC score implies a better classifier.

\textbf{Comparison with the reference model.}
Table \ref{tbl:classification_comparison} gives the AUC scores for all the classes.
We compare our results with the reference baseline fairly:
we, as the reference, use ImageNet pre-trained ResNet-50
\footnote{Using ResNet-v2~\cite{he2016identity} shows marginal performance difference for our network compared to ResNet-v1~\cite{he2016deep} used in the reference baseline.}
, after which a convolution layer follows; both works use $70\%$ images for training and $20\%$ for testing, and we also conduct a 5-fold cross-validation to show the robustness of our model.

Compared to the reference model, our proposed model achieves better AUC scores for most diseases. The overall improvement is remarkable and the standard errors are small.
The large objects, such as ``Cardiomegaly'', ``Emphysema'', and ``Pneumothorax'', are as well recognized as the reference model.
Nevertheless, for small objects like ``Mass'' and ``Nodule'', the performance is significantly improved. Because our model slices the image into small patches and uses bounding boxes to supervise the training process, the patch containing small object stands out of all the patches to represent the complete image.
For ``Hernia'', there are only $227$ (0.2\%) samples in the dataset. These samples are not annotated with bounding boxes. Thus, the standard error is relatively larger than other diseases.

\textbf{Bounding box supervision improves classification performances.}
We consider using $0\%$ annotated images as our own baseline (right groups in Figure~\ref{fig:classification}).
We use $80\%$ annotated images (left groups in Figure~\ref{fig:classification}) to compare with the our own baseline.
We plot the mean performance for the cross-validation in Figure~\ref{fig:classification}, the standard errors are not plotted but similar to the numbers reported in Table~\ref{tbl:classification_comparison}.
The number of $80\%$ annotated images is just $704$, which is quite small compared to the number of $20\%$ unannotated images ($22,248$).
We observe in Figure \ref{fig:classification} that for almost all the disease types, using $80\%$ annotated images to train the model improves the prediction performance (by comparing the bars with the same color in two groups for the same disease).
For some disease types, the absolute improvement is significant ($>5\%$).
We believe that this is because all the disease classifiers share the same underlying image model; a better-trained image model using eight disease annotations can improve all 14 classifiers' performance. Specifically, some diseases, annotated and unannotated, share similar visual features. For example, ``Consolidation'' and ``Effusion'' both appear as fluid accumulation in the lungs, but only ``Effusion'' is annotated. The feature sharing enables supervision for ``Effusion'' to improve ``Consolidation'' performance as well.

\textbf{Bounding box supervision reduces the demand of the training images.}
Importantly, it requires less unannotated images to achieve the similar AUC scores by using a small set of annotated images for training.
As denoted with red circles in Figure~\ref{fig:classification}, taking ``Edema'' as an example, using $40\%$ ($44,496$) unannotated images with $80\%$ ($704$) annotated images ($45,200$ in total) outperforms the performance of using only $80\%$ ($88,892$) unannotated images.

\textbf{Discussion.}
Generally, decreasing the amount of unannotated images (from left to right in each bar group) will degrade AUC scores accordingly in both groups of $0\%$ and $80\%$ annotated images.
Yet as we decrease the amount of unannotated images, using annotated images for training gives smaller AUC degradation or even improvement.
For example, we compare the ``Cardiomegaly'' AUC degradation for two pairs of experiments: \{annotated:80\%, unannotated:80\% and 20\%\} and \{annotated:0\%, unannotated:80\% and 20\%\}. The AUC degradation for the first group is just $0.07$ while that for the second group is $0.12$ (accuracy degradation from blue to yellow bar).

When the amount of unannotated images is reduced to $0\%$, the performance is significantly degraded.
Because under this circumstance, the training set only contains positive samples for eight disease types and lacks the positive samples of the other six.
Interestingly, ``Cardiomegaly'' achieves the second best score (AUC = $0.8685$, the second green bar in Figure~\ref{fig:classification}) when only annotated images are trained. The possible reason is that the location of cardiomegaly is always fixed to the heart covering a large area of the image and the feature distributions for enlarged hearts are similar to normal ones.
Without unannotated samples, the model easily distinguishes the enlarged hearts from normal ones given supervision from bounding boxes.
When the model sees hearts without annotation, the enlarged ones are disguised and fail to be recognized. As more unannotated samples are trained, the enlarged hearts are recognized again by image-level supervision (AUC from $0.8086$ to $0.8741$).

\begin{figure*}[t]
  \centering
  \vspace{-1em}
  \includegraphics[width=1.9\columnwidth]{figures/loc_ior_1.pdf}
  \caption{Disease localization accuracy using IoR where T(IoR)=0.1. Training set:  annotated samples, \{$0\%$ ($0$), $20\%$ ($176$), $40\%$ ($352$), $60\%$ ($528$), $80\%$ ($704$)\} from left to right for each disease type; unannotated samples, $100\%$ ($111,240$ images). The evaluation set is $20\%$ annotated samples which are not included in the training set. For each disease, the accuracy is increased from left to right, as we increase the amount of annotated samples, because more annotated samples bring more bounding box supervision to the joint model.}
  \label{fig:localization_ior_1}
  \vspace{-0.5em}
\end{figure*}
\begin{figure*}
  \centering
  \includegraphics[width=1.9\columnwidth]{figures/loc_ior_2.pdf}
  \caption{Disease localization accuracy using IoR where T(IoR)=0.1. Training set:  annotated samples, $80\%$ ($704$ images); unannotated samples,  \{$0\%$ ($0$), $20\%$ ($22,248$), $40\%$ ($44,496$), $60\%$ ($66,744$), $80\%$ ($88,892$), $100\%$ ($111,240$)\} from left to right for each disease type. The evaluation set is $20\%$ annotated samples which are not included in the training set. Using annotated samples only can produce a model which localizes some diseases. As the amount of unannotated samples increases in the training set, the localization accuracy is improved and all diseases can be localized. The joint formulation for both types of samples enables unannotated samples to improve the performance with weak supervision.  }
  \label{fig:localization_ior_2}
  \vspace{-1.5em}
\end{figure*}

\begin{table*}[t]
    \centering
    \vspace{-1em}
	\resizebox{1.85\columnwidth}{!}{
		\begin{tabular}{c|ccccccccc}
		\hline
T(IoU)               & \multicolumn{1}{c|}{Model} & Atelectasis               & Cardiomegaly              & Effusion                 & Infiltration             & Mass                     & Nodule                  & Pneumonia              & Pneumothorax            \\ \hline
\multirow{2}{*}{0.1} & \multicolumn{1}{c|}{ref.}  & 0.69                      & 0.94                      & 0.66                     & 0.71                     & 0.40                     & 0.14                    & \textbf{0.63}                   & 0.38                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.71}\pm 0.05 $  & $\textbf{0.98} \pm 0.02 $ & $\textbf{0.87} \pm 0.03$ & $\textbf{0.92} \pm 0.05$ & $\textbf{0.71} \pm0.10$  & $\textbf{0.40} \pm0.10$ & $0.60\pm 0.11$        & $\textbf{0.63}\pm0.09$ \\ \hline
\multirow{2}{*}{0.2} & \multicolumn{1}{c|}{ref.}  & 0.47                      & 0.68                      & 0.45                     & 0.48                     & 0.26                     & 0.05                    & 0.35                   & 0.23                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.53}\pm 0.05 $  & $\textbf{0.97} \pm 0.02 $ & $\textbf{0.76} \pm 0.04$ & $\textbf{0.83}\pm 0.06 $ & $\textbf{0.59} \pm 0.10$ & $\textbf{0.29}\pm 0.10$ & $\textbf{0.50}\pm0.12$ & $\textbf{0.51}\pm0.08$  \\ \hline
\multirow{2}{*}{0.3} & \multicolumn{1}{c|}{ref.}  & 0.24                      & 0.46                      & 0.30                     & 0.28                     & 0.15                     & 0.04                    & 0.17                   & 0.13                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.36}\pm 0.08$   & $\textbf{0.94} \pm 0.01$  & $\textbf{0.56}\pm 0.04$  & $\textbf{0.66}\pm 0.07$  & $\textbf{0.45}\pm 0.08$  & $\textbf{0.17}\pm0.10$  & $\textbf{0.39}\pm0.09$ & $\textbf{0.44}\pm0.10$  \\ \hline
\multirow{2}{*}{0.4} & \multicolumn{1}{c|}{ref.}  & 0.09                      & 0.28                      & 0.20                     & 0.12                     & 0.07                     & 0.01                    & 0.08                   & 0.07                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.25}\pm 0.07$   & $\textbf{0.88}\pm 0.06$   & $\textbf{0.37}\pm 0.06$  & $\textbf{0.50}\pm 0.05$  & $\textbf{0.33}\pm 0.08$  & $\textbf{0.11}\pm0.02$  & $\textbf{0.26}\pm0.07$ & $\textbf{0.29}\pm0.06$  \\ \hline
\multirow{2}{*}{0.5} & \multicolumn{1}{c|}{ref.}  & 0.05                      & 0.18                      & 0.11                     & 0.07                     & 0.01                     & 0.01                    & 0.03                   & 0.03                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.14}\pm 0.05$   & $\textbf{0.84}\pm 0.06$   & $\textbf{0.22}\pm 0.06$  & $\textbf{0.30}\pm 0.03$  & $\textbf{0.22}\pm 0.05$  & $\textbf{0.07}\pm0.01$  & $\textbf{0.17}\pm0.03$ & $\textbf{0.19}\pm0.05$  \\ \hline
\multirow{2}{*}{0.6} & \multicolumn{1}{c|}{ref.}  & 0.02                      & 0.08                      & 0.05                     & 0.02                     & 0.00                     & 0.01                    & 0.02                   & 0.03                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.07}\pm 0.03$   & $\textbf{0.73}\pm 0.06$   & $\textbf{0.15}\pm 0.06$  & $\textbf{0.18}\pm 0.03$  & $\textbf{0.16}\pm 0.06$  & $\textbf{0.03}\pm0.03$  & $\textbf{0.10}\pm0.03$ & $\textbf{0.12}\pm0.02$  \\ \hline
\multirow{2}{*}{0.7} & \multicolumn{1}{c|}{ref.}  & 0.01                      & 0.03                      & 0.02                     & 0.00                     & 0.00                     & 0.00                    & 0.01                   & 0.02                    \\ \cline{2-10}
                     & \multicolumn{1}{c|}{ours}  & $\textbf{0.04}\pm 0.01$   & $\textbf{0.52}\pm 0.05$   & $\textbf{0.07}\pm 0.03$  & $\textbf{0.09}\pm0.02$   & $\textbf{0.11}\pm0.06 $  & $\textbf{0.01}\pm 0.00$ & $\textbf{0.05}\pm0.03$ & $\textbf{0.05}\pm0.03$  \\ \hline
		\end{tabular}
	}
\caption{Disease localization accuracy comparison using IoU where T(IoU)=\{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7\}.The bold values denote the best results. Note that we round the results to two decimal digits for table readability. Using different thresholds, our model outperforms the reference baseline in most cases and remains capability of localizing diseases when the threshold is big.}
\label{tbl:localization_iou_comparison}
\vspace{-0.5em}
\end{table*}

 \begin{figure*}
    \centering
    \includegraphics[width=1.95\columnwidth]{figures/qualitative.pdf}
    \caption{Example localization visualization on the test images. The visualization is generated by rendering the final output tensor as heatmaps and overlapping on the original images. We list some thoracic diseases as examples. The left image in each pair is the original chest X-ray image and the right one is the localization result. All examples are positive for corresponding labels. We also plot the ground-truth bounding boxes in yellow on the results when they are provided in the dataset.}
    \label{fig:qualitative}
    \vspace{-1.5em}
 \end{figure*}

\subsection{Disease localization}
Similarly, we conduct a 5-fold cross-validation.
For each fold, we have done three experiments.
In the first experiment, we investigate the importance of bounding box supervision by using all the unannotated images and increasing the amount of annotated images from $0\%$ to $80\%$ by the step of $20\%$ (Figure~\ref{fig:localization_ior_1}).
In the second one, we fix the amount of annotated images to $80\%$ and increase the amount of unannotated images from $0\%$ to $100\%$ by the step of $20\%$ to observe whether unannotated images are able to help annotated images to improve the performance (Figure~\ref{fig:localization_ior_2}).
At last, we train the model with $80\%$ annotated images and half ($50\%$) unannotated images to compare localization accuracy with the reference baseline~\cite{wang2017chestx} (Table~\ref{tbl:localization_iou_comparison}).
For each experiment, the model is always evaluated on the fixed $20\%$ annotated images for this fold.

\textbf{Evaluation metrics.}
We evaluate the detected regions (which can be non-rectangular and discrete) against the annotated ground truth (GT) bounding boxes,
using two types of measurement: intersection over union ratio (IoU) and intersection over the detected region (IoR)
\footnote{Note that we treat discrete detected regions as one prediction region, thus IoR is analogous to intersection over the detected bounding box area ratio (IoBB).}.
The localization results are only calculated for those eight disease types with ground truth provided.
We define a correct localization when either IoU $>$ T(IoU) or IoR $>$ T(IoR), where T(*) is the threshold.

\textbf{Bounding box supervision is necessary for localization.}
We present the experiments shown in Figure~\ref{fig:localization_ior_1}. The threshold is set as tolerable as T(IoR)=0.1 to show the training data combination effect on the accuracy. Please refer to the supplementary material for localization performance with T(IoU)=0.1, which is similar to Figure~\ref{fig:localization_ior_1}.
Even though the amount of the complete set of unannotated images is dominant compared with the evaluation set ($111,240$ v.s. $176$),
without annotated images (the most left bar in each group), the model fails to generate accurate localization for most disease types. Because in this situation, the model is only supervised by image-level labels and optimized using probabilistic approximation from patch-level predictions.
As we increase the amount of annotated images gradually from $0\%$ to $80\%$ by the step of $20\%$ (from left to right in each group), the localization accuracy for each type is increased accordingly.
We can see the necessity of bounding box supervision by observing the localization accuracy increase.
Therefore, the bounding box is necessary to provide accurate localization results and the accuracy is positively proportional to the amount of annotated images.
We have similar observations when T(*) varies.

\textbf{More unannotated data does not always mean better results for localization.}
In Figure~\ref{fig:localization_ior_2}, when we fix the amount of annotated images and increase the amount of unannotated ones for training (from left to right in each group), the localization accuracy does not increase accordingly. Some disease types achieve very high accuracy (even highest) without any unannotated images (the most left bar in each group), such as ``Pneumonia'' and ``Cardiomegaly''.
Similarly as described in the discussion of Section \ref{subsec:classification}, unannotated images and too many negative samples degrade the localization performance for these diseases.
All disease types experience an accuracy increase, a peak score, and then an accuracy fall (from orange to green bar in each group).
Therefore, with bounding box supervision, unannotated images will help to achieve better results in some cases and it is not necessary to use all of them.

\textbf{Comparison with the reference model.}
In each fold, we use $80\%$ annotated images and $50\%$ unannotated images to train the model and evaluate on the other $20\%$ annotated images in each fold. Since we use 5-fold cross-validation, the complete set of annotated images has been evaluated to make a relatively fair comparison with the reference model.
In Table \ref{tbl:localization_iou_comparison}, we compare our localization accuracy under varying T(IoU) with respect to the reference model in \cite{wang2017chestx}.
Please refer to the supplementary material for the comparison between our localization performance and the reference model with varying T(IoR).
Our model predicts accurate disease regions, not only for the easy tasks like ``Cardiomegaly'' but also for the hard ones like ``Mass'' and ``Nodule'' which have very small regions.
When the threshold increases, our model maintains a large accuracy lead over the reference model. For example, when evaluated by T(IoU)=0.6, our ``Cardiomegaly'' accuracy is still $73.42\%$ while the reference model achieves only $16.03\%$; our ``Mass'' accuracy is $14.92\%$ while the reference model fails to detect any ``Mass'' ($0\%$ accuracy).
In clinical practice, a specialist expects as accurate localization as possible so that a higher threshold is preferred.
Hence, our model outperforms the reference model with a significant improvement with less training data.
Please note that as we consider discrete regions as one predicted region, the detected area and its union with GT bboxs are usually larger than the reference work which generates multiple bounding boxes. Thus for some disease types like ``Pneumonia'', when the threshold is small, our result is not as good as the reference.

\subsection{Qualitative results}
Figure \ref{fig:qualitative} shows exemplary localization results of the unified diagnosis model.
The localization enables the explainability of chest X-ray images.
It is intuitive to see that our model produces accurate localization for the diseases compared with the given ground-truth bounding boxes.
Please note for ``Infiltration'' (3th and 4th images in the 3rd row of Figure \ref{fig:qualitative}), both sides of lungs for this patient are infiltrated. Since the dataset only has one bounding box for one disease per image, it misses annotating other bounding boxes for the same disease. Our model gives the remedy.
Even though the extra region decreases the IoR/IoU score in the evaluation, but in clinical practice, it provides the specialist with suspicious candidate regions for further examination.
When the localization results have no ground-truth bounding boxes to compare with, there is also a strong consistency between our results and radiological signs.
For example, our model localizes the enlarged heart region (1st and 2nd images in the 2nd row) which implies ``Cardiomegaly'', and the lung peripheries is highlighted (5th and 6th images in the 2nd row) implying ``Fibrosis'' which is in accordance with the radiological sign of the net-like shadowing of lung peripheries.
The ``Edema'' (7th and 8th images in the 1st row) and ``Consolidation'' (7th and 8th images in the 2nd row) are accurately marked by our model. ``Edema'' always appears in an area that is full of small liquid effusions as the example shows. ``Consolidation'' is usually a region of compressible lung tissue that has filled with the liquid which appears as a big white area. The model successfully distinguishes both diseases which are caused by similar reason.

\section{Conclusion}
We propose a unified model that jointly models disease identification and localization with limited localization annotation data. This is achieved through the same underlying prediction model for both tasks.
Quantitative and qualitative results demonstrate that our method significantly outperforms the state-of-the-art algorithm.

{\small
\bibliographystyle{ieee}