# Arbitrarily high order implicit ODE integration by correcting a neural
network approximation with Newton’s method

Daniel W. Crews  Computational Plasma Dynamics Lab, University of Washington,
Seattle, WA

###### Abstract

As a method of universal approximation deep neural networks (DNNs) are capable
of finding approximate solutions to problems posed with little more
constraints than a suitably-posed mathematical system and an objective
function. Consequently, DNNs have considerably more flexibility in
applications than classical numerical methods. On the other hand they offer an
uncontrolled approximation to the sought-after mathematical solution. This
suggests that hybridization of classical numerical methods with DNN-based
approximations may be a desirable approach. In this work a DNN-based
approximator inspired by the physics-informed neural networks (PINNs)
methodology is used to provide an initial guess to a Newton’s method iteration
of a very-high order implicit Runge-Kutta (IRK) integration of a nonlinear
system of ODEs, namely the Lorenz system. In the usual approach many explicit
timesteps are needed to provide a guess to the implicit system’s nonlinear
solver, requiring enough work to make the IRK method infeasible. The DNN-based
approach described in this work enables large implicit time-steps to be taken
to any desired degree of accuracy for as much effort as it takes to converge
the DNN solution to within a few percent accuracy. This work also develops a
general formula for the matrix elements of the IRK method for an arbitrary
quadrature order.

##  1 Introduction

Neural networks are universal approximators, so it seems inevitable that they
take their place in the toolbox of the working applied mathematician and
engineer. Setting up a neural net approximation for a given problem is fairly
easy compared to classical methods of discretization. The result looks like
the correct answer to the eye, which for many purposes may be good enough. Yet
at their current stage of development neural networks seem to provide an
uncontrolled approximation in the sense that there is no obvious relationship
between error, iterations, and the nature of the input. In contrast, what one
might call classical methods, like Newton’s root-finding algorithm, are a more
controlled approximation, but require a good guess to get started. This
project aims to work with the strength of both techniques by hybridizing a
neural network with classical numerical analysis.

Recently, neural networks have been used to solve PDEs in a framework called
physics-informed neural networks (PINNs) [1]. The authors constrained the PINN
loss function to satisfy a PDE via automatic differentiation as a continuous
approximation, and experimented with discrete approximations by informing the
PINN with classical numerical methods like implicit Runge-Kutta (IRK) time-
steppers. However, the output only looks like the correct solution, as the
neural network loss does not converge beyond an L2 error of $\approx 1\%$ [2].
To solve the problem of convergence, this project uses Newton’s method, a
beloved root-finding algorithm, to correct the neural network prediction for
an arbitrarily high order IRK method. While Newton’s method requires
calculation of a system’s Jacobian matrix, if one knows physics with which to
inform a neural network, then one can also calculate the Jacobian of said
physics!

Usually a less-accurate explicit time-step calculates the guess for an
implicit solution. However, explicit time-steps can be unstable for large
steps (e.g. the famous CFL condition), thereby limiting the implicit step-
size. Yet the IRK method on Gauss-Legendre points is A-stable, in Butcher’s
terminology, to any order [3]. The implicit time-stepper can take very large
steps provided that it has a good guess. The advantage of a PINN approach in
providing the guess for Newton’s method is to overcome the step-size
limitation of explicitly-predicted stages.

Section 2 describes the hybrid IRK-integrator method applied to the Lorenz
system, which is chosen for the chaoticity of its trajectories. High-order
methods are desirable to ensure that one remains on the correct trajectory
when taking a large step. Appendix A develops a closed-form expression for the
matrix elements of the IRK method to arbitrary order, and Appendix B then
develops the IRK Jacobian for Newton’s method. If desired, the implementation
of all work presented in this report may be found on the author’s Github page
using Python and Tensorflow 2.

##  2 Prediction of IRK stages with a neural network and correction

In the following, $m$ denotes the step number, Greek characters
$\alpha,\beta,\cdots$ denote vector components, and Roman characters
$i,j,\cdots$ the RK stages. Following Raissi et al., a physics-informed neural
net (PINN) can learn the IRK stages given the input and an appropriate
constraint [1]. In the context of the system of ODEs
$dy^{\alpha}/dt=f(\bm{y},t)$ the neural network takes the state $y^{m,\alpha}$
as input and for each input has $(n+1)$ outputs, the $n$ RK stages
$y^{\alpha}_{i}$ plus the full step output $y^{m+1,\alpha}$. To formulate the
loss function, the IRK system for the step $y^{m,\alpha}\to y^{m+1,\alpha}$ is
rearranged into the following form,

$\displaystyle y^{m,\alpha}$
$\displaystyle=y^{\alpha}_{i}-h\sum_{j=0}^{n-1}A^{j}_{i}f_{j}^{\alpha}(\bm{y}),$
(1) $\displaystyle y^{m,\alpha}$
$\displaystyle=y^{m+1,\alpha}-h\sum_{i=0}^{n-1}\frac{w_{i}}{2}f_{i}^{\alpha}(\bm{y}).$
(2)

The neural net is then trained on its own input $y^{m,\alpha}$ in an
autoencoder-like fashion so that the outputs
$\\{y_{i}^{\alpha},y^{m+1,\alpha}\\}$ satisfy the system of Eqs. 1 and 2 under
a mean-square error minimization objective. In this way the only data required
is that of the initial condition, and not the solution at any later time.
Requiring only the initial data is a desirable property in an integrator.
After training to a suitable loss, the network output is given to Newton’s
method as its initial guess.

(a) Activation function tanh.

(b) Activation function ELU.

Figure 1: An orbit in the Lorenz system integrated using a fifty-stage IRK
step of step-size $h=0.75$, solved by Newton iteration based on the orbit
predicted by a neural network approximation. Shown is a comparison of neural
net predictions using hyperbolic tangent (tanh) and exponential linear unit
(ELU) given 10,000 training epochs, and the result of Newton iteration using
the neural net output as the initial guess. Both guesses allow Newton’s method
to take over and produce an overall L2 error for the implicit system of
$10^{-10}$. Clearly, ELU works better here. Newton’s method is not observed to
converge given a worse guess, such as explicit forward Euler to each RK time-
stage.

###  2.1 Application to the Lorenz system

For an application, consider the famous Lorenz system of ODEs,

$\displaystyle\frac{dx}{dt}$ $\displaystyle=\sigma(y-x),$ (3)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=x(\rho-z)-y,$ (4)
$\displaystyle\frac{dz}{dt}$ $\displaystyle=xy-\beta z,$ (5)

with “original” parameters $\sigma=10$, $\beta=8/3$, and $\rho=28$. The Lorenz
system’s Jacobian matrix is

$J^{\alpha}_{\beta}=\begin{bmatrix}-\sigma&\sigma&0\\\ \rho&-1&-x\\\
y&x&-\beta\end{bmatrix}.$ (6)

Choosing an initial condition $q_{0}=[10.54,4.112,35.82]$ [4], a neural net is
trained on the initial condition vector of shape $(3,1)$. The net has three
hidden layers of three neurons each, and an output layer of shape $(3,n+1)$.
That such a small number of neurons in each layer ends up working is, to this
novice at least, quite remarkable!

To choose the activation function, Fig. 1 considers fitting the IRK system for
a time-step of $h=0.75$ on a 50-stage IRK method using both tanh and ELU
activations. The function ELU is deemed the winner, at least for this
application, and is used subsequently.

(a) Training epochs: 10,000.

(b) Training epochs: 40,000.

Figure 2: One-hundred stage IRK solution by Newton-iteration of a neural net
prediction of 10,000 and 40,000 epochs. The Newton solver converges in each
case, though a damping factor of $\gamma=0.25$ is necessary. Convergence of
the neural net is slow because the optimization algorithm becomes stuck in
local minima, possibly that of a nearby trajectory. The neural net prediction
in each case is good enough for the Newton solver to converge, which does not
occur if an explicit update is used to guess. The Newton solver is much faster
than waiting for the neural net to converge, which suggests the utility of
methods which hybridize ML techniques with classical numerical methods.

####  2.1.1 Pushing the limit with 100-stage implicit RK advances

Theoretically, a one-hundred stage implicit RK advance on Gauss-Legendre nodes
is $\mathcal{O}(h^{200})$-accurate. Thus, one can take very large steps and
still obtain a quite accurate solution for all times $t\in[t^{m},t^{m+1}]$ as
defined by the interpolating polynomial of the RK stages. However, it is not
known by this author to what extent the floating-point precision of the
quadrature nodes, weights, and RK matrix elements allow machine precision to
be actually obtained in the numerical solution. This will hopefully be the
subject of a future study.

Increasing the stages to one hundred, and the time-step towards $h=1$, is
observed to be difficult for both training the neural net and for iterating
the net output using Newton’s method. However, solutions are obtained with
sufficient training epochs on the neural net and by cranking down the damping
factor $\gamma$ on Newton’s method. A reasonably large time-step is $h=0.80$,
with a truncation error of order $\mathcal{O}((0.8)^{200})\approx 10^{-20}$,
far below machine precision.

Yet, just to push the limits, Fig. 2 investigates the case $h=0.95$. The
neural net did converge to the correct solution, but required a large number
of steps. Even with the good guess provided by the neural net, a significant
damping factor of $\gamma=0.25$ was necessary for Newton’s method to converge
at this large step. A 100-stage step of step-size $h=1$ was also obtained as
an experiment, but required a large amount of computational effort. Further,
any $h>1$ should amplify the error rather than decrease it, so this was not
investigated.

Figure 3 shows the uncorrected neural net outputs along with their combination
into the objective function (the initial condition) as an aid to visualize the
method. For illustration, this is shown for a case which is slowly converging,
as it is stuck in a local minimum. When converged, all the blue dots will
coincide with the initial condition $q_{0}$.

Figure 3: The neural net objective is to linearly combine the output stages
through the RK matrix to match the initial condition. This figure visualizes
the technique for a 100-stage IRK method which is stuck in a local minimum,
with the output stages shown in green and the objective points in blue. The
method is convergent when the blue points all match the initial condition.

Finally, the neural net predictor-Newton’s method corrector 100-stage IRK
time-stepper was applied to multiple time-steps. Figure 4 shows the Lorenz
attractor traced by the trajectory integrated out to ten steps each of step-
size $h=0.8$. While the IRK method with step-size $h=0.8$ has theoretical
truncation error $\mathcal{O}(10^{-20})$, the Newton iteration L2 error
threshold is set to $10^{-10}$.

(a) Perspective on Lorenz attractor.

(b) Another attractor perspective.

Figure 4: A trajectory of the Lorenz system integrated with ten time-steps of
a 100-stage implicit RK Gauss-Legendre method of step-size $h=0.8$, solved by
Newton iteration from a neural net prediction of the RK stages. Each color
shown corresponds to the stages of a single step.

##  3 Conclusions and extensions

Neural nets have immense promise as a universal approximation system,
broadened by their being relatively easy to set-up for a wide variety of
problems with automatic differentiation baked-in. However, the objective
function landscape is full of local minima. While the result of training a
neural net generally looks like the correct answer, there are usually
convergence issues due to lack of constraints in the posed optimization
problem.

This project combined the neural net output with an approximation of classical
numerical analysis, namely Newton’s method. The neural network prediction is
used as a guess in order to solve the nonlinear system associated with an
arbitrarily high-order implicit time-integration of a system of ODEs. Assisted
by the good guesses of the neural network, Newton’s method was found to
converge for an $\mathcal{O}(h^{200})$-accurate method.

While this study looked at the Lorenz system, there’s no reason that this ODE
couldn’t correspond to the semi-discrete ODE system of a spatially-discretized
PDE. The only difficulty of such extensions is to calculate the system
Jacobian, which is shown in Appendix B to consist of the Runge-Kutta
coefficient matrix and the Jacobian of the ODE system. As Appendix A
establishes a closed-form expression for the RK coefficient matrix, to apply
the method to any system of ODEs one need only calculate its Jacobian. Further
extensions of this project can look at more sophisticated iterative methods
which incorporate information on second derivatives, such as BFGS.

##  4 Acknowledgements

The author would like to thank J. Bakarji, S.L. Brunton, J.N. Kutz, U.
Shumlak, and A.D. Stepanov for teaching and helpful discussions.

## References

  * [1] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving partial differential equations.  J. Comp. Phys., 328, 2019.
  * [2] A.D. Jagtap and G.E. Karniadakis.  Extended physics-informed neural networks (XPINNs): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations.  Commun. Comput. Phys., 28, 2020.
  * [3] J.C. Butcher.  Numerical Methods for Ordinary Differential Equations.  John Wiley & Sons, Second edition, 2008.
  * [4] MGolden96.  Gauss–Legendre method; A Practical Example.  Wikipedia, June 2021\.
  * [5] S. Hassani.  Mathematical Physics: A Modern Introduction to its Foundations.  Springer, Second edition, 2013.
  * [6] I. Stakgold.  Boundary Value Problems of Mathematical Physics, volume I.  The Macmillan Company, 1967.
  * [7] C. Dunkl and Y. Xu.  Orthogonal polynomials of several variables.  Cambridge University Press, Second edition, 2014.
  * [8] R. Bulirsch and J. Stoer.  Introduction to numerical analysis.  Springer, New York, 1991.
  * [9] Nodes and Weights of Gauss-Legendre Calculator.  https://kesian.casio.com/, June 2021.
  * [10] A. Townsend.  The race for high order Gauss-Legendre quadrature.  SIAM News, 2015.
  * [11] I. Bogaert.  Iteration-free computation of Gauss-Legendre quadrature nodes and weights.  SIAM J. Sci. Comput., 36, 2014.

##  Appendix A Implicit RK integration using Gauss-Legendre quadrature

This section considers the implicit Runge-Kutta (RK) method using Gauss-
Legendre (GL) quadrature and develops a closed-form expression for the RK
coefficient matrix. First, recall that

$\int_{-1}^{1}f(\xi)d\xi\approx\sum_{i=0}^{n-1}w_{i}f(\xi_{i})$ (7)

with $\\{w_{i}\\}_{i=0}^{n-1}$ quadrature weights and
$\\{\xi_{i}\\}_{i=0}^{n-1}$ the roots of $P_{n}(\xi)$, the $n$’th Legendre
polynomial. The integration is exact for polynomials of order $\leq 2n-1$.
There is a beautiful theory of the interpolation polynomials through GL
quadrature points; it is best understood in terms of the Legendre polynomial
completeness relation [5]

$\sum_{s=0}^{\infty}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(x)P_{s}(y)=\delta(y-x),\quad\quad-1\leq
x,y\leq 1.$ (8)

######  Theorem 1.

The interpolation polynomials $\ell_{j}(\xi)$ through $n$ Gauss-Legendre nodes
satisfy

$\ell_{j}(\xi)=w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{j})P_{s}(\xi),\quad\text{
and }\quad\ell_{j}(\xi_{i})=\delta_{ij}$ (9)

with $\\{w_{j}\\}_{j=0}^{n-1}$, $\\{\xi_{j}\\}_{j=0}^{n-1}$ the quadrature
weights and nodes respectively111Considering the summation
$\lim_{j\to\infty}\sum_{j=0}^{n-1}\ell_{j}(\xi)f(\xi)$ suggests that the the
interpolation polynomial plays a role in discrete integration as if it were a
$\delta$-sequence on $[-1,1]$ in the sense of distributions [6]..

Proof: The interpolation polynomials are defined in Lagrange form as

$\ell_{j}(\xi)=\prod_{\begin{subarray}{c}k=0\\\ k\neq
j\end{subarray}}^{n-1}\frac{\xi-\xi_{k}}{\xi_{j}-\xi_{k}}.$ (10)

Now, note that in addition to continuous orthogonality

$\int_{-1}^{1}P_{k}(\xi)P_{\ell}(\xi)d\xi=\frac{1}{k+\frac{1}{2}}\delta_{k\ell}$
(11)

the Legendre polynomials are discretely orthogonal [7] [8],

$\sum_{s=0}^{n-1}w_{s}P_{k}(\xi_{s})P_{\ell}(\xi_{s})=\frac{1}{k+\frac{1}{2}}\delta_{k\ell}$
(12)

with $w_{s}$ the quadrature weights. This follows from the order $2n-1$
quadrature property as the sum corresponds to the continuous integral
$\int_{-1}^{1}P_{k}(x)P_{\ell}(x)dx$. Expansions in the Lagrange basis of Eqn.
10 and in the Legendre basis are linearly related due to the interpolation
property $\ell_{j}(\xi_{i})=\delta_{ij}$,

$f(\xi)=\sum_{j=0}^{n-1}f_{j}\ell_{j}(\xi)=\sum_{k=0}^{n-1}c_{k}P_{k}(\xi),\quad\implies\quad
f_{j}=\mathcal{V}_{j}^{k}c_{k}$ (13)

where $\mathcal{V}^{k}_{j}=P_{k}(\xi_{j})$ is termed the (generalized)
Vandermonde matrix. The inverse transform follows directly from Eqn. 12 as
$(\mathcal{V}^{-1})_{j}^{k}=w_{j}P_{k}(\xi_{j})(k+\frac{1}{2})$. In
particular, the Lagrange interpolant itself is

$\ell_{j}(\xi)=w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{j})P_{s}(\xi)\approx
w_{j}\delta(\xi-\xi_{j}).$ (14)

###  A.1 The Runge-Kutta method for Gauss-Legendre quadrature

The IRK coefficient matrix for GL stage-points is developed to arbitrary order
using the interpolation polynomial of Eq. 9. To recall the Runge-Kutta method,
consider the initial value problem

$\frac{dy}{dt}=f(t,y),\quad y(t_{m})\equiv y^{m}.$ (15)

The journey to a new step $y^{m}\to y^{m+1}$ with step-size $h$ is calculated
via a sequence of $n$ stages $\\{y^{m,c_{i}}\\}_{i=0}^{n-1}$ at the stage
times $\\{t_{i}\\}_{i=0}^{n-1}$ according to the system

$\displaystyle y^{m,c_{i}}$
$\displaystyle=y^{m}+h\sum_{j=0}^{n-1}A_{i}^{j}f(y^{m,c_{j}}),$ (16)
$\displaystyle y^{m+1}$
$\displaystyle=y^{m}+h\sum_{i=0}^{n-1}b_{i}f(y^{m,c_{i}})$ (17)

where $A_{i}^{j}$ is the RK coefficient matrix and $\\{b_{i}\\}_{i=0}^{n-1}$
the stage combination coefficients. The stage locations $\\{t_{i}\\}$,
combination coefficients $\\{b_{i}\\}$, and RK matrix $A_{i}^{j}$ are said to
form the Butcher tableau [3].

######  Theorem 2.

The Butcher tableau for the Gauss-Legendre IRK method consists of:

  * •

Stage times: the GL nodes $\\{\xi_{i}\\}_{i=0}^{n-1}$ transformed from
$[-1,1]\to[t^{m},t^{m+1}]$ according to

$\xi=\frac{h}{2}(t-\bar{t}^{m})$ (18)

with step-size $h\equiv t^{m+1}-t^{m}$ and mid-point
$\bar{t}^{m}\equiv\frac{1}{2}(t^{m+1}+t^{m})$,

  * •

Combination coefficients: the GL quadrature weig hts, $b_{i}\equiv w_{i}/2$
(with $\frac{1}{2}$ from the transform).

  * •

RK coefficient matrix: the matrix elements are given by

$A_{i}^{j}=\frac{w_{j}}{2}\Big{(}1+\sum_{s=0}^{n-1}P_{s}(\xi_{j})\frac{P_{s+1}(\xi_{i})-P_{s-1}(\xi_{i})}{2}\Big{)}$
(19)

where $P_{-1}(\xi)\equiv 1$, and by definition of the quadrature nodes
$P_{n}(\xi_{j})=0$.

Proof: Considering the ODE $\frac{dy}{dt}=f(t,y)$, project the RHS onto the
basis defined by the interpolants $\ell_{i}(t)$ through the GL quadrature
points $\\{t_{i}\\}_{i=1}^{n}$ defined by the affine transform of Eq. 18,

$f(y,t)\approx\sum_{j=0}^{n-1}f(y^{m,t_{j}})\ell_{j}(t).$ (20)

Solving the ODE by integration, one has

$\displaystyle y(t^{m,c_{i}})-y(t^{m})$
$\displaystyle=\int_{t^{m}}^{t^{m,c_{i}}}f(y,t)dt$ (21)
$\displaystyle\approx\sum_{j=0}^{n-1}f(y^{m,t_{j}})\int_{t_{m}}^{t_{m,c_{i}}}\ell_{j}(t)dt$
(22)

so the calculation reduces to integration of the interpolating polynomial.
Using Eq. 9 and the affine transform, this further reduces to integrating the
Legendre polynomial,

$\int_{t^{m}}^{t^{m,c_{i}}}\ell_{j}(t)dt=\frac{h}{2}w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{j})\int_{-1}^{\xi_{i}}P_{s}(\xi)d\xi.$
(23)

Recalling the identity

$P_{s}(\xi)=\frac{1}{2s+1}\frac{d}{d\xi}\Big{(}P_{s+1}(\xi)-P_{s-1}(\xi)\Big{)}$
(24)

the Legendre polynomial integrates to

$\int_{-1}^{\xi_{i}}P_{s}(\xi_{i})d\xi=\frac{1}{2}\frac{1}{s+\frac{1}{2}}\Big{[}(P_{s+1}(\xi_{i})-P_{s-1}(\xi_{i}))-(P_{s+1}(-1)-P_{s-1}(-1))\Big{]}.$
(25)

The boundary term $\xi=-1$ cancels for all $s>0$, but for $s=0$ an extra
constant factor is picked up. This establishes Eq. 19. Lastly, the combination
coefficients are found by integrating all the way to $t^{m+1}$, for an
integral $\int_{-1}^{1}P_{s}(\xi)d\xi=2\delta_{s0}$ by orthogonality. Thus,

$y^{m+1}=y^{m}+h\sum_{i=1}^{n}\frac{w_{i}}{2}f(y^{m,c_{i}}).$ (26)

The class of RK methods is very wide, and a large body of research exists on
developing general methods. One might get the impression that high-order
methods are difficult to construct. However, the preceding shows that when
Gauss-Legendre points are utilized the method may be developed to arbitrarily
high order, where only the quadrature nodes and weights are needed. For this
project methods up to $n=100$ are acquired using an online calculator [9],
though recent research has developed explicit formulas via asymptotic methods
to calculate the nodes and weights to floating point precision for any $n\geq
20$ [10] [11].

##  Appendix B The Newton-Raphson method for vector systems

The Newton-Raphson root-finding method is a beloved approximation and is close
to the hearts of many. To quickly review, given a system of equations in
vector form

$f^{\alpha}(\bm{x})=0$ (27)

its linearization about a point $x_{0}^{\alpha}$ is, using the Jacobian matrix
$J^{\alpha}_{\beta}\equiv\frac{\partial f^{\alpha}}{\partial x^{\beta}}$,

$f^{\alpha}(\bm{x}_{0})+J^{\alpha}_{\beta}(\bm{x}_{0})(x^{\beta}-x_{0}^{\beta})=0$
(28)

where summation convention is used. Based on an initial guess $x_{0}^{\beta}$,
one solves the linear system,

$J^{\alpha}_{\beta}(\bm{x}_{0})(x^{\beta}-x_{0}^{\beta})=-f^{\alpha}(\bm{x}_{0})$
(29)

for unknown $q^{\beta}\equiv x^{\beta}-x_{0}^{\beta}$. An iteration is then
obtained as

$x_{1}^{\beta}=q^{\beta}+x_{0}^{\beta}.$ (30)

One then repeats the process, sometimes adding a damping factor $|\gamma|\leq
1$ to the iteration to ensure convergence. Evidently, to use the method one
must find the Jacobian matrix of the system.

###  B.1 Newton-Raphson iteration for the implicit Runge-Kutta method

Considering a vector unknown $y^{\alpha}(t)$, the RK method integrates the
system of ODEs

$\frac{dy^{\alpha}}{dt}=f^{\alpha}(\bm{y},t).$ (31)

Denoting the right-hand side evaluation for the RK stages as the vector
$\bm{k}^{i}\equiv k^{\alpha,i}$, for a vector unknown the RK system may be
written as

$\displaystyle y^{m+1,\alpha}$
$\displaystyle=y^{m,\alpha}+hb_{i}k^{\alpha,i},$ (32) $\displaystyle
k^{i,\alpha}$
$\displaystyle=f^{\alpha}(t^{m}+c_{i}h,y^{m}+hA_{i}^{j}k^{\alpha}_{j})$ (33)

with RK matrix $A_{i}^{j}$ and coefficients $c_{i}$, $b_{i}$ found in Appendix
A.1. The second of these equations forms the implicit system of equations to
be solved, so it’s sufficient to apply Newton’s method to

$k^{i,\alpha}-f^{\alpha}(\bm{k}^{i})=0.$ (34)

Having linearized, for each iteration one solves the system

$\mathbb{J}^{i,\alpha}_{j,\beta}q^{j,\beta}=f^{\alpha}(\bm{k}^{i})-k^{i,\alpha}.$
(35)

The iteration update is then $k\texttt{ += }\gamma q$ with $\gamma$ an iterate
damping factor. The Jacobian $\mathbb{J}$ is a linear operator
$\mathbb{J}:\mathbb{R}^{\Lambda,n}\to\mathbb{R}^{\Lambda,n}$ for $\Lambda$
vector components and the $n$ RK stages. As a tensor-like object it is found
most easily in index notation.

######  Theorem 3.

The Jacobian tensor $\mathbb{J}^{i,\alpha}_{j,\beta}$ of the implicit Runge-
Kutta method is

$\mathbb{J}^{i,\alpha}_{j,\beta}=\delta^{i}_{j}\delta^{\alpha}_{\beta}-hA^{i}_{j}J^{\alpha,i}_{\beta}$
(36)

where $J^{\alpha,i}_{\beta}\equiv\partial_{y^{\beta}}f^{\alpha}(\bm{k}^{i})$
is the Jacobian matrix of the original ODE evaluated at the i’th RK stage.

Proof: Using the symbol $\partial_{k^{j,\beta}}$, one has

$\displaystyle\mathbb{J}^{i,\alpha}_{j,\beta}$
$\displaystyle=\partial_{k^{j,\beta}}(k^{i,\alpha}-f^{\alpha}(\bm{k}^{i}))$
(37)
$\displaystyle=\delta^{i}_{j}\delta^{\alpha}_{\beta}-\partial_{k^{j,\beta}}f^{\alpha}\Big{(}\bm{k}(\bm{y})\Big{)}$
(38)
$\displaystyle=\delta^{i}_{j}\delta^{\alpha}_{\beta}-(\partial_{k^{j,\beta}}y^{i,\gamma})\partial_{y^{i,\gamma}}f^{\alpha}(\bm{y})$
(39)
$\displaystyle=\delta^{i}_{j}\delta^{\alpha}_{\beta}-hA^{i}_{j}\delta_{\beta}^{\gamma}J^{\alpha,i}_{\gamma}$
(40)
$\displaystyle=\delta^{i}_{j}\delta^{\alpha}_{\beta}-hA^{i}_{j}J^{\alpha,i}_{\beta}$
(41)

In summary, the Jacobian for an IRK method of arbitrarily high order can be
calculated using: i) the RK coefficient matrix, and ii) the Jacobian of the
original ODE system. The products may be done numerically using
numpy.einsum(), for example. Of course, one needs a good guess for the RK
stages $y^{\alpha,i}$ to get started. A guess for large time-steps can be
furnished with a neural network.


11institutetext: Medical Imaging Center, Vingroup Big Data Institute, Hanoi,
Vietnam 22institutetext: School of Information and Communication Technology,
Hanoi University of Science and Technology, Hanoi, Vietnam 33institutetext:
College of Engineering & Computer Science, VinUniversity, Hanoi, Vietnam
44institutetext: Department of Mathematics, Yale University, New Heaven, USA
All correspondence should be addressed to Hieu H. Pham at 44email:
<EMAIL_ADDRESS>
# VinDr-SpineXR: A deep learning framework for spinal lesions detection and
classification from radiographs

Hieu T. Nguyen 1122 Hieu H. Pham 1133 Nghia T. Nguyen 11 Ha Q. Nguyen 1133
Thang Q. Huynh 22 Minh Dao 11 Van Vu 1144

###### Abstract

Radiographs are used as the most important imaging tool for identifying spine
anomalies in clinical practice. The evaluation of spinal bone lesions,
however, is a challenging task for radiologists. This work aims at developing
and evaluating a deep learning-based framework, named VinDr-SpineXR, for the
classification and localization of abnormalities from spine X-rays. First, we
build a large dataset, comprising 10,468 spine X-ray images from 5,000
studies, each of which is manually annotated by an experienced radiologist
with bounding boxes around abnormal findings in 13 categories. Using this
dataset, we then train a deep learning classifier to determine whether a spine
scan is abnormal and a detector to localize 7 crucial findings amongst the
total 13. The VinDr-SpineXR is evaluated on a test set of 2,078 images from
1,000 studies, which is kept separate from the training set. It demonstrates
an area under the receiver operating characteristic curve (AUROC) of 88.61%
(95% CI 87.19%, 90.02%) for the image-level classification task and a mean
average precision<EMAIL_ADDRESS>of 33.56% for the lesion-level localization task.
These results serve as a proof of concept and set a baseline for future
research in this direction. To encourage advances, the dataset, codes, and
trained deep learning models are made publicly available.

###### Keywords:

Spine X-rays Classification Detection Deep learning.

##  1 Introduction

###  1.1 Spine X-ray interpretation

Conventional radiography or X-ray has the ability to offer valuable
information than many other imaging modalities in the assessment of spinal
lesions [20, 29]. It has been the primary tool widely used to identify and
monitor various abnormalities of the spine. A wide range of spine conditions
can be observed from the X-rays like fractures, osteophytes, thinning of the
bones, vertebral collapse, or tumors [2, 24]. In clinical practice,
radiologists usually interpret and evaluate the spine on X-ray scans stored in
Digital Imaging and Communications in Medicine (DICOM) standard. Abnormal
findings could be identified based on differences in density, intensity, and
geometry of the lesions in comparison with normal areas. In many cases, those
differences might be subtle and it requires an in-depth understanding of
diagnostic radiography to spot them out. Large variability in the number,
size, and general appearance of spine lesions makes the interpretation of
spinal X-rays a complex and time-consuming task. These factors could lead to
the risk of missing significant findings [19], resulting in serious
consequences for patients and clinicians.

The rapid advances in machine learning, especially deep neural networks, have
demonstrated great potential in identifying diseases from medical imaging data
[26]. The integration of such systems in daily clinical routines could lead to
a much more efficient, accurate diagnosis and treatment [12]. In this study,
we aim to develop and validate a deep learning-based computer-aided diagnosis
(CAD) framework called VinDr-SpineXR that is able to classify and localize
abnormal findings from spine X-rays. Both the training and validation of our
proposed system are performed on our own dataset where radiologists’
annotations serve as a strong ground truth.

###  1.2 Related work

Developing CAD tools with a high clinical value to support radiologists in
interpreting musculoskeletal (MSK) X-ray has been intensively studied in
recent years [6]. Early approaches to the analysis of spine X-rays focus on
using radiographic textures for the detection of several specific pathologies
such as vertebral fractures [10], osteoporosis [9], or osteolysis [31].
Currently, deep convolutional networks [14] (CNNs) have shown their
significant improvements for the MSK analysis from X-rays [17, 25, 16, 11, 30,
33, 21]. Most of these studies focus on automated fracture detection and
localization [16, 11, 30, 33, 21]. To the best of our knowledge, no existing
studies devoted to the development and evaluation of a comprehensive system
for classifying and localizing multiple spine lesions from X-ray scans. The
lack of large datasets with high-quality images and human experts’ annotations
is the key obstacle. To fill this gap, this work focuses on creating a
significant benchmark dataset of spine X-rays that are manually annotated at
the lesion level by experienced radiologists. We also propose to develop and
evaluate, based on our dataset, a deep learning-based framework that includes
a normal versus abnormal classifier followed by a lesion detector that
localizes multiple categories of findings with bounding boxes. Both of these
tasks can be beneficial in clinical practice: the normal versus abnormal
classifier helps triage patients, while the lesion detector helps speed up the
reading procedure and complements radiologist’s observations.

###  1.3 Contribution

Our contributions in this paper are two folds. First, we present a new large-
scale dataset of 10,469 spine X-ray images from 5,000 studies that are
manually annotated with 13 types of abnormalities by radiologists. This is the
largest dataset to date that provides radiologist’s bounding-box annotations
for developing supervised-learning _object detection_ algorithms. Table 1
provides a summary of publicly available MSK datasets, in which two previous
spine datasets [5, 11] are significantly smaller than ours in size.
Furthermore, human expert’s annotations of spinal abnormalities are not
available in those datasets. Second, we develop and evaluate VinDr-SpineXR – a
deep learning framework that is able to classify and localize multiple spine
lesions. Our main goal is to provide a nontrivial baseline performance of
state-of-the-art deep learning approaches on the released dataset, which could
be useful for further model development and comparison of novel CAD tools. To
facilitate new advances, we have made the dataset available for public access
on our project’s webpage111https://vindr.ai/datasets/spinexr. The codes used
in the experiments are available at Github222https://github.com/vinbigdata-
medical/vindr-spinexr.

Table 1: Overview of publicly available MSK image datasets.  Dataset | Year | Study type | Label | # Im.
---|---|---|---|---
Digital Hand Atlas [5] | 2007 | Left hand | Bone age | 1,390
Osteoarthritis Initiative [1] | 2013 | Knee | K&L Grade | 8,892
MURA [21] | 2017 | Upper body | Abnormalities | 40,561
RSNA Pediatric Bone Age [7] | 2019 | Hand | Bone age | 14,236
Kuok et al. [13] | 2018 | Spine | Lumbar vertebrae mask | 60
Kim et al. [11] | 2020 | Spine | Spine position | 797
Ours | 2021 | Spine | Multiple abnormalities | 10,469

##  2 Proposed Method

###  2.1 Overall framework

The proposed deep learning framework (see Figure 1) includes two main parts:
(1) a classification network, which accepts a spine X-ray as input and
predicts if it could be a normal or abnormal scan; (2) a detection network
receives the same input and predicts the location of abnormal findings. To
maximize the framework’s detection performance, we propose a decision rule to
combine the outputs of the two networks.

Figure 1: Overview of the VinDr-SpineXR for spine abnormalities classification
and localization. A binary classifier takes as input one spine scan and
predicts its probability of abnormality. A detector takes the same scan as
input and provides bounding boxes along with probabilities of abnormal
findings at lesion-level. A decision rule is then proposed to combine the two
outputs and maximize the detection performance.

###  2.2 Dataset

####  2.2.1 Data collection

Spine X-rays with lesion-level annotations are needed to develop automated
lesion detection systems. In this work, more than 50,000 raw spine images in
DICOM format were retrospectively collected from the Picture Archive and
Communication System (PACS) of different primary hospitals between 2010–2020.
The data collection process was conducted under our cooperation with
participating hospitals. Since this research did not impact clinical care,
patient consent was waived. To keep patient’s Protected Health Information
(PHI) secure, all patient-identifiable information associated with the images
has been removed. Several DICOM attributes that are important for evaluating
the spine conditions like patient’s age and sex were retained.

####  2.2.2 Data annotation

To annotate data, we developed an in-house labeling framework called VinDr Lab
[18] which was built on top of a PACS. This is a web-based framework that
allows multiple radiologists to work remotely at the same time. Additionally,
it also provides a comprehensive set of annotation tools that help maximize
the performance of human annotators. A set of 5,000 spine studies were
randomly selected from the raw data after removing outliers (e.g. scans of
other body parts or low-quality). All these scans were then uploaded to the
labeling framework and assigned to 3 participating radiologists, who have at
least 10 years of experience, such that each scan is annotated by exactly one
radiologist. In this step, the radiologists decide whether a scan was abnormal
or normal as well as marked the exact location of each abnormality on the
scan. The labels and annotations provided by the radiologists were served as
ground truth for model development and validation later. During this process,
the radiologists were blinded to relevant clinical information except
patient’s age and sex. Finally, a total of 10,468 spine images from 5,000
studies were annotated for the presence of 13 abnormal findings: ankylosis,
disc space narrowing, enthesophytes, foraminal stenosis, fracture,
osteophytes, sclerotic lesion, spondylolysthesis, subchondral sclerosis,
surgical implant, vertebral collapse, foreign body, and other lesions. The “no
finding” label was intended to represent the absence of all abnormalities. For
the development of the deep learning algorithms, we randomly stratified the
labeled dataset, at study-level, into a development set of 4,000 studies and a
test set of 1,000 studies. Table 2 summarizes the data characteristics,
including patient demographic and the prevalence of each label via the number
of bounding boxes. Figure 2 shows several representative samples with and
without abnormal findings from our dataset.

Table 2: Characteristics of patients in the training and test datasets.

|  Characteristic |  Training set |  Test set |  Total
---|---|---|---|---
Statistics |  Years |  2011 to 2020 |  2015 to 2020 |  2011 to 2020
Number of studies |  4000 |  1000 |  5000
Number of images |  8390 |  2078 |  10468
Number of normal images |  4257 |  1068 |  5325
Number of abnormal images |  4133 |  1010 |  5143
Image size (pixel$\times$pixel, mean) |  2455 $\times$ 1606 |  2439 $\times$ 1587 |  2542 $\times$ 1602
Age (mean, years [range])* |  49 [6 – 94] |  49 [13 – 98] |  49 [6 – 98]
Male (%)* |  37.91 |  37.53 |  37.83
Female (%)* |  62.08 |  62.47 |  62.17
Data size (GiB) |  29.04 |  7.17 |  36.21
Lesion type |  1\. Ankylosis |  6 |  1 |  7
2\. Disc space narrowing |  924 |  231 |  1155
3\. Enthesophytes |  54 |  13 |  67
4\. Foraminal stenosis |  387 |  95 |  482
5\. Fracture |  10 |  2 |  12
6\. Osteophytes |  12562 |  3000 |  15562
7\. Sclerotic lesion |  24 |  5 |  29
8\. Spondylolysthesis |  280 |  69 |  349
9\. Subchondral sclerosis |  87 |  23 |  110
10\. Surgical implant |  429 |  107 |  536
11\. Vertebral collapse |  268 |  69 |  337
12\. Foreign body |  17 |  4 |  21
13\. Other lesions |  248 |  63 |  311

(*) These calculations were performed on studies where sex and age were
available.

Figure 2: Examples of spine X-rays with radiologist’s annotations, in which
abnormal findings are marked by rectangular bounding boxes.

###  2.3 Model development

####  2.3.1 Network architecture and training methodology

To classify a spine X-ray image as either normal or abnormal, three CNNs
(DenseNet-121, DenseNet-169, DenseNet-201) have been deployed. These networks
[8] are well-known to be effective for X-ray interpretation [22, 21]. Each
network accepts an image of the spine as input and outputs a binary label
corresponding to normal or abnormal scan. A total of 4,257 spine images with
normal findings and 4,133 spine images with abnormal findings (reflecting any
of the lesions) from the training set was used to optimize networks’ weights.
During this stage, we optimized cross-entropy loss between the image-level
labels and network’s outputs using SGD optimizer. The average ensemble of
three models serves as the final classifier. An image is considered abnormal
in the inference stage if its probability of abnormality is greater than an
optimal threshold. In particular, we determine the optimal threshold $c^{*}$
for the classifier by maximizing Youden’s index [32], $J(c)=q(c)+r(c)-1$,
where the sensitivity $q$ and the specificity $r$ are functions of the cutoff
value $c$. For the detection task, we aim to localize 7 important lesions:
osteophytes, disc space narrowing, surgical implant, foraminal stenosis,
spondylolysthesis, vertebral collapse, and other lesions. Due to the limited
number of positive samples, the rest lesions were considered “other lesions”.
State-of-the-art detectors, namely Faster R-CNN [23], RetinaNet [15],
EfficientDet [28], and Sparse R-CNN [27], have been deployed for this task.
Faster R-CNN [23] was chosen as a representative for anchor-based two-stage
detectors, which firstly propose a set of candidate regions then categorize
and refine their locations. RetinaNet [15] and EfficientDet [28] are both one-
state detectors that directly localize and classify the densely proposed
anchor boxes. Different from previous detectors, Sparse R-CNN [27] starts from
a set of initial proposals, then repeatedly refines and classifies these boxes
using attention mechanism. All the networks were trained to localize spine
lesions using the stochastic gradient descent (SGD) optimizer. During the
learning phase, the bounding box regression loss and region-level
classification loss were jointly minimized. To improve the generalization
performance of the detectors, learned data augmentation strategies for object
detection [34] were incorporated in the training procedure.

####  2.3.2 Decision fusion rule

Given an input image x, we denote $\hat{p}(\verb|abnormal||\textbf{x})$ as the
classifier’s output that reflects the probability of the image being abnormal.
To maximize the performance of the lesion detector, we propose a heuristic
fusion rule as follows. For any x with the prediction
$\hat{p}(\verb|abnormal||\textbf{x})\geq c^{*}$, all lesion detection results
are retained. For the case $\hat{p}(\verb|abnormal||\textbf{x})<c^{*}$, only
predicted bounding boxes with confidence higher than $0.5$ are kept.

##  3 Experiments and Results

###  3.1 Experimental setup & Implementation details

Several experiments were conducted to evaluate the performance of the VinDr-
SpineXR. First, we evaluated the classification and detection networks
independently on the test set. We then investigated the effect of the fusion
rule on the detection performance of the whole framework. All networks were
implemented and trained using PyTorch framework (version 1.7.1) on a server
with one NVIDIA V100 32GiB GPU. For the classification task, all training
images were rescaled to 224 $\times$ 224 pixels and normalized by the mean and
standard deviation of images from the ImageNet dataset. The networks were
initialized with pre-trained weights on ImageNet. Each network was trained
end-to-end using SGD for 10,000 iterations, approximately 2 hours. We used
mini-batches of size 32 and set the learning rate to 0.001. An ensemble of the
three best models served as the final classifier. For the detection task, all
training images were randomly downsampled such that the shorter edge ranging
between 640 and 800 pixels, then random augmentation transforms were applied
before grouping into mini-batches of 16 samples. The detectors were
initialized with weights pre-trained on COCO dataset then trained for 50,000
iterations (about 24 hours) using SGD with the learning rate reduced by
$10\times$ at the 30,000-th iteration.

###  3.2 Evaluation metrics

We report the classification performance using AUROC, sensitivity,
specificity, and F1 score. We also estimate the 95% confidence interval (CI)
by bootstrapping with 10,000 bootstraps for each measure. In each bootstrap, a
new dataset is constructed by sampling with replacement [3]. For the detection
task, we followed PASCAL VOC metric, mean average precision<EMAIL_ADDRESS>[4]. A
predicted finding is a true positive if it has an intersection over union
(IoU) of at least 0.5 with a ground truth lesion of the same class. For each
class, average precision (AP) is the mean of 101 precision values,
corresponding to recall values ranging from 0 to 1 with a step size of 0.01.
The final metric<EMAIL_ADDRESS>is the mean of AP over all lesion categories.

###  3.3 Experimental results

####  3.3.1 Classification performance

On the test set of 2,078 images, our ensemble model reported an AUROC of
88.61% (95% CI 87.19%, 90.02%). At the optimal operating point of 0.2808 (see
Supplementary Material), the network achieved a sensitivity of 83.07% (95% CI
80.64%, 85.34%), a specificity of 79.32% (95% CI 77.84%, 81.71%), and an F1
score of 81.06% (95% CI 79.20%, 82.85%), respectively. Table 3 shows
quantitative results over all the classification networks.

Table 3: Classification performance on the test set (in percent with 95% CI). Classifier |  AUROC |  F1 score |  Sensitivity |  Specificity
---|---|---|---|---
DenseNet121 |  86.93 (85.4,88.4) |  79.55 (77.6,81.4) |  80.39 (77.9,82.8) |  79.32 (76.9,81.7)
DenseNet169 |  87.29 (85.8,88.8) |  80.25 (78.3,82.1) |  81.74 (79.3,84.1) |  79.02 (76.6,81.4)
DenseNet201 |  87.14 (85.6,88.6) |  79.03 (77.1,80.9) |  77.97 (75.4,80.5) |  81.46 (79.1,83.8)
Ensemble |  88.61 (87.2,90.0) |  81.06 (79.2,82.9) |  83.07 (80.6,85.3) |  79.32 (77.8,81.7)

####  3.3.2 Detection performance

The AP of each abnormal finding detected by 4 detectors is shown in Table 4.
Sparse R-CNN [27] showed it as the best performing detector for this task with
a<EMAIL_ADDRESS>of 33.15% (see Supplementary Material). Meanwhile, RetinaNet [15]
showed the worst performance with a<EMAIL_ADDRESS>of 28.09%. We observed that the
reported performances varied over the target lesions, e.g. all the detectors
performed best on the LT10 label (surgical implant) and worst on the LT13
(other lesions).

Table 4: Spine X-ray detection performance on the test set.

Detector |  LT2(∗) |  LT4 |  LT6 |  LT8 |  LT10 |  LT11 |  LT13 |<EMAIL_ADDRESS>
---|---|---|---|---|---|---|---|---
Faster R-CNN [23] |  22.66 |  35.99 |  49.24 |  31.68 |  65.22 |  51.68 |  2.16 |  31.83
RetinaNet [15] |  14.53 |  25.35 |  41.67 |  32.14 |  65.49 |  51.85 |  5.30 |  28.09
EfficientDet [28] |  17.05 |  24.19 |  42.69 |  35.18 |  61.85 |  52.53 |  2.45 |  28.73
Sparse R-CNN [27] |  20.09 |  32.67 |  48.16 |  45.32 |  72.20 |  49.30 |  5.41 |  33.15

(*) LT2, LT4, LT6, LT8, LT10, LT11, LT13 denotes for disc space narrowing,
foraminal stenosis, osteophytes, spondylolysthesis, surgical implant,
vertebral collapse and other lesions, respectively, following the same
indexing in Table 2.

####  3.3.3 Effect of the decision fusion rule

Table 5 provides a comparison of the detection performance between the
detector network and the whole framework. We observed that by combining the
detection and classification networks via the proposed decision rule, the
whole framework leads to a higher performance over all lesions, except LT4 and
LT8. Visualization of predicted lesions by the VinDr-SpineXR is provided in
the Supplementary Material.

Table 5: Spine X-ray detection performance of the whole framework on the test set. Method |  LT2 |  LT4 |  LT6 |  LT8 |  LT10 |  LT11 |  LT13 |<EMAIL_ADDRESS>
---|---|---|---|---|---|---|---|---
Detector only |  18.06 |  28.89 |  34.23 |  41.52 |  62.45 |  42.85 |  4.03 |  33.15
Whole framework |  21.43 |  27.36 |  34.78 |  41.29 |  62.53 |  43.39 |  4.16 |  33.56

##  4 Discussion

The proposed decision fusion uses the result of the classifier to influence
the detector. As shown in Table 5, this rule, although very simple, helps
improve the mAP of the detector by 0.41%. We have also experimented with a
counterpart fusion rule to use the result of the detector to influence the
classifier. Specifically, we averaged classifier’s output with the highest box
score from detector’s outputs, boosting the AUROC and F1 score of the
classifier by 1.58% and 0.46%, respectively. These experiments have
highlighted the effectiveness of the proposed mutual ensemble mechanism.

##  5 Conclusions

In this work, we developed and validated VinDr-SpineXR – a deep learning-based
framework to detect and localize abnormalities in radiographs of spine. We
contributed a new labeled dataset and models for spine X-ray analysis. To the
best of our knowledge, this is the first effort to address the problem of
multiple lesions detection and localization in spine X-rays. Our experiments
on the dataset demonstrated the effectiveness of the proposed method. For
future work, we expect to extend our dataset and consider more lesion labels
for further experiments, including rare labels. We also plan to conduct more
experiments and evaluate the impact of the proposed framework in real-world
clinical scenarios.

## References

  * [1] Osteoarthritis initiative: A multi-center observational study of men and women. https://oai.epi-ucsf.org/datarelease/, accessed: 2021-02-22
  * [2] Baert, A.: Spinal imaging: Diagnostic imaging of the spine and spinal cord. Springer Science & Business Media (2007)
  * [3] Efron, B., Tibshirani, R.J.: An introduction to the bootstrap. No. 57 in Monographs on Statistics and Applied Probability, Chapman & Hall/CRC (1993)
  * [4] Everingham, M., Van Gool, L., Williams, C.K., Winn, J., Zisserman, A.: The Pascal visual object classes (VOC) challenge. International Journal of Computer Vision 88(2), 303–338 (2010)
  * [5] Gertych, A., Zhang, A., Sayre, J., Pospiech-Kurkowska, S., Huang, H.: Bone age assessment of children using a digital hand atlas. Computerized Medical Imaging and Graphics 31(4-5), 322–331 (2007)
  * [6] Gundry, M., Knapp, K., Meertens, R., Meakin, J.: Computer-aided detection in musculoskeletal projection radiography: A systematic review. Radiography 24(2), 165–174 (2018)
  * [7] Halabi, S.S., Prevedello, L.M., Kalpathy-Cramer, J., Mamonov, A.B., Bilbily, A., Cicero, M., Pan, I., Pereira, L.A., Sousa, R.T., Abdala, N., et al.: The RSNA pediatric bone age machine learning challenge. Radiology 290(2), 498–503 (2019)
  * [8] Huang, G., Liu, Z., Van Der Maaten, L., Weinberger, K.Q.: Densely connected convolutional networks. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp. 4700–4708 (2017)
  * [9] Kasai, S., Li, F., Shiraishi, J., Li, Q., Doi, K.: Computerized detection of vertebral compression fractures on lateral chest radiographs: Preliminary results with a tool for early detection of osteoporosis. Medical Physics 33(12), 4664–4674 (2006)
  * [10] Kasai, S., Li, F., Shiraishi, J., Li, Q., Nie, Y., Doi, K.: Development of computerized method for detection of vertebral fractures on lateral chest radiographs. In: Medical Imaging: Image Processing. vol. 6144, p. 61445D (2006)
  * [11] Kim, K.C., Cho, H.C., Jang, T.J., Choi, J.M., Seo, J.K.: Automatic detection and segmentation of lumbar vertebrae from x-ray images for compression fracture evaluation. Computer Methods and Programs in Biomedicine p. 105833 (2020)
  * [12] Krupinski, E.A., Berbaum, K.S., Caldwell, R.T., Schartz, K.M., Kim, J.: Long radiology workdays reduce detection and accommodation accuracy. Journal of the American College of Radiology 7(9), 698–704 (2010)
  * [13] Kuok, C.P., Fu, M.J., Lin, C.J., Horng, M.H., Sun, Y.N.: Vertebrae segmentation from X-ray images using convolutional neural network. In: International Conference on Information Hiding and Image Processing (IHIP). pp. 57–61 (2018)
  * [14] LeCun, Y., Bengio, Y., Hinton, G.: Deep learning. Nature 521(7553), 436–444 (2015)
  * [15] Lin, T.Y., Goyal, P., Girshick, R., He, K., Dollár, P.: Focal loss for dense object detection. In: IEEE International Conference on Computer Vision (ICCV). pp. 2980–2988 (2017)
  * [16] Lindsey, R., Daluiski, A., Chopra, S., Lachapelle, A., Mozer, M., Sicular, S., Hanel, D., Gardner, M., Gupta, A., Hotchkiss, R., et al.: Deep neural network improves fracture detection by clinicians. Proceedings of the National Academy of Sciences 115(45), 11591–11596 (2018)
  * [17] Mandal, I.: Developing new machine learning ensembles for quality spine diagnosis. Knowledge-based Systems 73, 298–310 (2015)
  * [18] Nguyen, N.T., Truong, P.T., Ho, V.T., Nguyen, T.V., Pham, H.T., Nguyen, M.T., Dam, L.T., Nguyen, H.Q.: VinDr Lab: A Data Platform for Medical AI. https://github.com/vinbigdata-medical/vindr-lab (2021)
  * [19] Pinto, A., Berritto, D., Russo, A., Riccitiello, F., Caruso, M., Belfiore, M.P., Papapietro, V.R., Carotti, M., Pinto, F., Giovagnoni, A., et al.: Traumatic fractures in adults: Missed diagnosis on plain radiographs in the emergency department. Acta Bio Medica: Atenei Parmensis 89(1), 111 (2018)
  * [20] Priolo, F., Cerase, A.: The current role of radiography in the assessment of skeletal tumors and tumor-like lesions. European Journal of Radiology 27, S77–S85 (1998)
  * [21] Rajpurkar, P., Irvin, J., Bagul, A., Ding, D., Duan, T., Mehta, H., Yang, B., Zhu, K., Laird, D., Ball, R.L., et al.: MURA: Large dataset for abnormality detection in musculoskeletal radiographs. arXiv preprint arXiv:1712.06957 (2017)
  * [22] Rajpurkar, P., Irvin, J., Zhu, K., Yang, B., Mehta, H., Duan, T., Ding, D., Bagul, A., Langlotz, C., Shpanskaya, K., et al.: ChexNet: Radiologist-level pneumonia detection on chest X-rays with deep learning. arXiv preprint arXiv:1711.05225 (2017)
  * [23] Ren, S., He, K., Girshick, R., Sun, J.: Faster r-cnn: Towards real-time object detection with region proposal networks. In: Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 28. Curran Associates, Inc. (2015), https://proceedings.neurips.cc/paper/2015/file/14bfa6bb14875e45bba028a21ed38046-Paper.pdf
  * [24] Rodallec, M.H., Feydy, A., Larousserie, F., Anract, P., Campagna, R., Babinet, A., Zins, M., Drapé, J.L.: Diagnostic imaging of solitary tumors of the spine: What to do and say. Radiographics 28(4), 1019–1041 (2008)
  * [25] Sa, R., Owens, W., Wiegand, R., Studin, M., Capoferri, D., Barooha, K., Greaux, A., Rattray, R., Hutton, A., Cintineo, J., et al.: Intervertebral disc detection in X-ray images using Faster R-CNN. In: Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). pp. 564–567 (2017)
  * [26] Shen, D., Wu, G., Suk, H.I.: Deep learning in medical image analysis. Annual Review of Biomedical Engineering 19, 221–248 (2017)
  * [27] Sun, P., Zhang, R., Jiang, Y., Kong, T., Xu, C., Zhan, W., Tomizuka, M., Li, L., Yuan, Z., Wang, C., Luo, P.: Sparse R-CNN: End-to-end object detection with learnable proposals. arXiv preprint arXiv:2011.12450 (2020)
  * [28] Tan, M., Pang, R., Le, Q.V.: EfficientDet: Scalable and efficient object detection. In: IEEE Conference on Computer Vision and Pattern Recognition (ECCV). pp. 10781–10790 (2020)
  * [29] Tang, C., Aggarwal, R.: Imaging for musculoskeletal problems. InnovAiT 6(11), 735–738 (2013)
  * [30] Thian, Y.L., Li, Y., Jagmohan, P., Sia, D., Chan, V.E.Y., Tan, R.T.: Convolutional neural networks for automated fracture detection and localization on wrist radiographs. Radiology: Artificial Intelligence 1(1), e180001 (2019)
  * [31] Wilkie, J.R., Giger, M.L., Pesce, L.L., Engh Sr, C.A., Hopper Jr, R.H., Martell, J.M.: Imputation methods for temporal radiographic texture analysis in the detection of periprosthetic osteolysis. In: Medical Imaging: Computer-Aided Diagnosis. vol. 6514, p. 65141L (2007)
  * [32] Youden, W.J.: Index for rating diagnostic tests. Cancer 3(1), 32–35 (1950)
  * [33] Zhang, X., Wang, Y., Cheng, C.T., Lu, L., Xiao, J., Liao, C.H., Miao, S.: A new window loss function for bone fracture detection and localization in X-ray images with point-based annotation. arXiv preprint arXiv:2012.04066 (2020)
  * [34] Zoph, B., Cubuk, E.D., Ghiasi, G., Lin, T.Y., Shlens, J., Le, Q.V.: Learning data augmentation strategies for object detection. In: IEEE European Conference on Computer Vision (ECCV). pp. 566–583 (2020)