Technical report

Neutron Transmission Strain Tomography for Non-Constant Stress-Free Lattice
Spacing

J.N. Hendriks, C. Jidling, T.B. Schön, A. Wills, C.M. Wensrich, E.H. Kisi

  * $\bullet$

Please cite this version:
J.N. Hendriks, C. Jidling, T.B. Schön, A. Wills, C.M. Wensrich, E.H. Kisi.
Neutron Transmission Strain Tomography for Non-Constant Stress-Free Lattice
Spacing. Nuclear instruments and methods in physics research section B,
456:64-73, 2019.

###### Abstract

Recently, several algorithms for strain tomography from energy-resolved
neutron transmission measurements have been proposed. These methods assume
that the stress-free lattice spacing $d_{0}$ is a known constant limiting
their application to the study of stresses generated by manufacturing and
loading methods that do not alter this parameter. In this paper, we consider
the more general problem of jointly reconstructing the strain and $d_{0}$
fields. A method for solving this inherently non-linear problem is presented
that ensures the estimated strain field satisfies equilibrium and can include
knowledge of boundary conditions. This method is tested on a simulated data
set with realistic noise levels, demonstrating that it is possible to jointly
reconstruct $d_{0}$ and the strain field.

##  1 Introduction

Energy-resolved neutron transmission methods can generate lower dimensional
(one- or two-dimensional) images of strain from a higher dimensional (two- or
three-dimensional) strain field within a polycrystalline material. The
‘tomographic’ reconstruction of an unknown strain field from these images can
be used to study the residual strain and stress within engineering components.
Residual stresses are those which remain after applied loads are removed (e.g.
due to heat treatment, plastic deformation, etc.), and may have significant
and unintended impact on a component’s effective strength and service life —
in particular its fatigue life. Measuring and quantifying these strains is
important for the validation of predictive design tools, such as Finite
Element Analysis, and to aid the development of novel manufacturing techniques
— i.e. additive manufacturing.

These strain images are generated by analysing features known as Bragg-edges
in the relative transmission of a neutron pulse through a sample. Bragg-edges
are sudden increases in the intensity as a function of wavelength and occur
when the scattering angle $2\vartheta$ reaches $180^{\circ}$, beyond which no
further coherent scattering can occur. The wavelength $\lambda$ at which these
Bragg-edges occur can be related to the lattice spacing $d$ within the sample
through Bragg’s law: $\lambda=2d\sin\vartheta$. Assuming minimal texture, this
can be used to provide a relative measure of strain;

$\langle\epsilon\rangle=\frac{d-d_{0}}{d_{0}},$ (1)

where $d_{0}$ is the stress-free lattice spacing and $\langle\epsilon\rangle$
is a through thickness average of the normal, elastic strain in the direction
of the beam.

The determination of $d_{0}$ is a problem inherent to diffraction and
transmission strain analysis. For specific cases where the loading mechanism
does not result in changes to the stress-free lattice parameter, its value can
be measured prior to loading and in the simplest case (e.g. for an annealed
sample) a constant value throughout the sample can be assumed. Several
algorithms for strain tomography assuming a known, constant stress-free
lattice spacing have been developed. Reconstruction of axisymmetric strain
fields is considered in (Abbey et al., 2009, 2012; Kirkwood et al., 2015;
Gregg et al., 2017) and more general two-dimensional strain fields in (Gregg
et al., 2018; Jidling et al., 2018; Hendriks et al., 2018a).

Many manufacturing techniques (e.g. welding and additive manufacturing) can
alter the lattice spacing; for example, as a result of inhomogeneously
distributed phase changes (such as the Martensite transformation), or due to
gradients in composition as a result of differing chemical states in the
starting materials. Since the lattice spacing (in this case $d_{0}$) is
sensitive to crystal structure and composition changes, the stress-free
lattice parameter may vary throughout the sample. Ignoring variations in
$d_{0}$ would cause severe degradation in the quality of a reconstructed
strain field. In such cases, measuring $d_{0}$ is more challenging and has
been achieved in neutron diffraction measurements by measuring additional
directions of strain (Luzin et al., 2011; Choi et al., 2007) and by
destructive methods where the strain is relieved by wire cutting the sample
into a grid allowing the stress-free lattice spacing to be measured throughout
the sample (Paradowska et al., 2005). Although the latter of these two options
could be applied to strain tomography it requires the destruction of the
sample and creates an additional tomography problem, requiring another set of
measurements to be acquired.

Here, we present a method capable of jointly reconstructing the strain field
and the $d_{0}$ field from a single set of neutron transmission images. To
achieve this both the strain and $d_{0}$ are modelled by a Gaussian process
(see for example Rasmussen and Williams (2006)) and equilibrium and boundary
conditions are built into the strain model (Jidling et al., 2017). This
extends the Gaussian process approach presented by Jidling et al. (2018);
Hendriks et al. (2018a) to handle the inherently non-linear nature of this
problem. A numerically tractable algorithm based on variational inference (see
for example Blei et al. (2017); Jordan et al. (1999)) is provided and the
method is validated on a simulated data set.

##  2 Problem Statement

This paper focuses on the joint reconstruction of the strain field
$\boldsymbol{\epsilon}(\mathbf{x})$ and a non-constant stress-free lattice
parameter $d_{0}(\mathbf{x})$ from a set of neutron transmission images.
Restricting the problem to two dimensions, gives the strain field as the
symmetric tensor

$\boldsymbol{\epsilon}(\mathbf{x})=\begin{bmatrix}\epsilon_{xx}(\mathbf{x})&\epsilon_{xy}(\mathbf{x})\\\
\epsilon_{xy}(\mathbf{x})&\epsilon_{yy}(\mathbf{x})\end{bmatrix},$ (2)

where $\mathbf{x}=\begin{bmatrix}x&y\end{bmatrix}^{\mathsf{T}}$. For brevity,
the unique components of strain will be written as
$\bar{\epsilon}=\begin{bmatrix}\epsilon_{xx}&\epsilon_{xy}&\epsilon_{yy}\end{bmatrix}^{\mathsf{T}}$
with the coordinate $\mathbf{x}$ omitted where appropriate.

Here, we consider the lattice spacings $d$ as the measurements rather than the
standard approach which considers the relative strain of the form form (1).
This allows the measurements to be explicitly related to both the strain and
the stress-free lattice parameter through the Longitudinal Ray Transform (LRT)
(Lionheart and Withers, 2015):

$y(\boldsymbol{\eta})=d(\boldsymbol{\eta})+e=\frac{1}{L}\int\limits_{0}^{L}\bar{\mathbf{n}}\bar{\boldsymbol{\epsilon}}(\mathbf{p}+\hat{\mathbf{n}}s)d_{0}(\mathbf{p}+\hat{\mathbf{n}}s)+d_{0}(\mathbf{p}+\hat{\mathbf{n}}s)\,\mathrm{d}s+e.$
(3)

where $e\sim\mathcal{N}(0,\sigma_{n}^{2})$ and the geometry of each
measurement is given by the parameter set
$\boldsymbol{\eta}=\\{\hat{\mathbf{n}},L,\mathbf{p}\\}$; with
$\hat{\mathbf{n}}$ as the beam direction, $L$ as the irradiation length,
$\mathbf{p}=\begin{bmatrix}x_{0}&y_{0}\end{bmatrix}^{\mathsf{T}}$ as the point
of initial intersection between the ray and the sample, and
$\bar{\mathbf{n}}=\begin{bmatrix}\hat{n}_{1}^{2}&2\hat{n}_{1}\hat{n}_{2}&\hat{n}_{2}^{2}\end{bmatrix}$.
See Figure 1 for the measurement geometry. These measurements are a non-linear
function of the two phenomena we wish to estimate; $\boldsymbol{\epsilon}$ and
$d_{0}$.

For details on the analysis of neutron transmission data to determine these
lattice spacings the reader is referred to Santisteban et al. (2002b, a);
Tremsin et al. (2011, 2012). It is also worth noting that the standard
deviation $\sigma_{n}$ of these measurements is available.

Figure 1: LRT measurement geometry. Each measurement made by a detector pixel
is associated with a ray of direction $\hat{\mathbf{n}}$ that enters the
sample at $\mathbf{p}$ and has a total irradiated length of $L$.

Furthermore, the strain field inside a sample is a physical property and as
such it is subject to equilibrium and boundary conditions (Sadd, 2009).
Therefore, it is natural to constrain estimates of the strain field to satisfy
these conditions. Using Hooke’s law the equilibrium conditions can be written
directly in terms of strain. In two dimensions, this relies on an assumption
of plane strain or plane stress. Plane stress is assumed for the remainder of
this work, giving the equilibrium conditions as

$\begin{split}\frac{\partial}{\partial
x}(\epsilon_{xx}+\nu\epsilon_{yy})+\frac{\partial}{\partial
y}(1-\nu)\epsilon_{xy}&=0,\\\ \frac{\partial}{\partial
x}(\epsilon_{yy}+\nu\epsilon_{xx})+\frac{\partial}{\partial
y}(1-\nu)\epsilon_{xy}&=0,\\\ \end{split}$ (4)

where $\nu$ is Poisson’s ratio.

Boundary conditions, in particular the load free surfaces, may also be known.
For an unloaded surface, the distribution of forces known as tractions will be
zero. Through equilibrium this places additional linear constraints on the
strain field, which, assuming plane stress, can be written as

$\mathbf{0}=\begin{bmatrix}n_{\perp 1}&n_{\perp 2}&0\\\ 0&n_{\perp 1}&n_{\perp
2}\end{bmatrix}\begin{bmatrix}1&0&-\nu\\\ 0&1+\nu&0\\\
-\nu&0&1\end{bmatrix}\bar{\boldsymbol{\epsilon}}(\mathbf{x}_{b}),$ (5)

where $\mathbf{x}_{b}$ is a point on an unloaded surface and
$\mathbf{n}_{\perp}$ is the normal to the surface at this point.

An approach to enforcing equilibrium in the estimated strain field is to
define a Gaussian process for the Airys stress function from which strain can
be derived (Jidling et al., 2018). This non-parametric approach was
demonstrated experimentally by Jidling et al. (2018) and compared to other
parametric approaches by Hendriks et al. (2018a) with promising results.
Boundary conditions in the form of (5) can be included in the estimation
process as artificial measurements of zero traction (Hendriks et al., 2018a).

We wish to extend this approach so that both the strain field and the stress-
free lattice spacing can be estimated. As the measurements are a non-linear
function of the unknowns we cannot directly apply the standard Gaussian
process regression methods (Rasmussen and Williams, 2006). There exists
several approaches to approximate Gaussian processes for non-linear functions;
the Laplace approximation (Rasmussen and Williams, 2006; Bishop et al., 1995),
GP variational inference (Steinberg and Bonilla, 2014), and Markov Chain Monte
Carlo methods (such as Elliptical Slice Sampling (Murray and Adams, 2010)).
For these methods, the measurements are modelled as non-linear functions of
the GP sampled at the measurement locations (known as latent function values).
The latent function values that best111For a given criterion of best fit,
whether it be marginal log likelihood, cross-validation, etc. match the data
are determined by one of the above methods. Then, Gaussian process regression
is applied with the latent function values taking the place of measurements to
determine the function values at the new locations of interest.

The non-static nature of the integral measurement model (3) makes it unclear
how to express the measurements as a function of a finite set of latent
function values, and hence the above approaches to approximating the GP for
non-linear measurements cannot be applied directly. In the following section,
we utilise an finite basis function approximation to the GP, and by viewing
the problem from an alternate perspective we show how variational inference
can be used to solve this non-linear problem.

##  3 Method

The method presented here is to define a Gaussian process model for the strain
field and the stress-free lattice spacing. This Gaussian process model is then
approximated using a Hilbert space approximation (Solin and Särkkä, 2014;
Jidling et al., 2018). This has two benefits; firstly it removes the need for
numerical integration of the covariance function (as discussed by Jidling et
al. (2018, 2017)), and secondly it allows us to reformulate the problem as a
set of basis functions with unknown coefficients. Variational inference can
then be used to learn the coefficients from the LRT measurements and
artificial measurements of zero traction.

###  3.1 Gaussian Process model

The Gaussian process (GP) is a Gaussian distribution of spatially correlated
functions;

$f(\mathbf{x})\sim\mathcal{GP}\left(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime})\right).$
(6)

The characteristics of the functions belonging to this distribution are
governed by a mean function $m(\mathbf{x})$ and a covariance function
$k(\mathbf{x},\mathbf{x}^{\prime})$. The covariance function describes the
correlation between the function values $f(\mathbf{x})$ and
$f(\mathbf{x}^{\prime})$ at any two points $\mathbf{x}$ and
$\mathbf{x}^{\prime}$. Careful design of the covariance function can ensure
that only functions satisfying desired characteristics belong to the
distribution.

Here, we wish to design the covariance function to ensure that only strain
fields satisfying equilibrium are contained in the GP prior. Following the
formulation in (Jidling et al., 2018) a GP model for the Airys stress
functions is defined;
$\varphi(\mathbf{x})\sim\mathcal{GP}\left(0,k_{\varphi}(\mathbf{x},\mathbf{x}^{\prime})\right)$.
Under the assumption that the sample is plane stress, isotropic, and
contiguous, the Airy’s stress functions can be related to strain through the
mapping

$\bar{\boldsymbol{\epsilon}}(\mathbf{x})=\mathcal{V}^{\mathbf{x}}\varphi(\mathbf{x}),\qquad\mathcal{V}^{\mathbf{x}}=\begin{bmatrix}\frac{\partial^{2}}{\partial
y^{2}}-\nu\frac{\partial^{2}}{\partial x^{2}}\\\
-(1+\nu)\frac{\partial^{2}}{\partial x\partial y}\\\
\frac{\partial^{2}}{\partial x^{2}}-\nu\frac{\partial^{2}}{\partial
y^{2}}\end{bmatrix},$ (7)

where $\mathcal{V}$ is a linear operator, and the superscript denotes which
variable the operator acts on. As GPs are closed under linear operators
(Papoulis and Pillai, 2002; Hennig and Kiefel, 2013; Wahlström, 2015; Jidling
et al., 2018) a GP model for strain that satisfies equilibrium can now be
defined;

$\bar{\boldsymbol{\epsilon}}(\mathbf{x})\sim\mathcal{GP}\left(0,\mathcal{V}^{\mathbf{x}}k_{\varphi}(\mathbf{x},\mathbf{x}^{\prime})\mathcal{V}^{\mathbf{x}^{\prime}}{}^{\mathsf{T}}\right),$
(8)

where a prior mean function of zero has been chosen.

Additionally, $d_{0}$ function is also modelled by a GP;

$d_{0}(\mathbf{x})\sim\mathcal{GP}(\bar{d}_{0},k_{d_{0}}(\mathbf{x},\mathbf{x}^{\prime})).$
(9)

where the prior mean $\bar{d}_{0}$ is chosen to be close to the expected
theoretical stress-free lattice spacing for the material used or a measured
average in a stress-free sample. The choice of prior mean function does not
mean that we believe the $d_{0}$ and $\bar{\boldsymbol{\epsilon}}$ functions
to be a particular value, but rather that we do not have any information to
suggest otherwise. After the inclusion of measurement information, the mean of
the posterior estimate will be updated.

There exists a number of options for the base covariance functions
$k_{\varphi}(\mathbf{x},\mathbf{x}^{\prime})$ and
$k_{d_{0}}(\mathbf{x},\mathbf{x}^{\prime})$, with both the squared-exponential
and the Matérn covariance functions having been successfully used for strain
estimation (Jidling et al., 2018; Hendriks et al., 2018a). For a more thorough
discourse on available covariance functions the reader is referred to
Rasmussen and Williams (2006).

Having defined suitable GP models for the strain and $d_{0}$ fields we now
wish to estimate these fields from the LRT and traction measurements. However,
the LRT is a non-linear function of these fields and consequently a closed
form solution does not exist. The following presents a method for obtaining
these estimates that approximates the GP by a finite number of basis functions
allowing variational inference to be applied.

###  3.2 Hilbert Space Approximation to the GP Prior

Here, we make use of the approximation method proposed by (Solin and Särkkä,
2014) and demonstrated to be suitable for the problem of strain tomography
(Jidling et al., 2018). This method approximates our covariance function by a
finite sum of $m$ basis functions;

$k(\mathbf{x},\mathbf{x}^{\prime})=\sum_{j=1}^{m}\phi_{i}(\mathbf{x})S(\boldsymbol{\lambda}_{j})\phi_{j}(\mathbf{x}^{\prime}),\qquad$
(10)

where $S$ is the spectral density of the covariance function. For a stationary
covariance function $k=k(\mathbf{r})$, where
$\mathbf{r}=\mathbf{x}-\mathbf{x}^{\prime}$, the spectral density and the
basis functions are given by;

$S(\boldsymbol{\omega})=\int
k(\mathbf{r})e^{-i\boldsymbol{\omega}^{\mathsf{T}}\mathbf{r}}\,\mathrm{d}\mathbf{r},\quad\phi_{j}=\frac{1}{\sqrt{L_{x}L_{y}}}\sin(\lambda_{x,j}(x+L_{x}))\sin(\lambda_{y,j}(y+L_{y})),$
(11)

where $L_{x}$ and $L_{y}$ control the domain size, and
$\boldsymbol{\lambda}=[\lambda_{x},\lambda_{y}]^{\mathsf{T}}$ encodes spatial
frequencies of the basis functions. The basis functions are chosen as a
solution to the Dirichlet boundary conditions on a rectangular domain, which
is a natural choice for the Laplace eigenvalue problem that needs to be solved
to approximate the GP (Jidling et al., 2018). The parameters
$\theta=\\{l_{x},l_{y},\sigma_{f}\\}$ are commonly called ‘hyperparameters’
and can be chosen by optimisation (as discussed in Section 5.2). For our
application the domain size and spatial frequencies are chosen such that the
basis functions spanned a region where their spectral densities, were greater
than a minimum threshold. This helps to ensure that the dominant frequencies
of the response are captured while maintaining numerical stability.

At this stage, the alternative view point of Bayesian linear regression can be
taken. This approach models the unknown function by a set of basis functions
with Gaussian coefficients;

$f(\mathbf{x})=\sum_{j=1}^{m}\phi_{j}(\mathbf{x})w_{j}=\boldsymbol{\phi}(\mathbf{x})\mathbf{w},\qquad
w_{j}\sim\mathcal{N}(\mu_{j},S(\boldsymbol{\lambda}_{j})),$ (12)

where $\boldsymbol{\phi}(\mathbf{x})$ and $\mathbf{w}$ have dimensions $[1,m]$
and $[m,1]$, respectively. This gives the following model for the strain field
$\bar{\boldsymbol{\epsilon}}(\mathbf{x})$ and the stress-free lattice spacing
$d_{0}(\mathbf{x})$;

$\begin{split}&\bar{\boldsymbol{\epsilon}}_{*}(\mathbf{x})=\boldsymbol{\phi}_{\epsilon}\mathbf{w}_{\varphi},\quad\phi_{\epsilon,j}(\mathbf{x})=\mathcal{V}^{\mathbf{x}}\phi_{\varphi,j}(\mathbf{x}),\quad\phi_{\varphi,j}=\frac{1}{\sqrt{L_{\varphi
x}L_{\varphi y}}}\sin(\lambda_{\varphi x,j}(x+L_{\varphi
x}))\sin(\lambda_{\varphi y,j}(y+L_{\varphi y})),\\\
&d_{0*}(\mathbf{x})=\boldsymbol{\phi}_{d_{0}}\mathbf{w}_{d_{0}},\quad\phi_{d_{0},k}(\mathbf{x})=\frac{1}{\sqrt{L_{d_{0}x}L_{d_{0}y}}}\sin(\lambda_{d_{0}x,k}(x+L_{d_{0}x}))\sin(\lambda_{d_{0}y,k}(y+L_{d_{0}y})),\end{split}$
(13)

where the unknown coefficients are independently normally distributed to
approximate our GP model;
$w_{\varphi,j}\sim\mathcal{N}(0,S_{\varphi}(\boldsymbol{\lambda}_{\varphi,j}))$
and
$w_{d_{0},j}\sim\mathcal{N}(\mu_{d_{0},j},S_{d_{0}}(\boldsymbol{\lambda}_{d_{0},j}))$.
Where the means $\mu_{d_{0},j}$ are chosen so that the prior has the constant
value $\bar{d}_{0}$. In this work, basis functions and parameters
corresponding to the $d_{0}$ field will be denoted by the subscript $d_{0}$
and the subscript $k$ will be used as an index. Likewise, basis functions and
parameters corresponding to the Airys stress function will be denoted by the
subscript $d_{0}$ and the subscript $j$ will be used as an index. The expanded
expressions for $\boldsymbol{\phi}_{\epsilon}$ are given in Appendix A.

Using the LRT (3) we can write a model for a predicted measurement as a non-
linear function of the unknown coefficients;

$\begin{split}y_{*}&=\frac{1}{L}\int\limits_{0}^{L}\bar{\mathbf{n}}\left(\sum_{j}\sum_{k}\phi_{\epsilon,j}(\mathbf{p}+\hat{\mathbf{n}}s)w_{\varphi,j}\phi_{d_{0},k}(\mathbf{p}+\hat{n}s)w_{d_{0},k}\right)+\left(\sum_{k}\phi_{d_{0},k}(\mathbf{p}+\hat{n}s)w_{d_{0},k}\right)\,\mathrm{d}s\\\
&=g_{y}(\mathbf{w}_{\varphi},\mathbf{w}_{d_{0}},\boldsymbol{\eta}),\end{split}$
(14)

where we have restricted ourselves to a single measurement to simplify the
notation. These integrals can be analytically evaluated and the equations are
given in Appendix A. Predictions of the boundary tractions $\mathbf{y}_{t}$ at
a boundary location $\mathbf{x}_{b}$ with surface normal $\mathbf{n}_{\perp}$
can be written as a linear function of the unknown coefficients;

$\begin{split}y_{t*}&=\underbrace{\begin{bmatrix}n_{\perp 1}&n_{\perp 2}&0\\\
0&n_{\perp 1}&n_{\perp
2}\end{bmatrix}}_{\mathbf{T}}\boldsymbol{\phi}_{\varphi}(\mathbf{x})\mathbf{w}_{\varphi}=\mathbf{T}(\mathbf{n}_{\perp})\boldsymbol{\phi}_{\varphi}(\mathbf{x}_{b})\mathbf{w}_{\varphi}\\\
&=g_{t}(\mathbf{w},\mathbf{x}_{b},\mathbf{n}_{\perp}).\end{split}$ (15)

The coefficients $\mathbf{w}_{\varphi}$ and $\mathbf{w}_{d_{0}}$ are random
variables; as such the predictions $\bar{\boldsymbol{\epsilon}}_{*}$, $y_{*}$,
and $y_{t*}$ are also random variables. The problem now is to determine the
distribution of the coefficients given a set of LRT and traction measurements.
This problem is now in a form allowing variational inference to be used to
approximate a solution to the non-linear problem.

###  3.3 Variational Inference

Variational inference (Blei et al., 2017; Jordan et al., 1999) provides an
approximation to the posterior distribution by assuming that it has a certain
functional form that contain unknown parameters. These unknown parameters are
found using optimization, where some distance measure is minimized. We will in
this section provide the details enabling the use of variational inference in
solving our problem.

Given $n$ transmission measurements and $n_{t}$ traction measurements, such
that a vector of all measurements is given by
$\mathbf{Y}=[y_{1},\dots,y_{n},y_{t,1},\dots,y_{t,n_{t}}]^{\mathsf{T}}$, the
problem can be written as having prior and likelihood

$\begin{split}p(\mathbf{w})\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}_{p})\
\ \text{and}\ \
p(\mathbf{Y}|\mathbf{w})\sim\mathcal{N}(\mathbf{Y}|g(\mathbf{w}),\boldsymbol{\Sigma}_{n}),\\\
\end{split}$ (16)

where
$\mathbf{w}=\begin{bmatrix}\mathbf{w}_{d_{0}}^{\mathsf{T}}&\mathbf{w}_{\varphi}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$,
$\boldsymbol{\mu}$ is a vector of all the prior means, and
$\boldsymbol{\Sigma}_{p}$ is a matrix with the coefficients prior variance on
the diagonals. Here, $g(\cdot)$ is the combined measurement model that
expresses the measurement vector $\mathbf{Y}$ as a function of the
coefficients. This function is constructed using both (14) and (15). Finally,
$\boldsymbol{\Sigma}_{n}=\text{diag}(\sigma_{n}^{2}I_{n\times
n},\sigma_{t}^{2}I_{n_{t}\times n_{t}})$, where $\sigma_{t}^{2}$ is a small
variance placed on the artificial traction measurements added for numerical
reasons.

The non-linear measurement function $g(\cdot)$ makes the likelihood
intractable as the prior and likelihood are no longer conjugate. Consequently,
the posterior $p(\mathbf{w}|\mathbf{Y})$ is also intractable and so we find an
approximate solution using variational inference (Jordan et al., 1999). The
idea is to approximate the true posterior by the Gaussian distribution
$q(\mathbf{w})\sim\mathcal{N}(\hat{\mathbf{w}},\mathbf{C})$, and find the mean
$\hat{\mathbf{w}}$ and covariance $\mathbf{C}$ for this distribution that
maximise the Free Energy $\mathcal{F}$. The Free Energy places a lower bound
on the log marginal likelihood and hence provides a measure of how well our
posterior fits the data;

$\log p(\mathbf{Y})\geq\mathbb{E}\left[\log
p(\mathbf{Y}|\mathbf{W})\right]-KL\left[q(\mathbf{w})||p(\mathbf{w}|\mathbf{Y})\right]=\mathcal{F}$
(17)

where, in this case, $\mathbb{E}[\cdot]$ is the expected value with respect to
the approximate posterior $q(\mathbf{w})$ and $KL[\cdot]$ is the Kullback
Leibler divergence which provides a measure of difference between the
approximate posterior and the true posterior. These terms can be evaluated as
(Steinberg and Bonilla, 2014);

$\begin{split}\mathbb{E}\left[\log
p(\mathbf{Y}|\mathbf{W})\right]&=\frac{1}{2}\left[N\log
2\pi+\log|\boldsymbol{\Sigma}_{n}|+(\mathbf{Y}-\mathbb{E}\left[g(\mathbf{w})\right])^{\mathsf{T}}\boldsymbol{\Sigma}_{n}^{-1}(\mathbf{Y}-\mathbb{E}\left[g(\mathbf{w})\right])\right],\\\
KL\left[q(\mathbf{w})||p(\mathbf{w}|\mathbf{Y})\right]&=\frac{1}{2}\left[\text{tr}(\boldsymbol{\Sigma}_{p}^{-1}\mathbf{C})+\left(\boldsymbol{\mu}-\hat{\mathbf{w}}\right)^{\mathsf{T}}\boldsymbol{\Sigma}_{p}^{-1}\left(\boldsymbol{\mu}-\hat{\mathbf{w}}\right)-\log|\mathbf{C}|+\log|\boldsymbol{\Sigma}_{p}|-N\right],\end{split}$
(18)

where $N=n+n_{t}$. Here, the expectation of the non-linear function
$\mathbb{E}\left[g(\mathbf{w})\right]$ is intractable (Steinberg and Bonilla,
2014) and so the expected maximum is used
$\hat{\mathbf{Y}}=g(\hat{\mathbf{w}})$;

$\begin{split}\mathcal{F}&\approx-\frac{1}{2}\big{[}N\log
2\pi+\log|\boldsymbol{\Sigma}_{n}|-\log|\mathbf{C}|+\log|\boldsymbol{\Sigma}_{p}|+(\mathbf{Y}-g(\hat{\mathbf{w}}))^{\mathsf{T}}\boldsymbol{\Sigma}_{n}^{-1}(\mathbf{Y}-g(\hat{\mathbf{w}}))+\left(\boldsymbol{\mu}-\hat{\mathbf{w}}\right)^{\mathsf{T}}\boldsymbol{\Sigma}_{p}^{-1}\left(\boldsymbol{\mu}-\hat{\mathbf{w}}\right)\big{]}\end{split}$
(19)

The optimal posterior mean is chosen to maximise the Free Energy. To perform
this optimisation a modified Newton’s method is used where the step direction
is $\mathbf{q}=-\mathbf{H}^{-1}\mathbf{g}$ and we can calculate the gradient,
$\mathbf{g}$, and Hessian, $\mathbf{H}$, of the cost as

$\begin{split}\mathbf{g}&=\mathbf{J}^{\mathsf{T}}\boldsymbol{\Sigma}_{n}^{-1}(\mathbf{Y}-g(\hat{\mathbf{w}}))-\boldsymbol{\Sigma}_{p}^{-1}\hat{\mathbf{w}}\\\
\mathbf{H}&=-\mathbf{J}^{\mathsf{T}}\boldsymbol{\Sigma}_{n}^{-1}\mathbf{J}+\frac{\partial\mathbf{J}^{\mathsf{T}}}{\partial\mathbf{w}}\Sigma_{n}^{-1}(\mathbf{Y}-g(\hat{\mathbf{w}}))-\boldsymbol{\Sigma}_{p}^{-1}\end{split}$
(20)

where
$\mathbf{J}=\begin{bmatrix}\frac{\partial\hat{\mathbf{Y}}}{\partial\hat{\mathbf{w}}}^{\mathsf{T}}&\frac{\partial\hat{\mathbf{Y}_{t}}}{\partial\hat{\mathbf{w}}}^{\mathsf{T}}\end{bmatrix}^{\mathsf{T}}$,
and the derivatives and second derivatives are given in Appendix B. At each
iteration we update the coefficients according to

$\hat{\mathbf{w}}_{k+1}=(1-\alpha)\hat{\mathbf{w}}_{k}+\alpha\mathbf{q}+\alpha\boldsymbol{\mu},$
(21)

A backwards line search is used to ensure that $\mathcal{F}$ is increased in
each iteration. Once the optimal posterior mean is found, the covariance can
be found by setting $\frac{\partial\mathcal{F}}{\partial\mathbf{C}}=0$ and
linearising about $\hat{\mathbf{w}}$ (Steinberg and Bonilla, 2014), giving;

$\mathbf{C}=\left[\boldsymbol{\Sigma}_{p}^{-1}+\mathbf{J}^{\mathsf{T}}\Sigma_{n}^{-1}\mathbf{J}\right]^{-1}.$
(22)

Pseudo-code for an algorithm to find approximate distribution of the
coefficients $q(\mathbf{w})\sim\mathcal{N}(\hat{\mathbf{w}},\mathbf{C})$ is
given in Algorithm 1. Once the coefficients are found, estimates of the strain
and $d_{0}$ fields can be estimated. The approximate poster mean and variance
for the strain field and stress-free lattice spacing can be computed as

$\begin{split}\begin{bmatrix}\hat{d}_{0}(\mathbf{x})\\\
\hat{\bar{\boldsymbol{\epsilon}}}\end{bmatrix}&=\begin{bmatrix}\boldsymbol{\phi}_{\epsilon}(\mathbf{x})&0\\\
0&\boldsymbol{\phi}_{\epsilon}(\mathbf{x})\end{bmatrix}\hat{\mathbf{w}},\\\
\hat{\Sigma}&=\begin{bmatrix}\boldsymbol{\phi}_{\epsilon}(\mathbf{x})&0\\\
0&\boldsymbol{\phi}_{\epsilon}(\mathbf{x})\end{bmatrix}\mathbf{C}\begin{bmatrix}\boldsymbol{\phi}_{\epsilon}(\mathbf{x})&0\\\
0&\boldsymbol{\phi}_{\epsilon}(\mathbf{x})\end{bmatrix}^{\mathsf{T}},\end{split}$
(23)

where $\hat{\Sigma}$ is the joint covariance of the strain and $d_{0}$
estimates. Next, this method is demonstrated on a set of measurements
simulated from a theoretical cantilever beam strain field and an artificial
$d_{0}$ field.

Algorithm 1 Variational inference algorithm for finding the coefficients
$q(\mathbf{w})\sim\mathcal{N}(\hat{\mathbf{w}},\mathbf{C})$. Requires the
hyperparameters $\theta$, the specified number of basis functions
$m_{\varphi}$ and $m_{d_{0}}$, the LRT measurement information
$\\{y_{i},\boldsymbol{\eta}_{i}|\forall i=1,\dots,n\\}$ and the boundary
traction information
$\\{y_{t,i}=0,\mathbf{x}_{b,i},\mathbf{n}_{\perp,i}|\forall
i=1,\dots,n_{t}\\}$.

1:procedure Find Coefficients

2: Compute the basis functions for the LRT measurements using Equation 14

3: Compute the basis functions for the traction measurements using Equation 15

4: Build prior variance $\boldsymbol{\Sigma}_{p}$

5: Initialise the coefficients $\hat{\mathbf{w}}_{1}$

6: set $k=1$

7: while Stopping criteria not met do

8: Compute the gradient $\mathbf{g}$ and Hessian $\mathbf{H}$ linearised about
$\hat{\mathbf{w}}_{k}$ according to Equation 20

9: Calculate $\hat{\mathbf{w}}_{k+1}$ using Equation 21 and a backward line
search

10: $k=k+1$

11: end while

12: Calculate the covariance $\mathbf{C}$ according to Equation 22

13: return $\hat{\mathbf{w}}_{k}$ and $\mathbf{C}$

14:end procedure

##  4 Simulation Results

The method’s ability to jointly reconstruct the strain field and a $d_{0}$
field is demonstrated using simulated measurements. Reconstructions from
measurements simulated through two strain fields is shown; the Saint-Venant
approximate strain field for a cantilver beam, and a Finite Element Analysis
(FEA) strain field from an in-situ loaded C-shape. Additionally, the
consequences of ignoring the $d_{0}$ variation on the reconstruction are shown
by using the linear measurement model and Gaussian process regression method
presented by Jidling et al. (2018) with the addition of traction constraints
as shown in Hendriks et al. (2018a). Matlab code to run both examples can be
found on Github (Hendriks, 2019).

###  4.1 Cantilever Beam Example

The method is first demonstrated for the theoretical Saint-Venant cantilever
beam as studied. Assuming plane stress, the Saint-Venant approximation to the
strain field is (Beer et al., 2010):

$\boldsymbol{\mathcal{E}}(\mathbf{x})=\begin{bmatrix}\frac{P}{EI}(L-x)y\\\
-\frac{(1+\nu)P}{2EI}\left(\left(\frac{h}{2}\right)^{2}-y^{2}\right)\\\
-\frac{\nu P}{EI}(L-x)y\end{bmatrix},$ (24)

where the geometry is defined in Figure 2. A synthetic stress-free lattice
spacing field is defined by

$d_{0}(\mathbf{x})=c_{0}\exp\left(-\frac{1}{2}(x-c_{1})^{2}/c_{2}^{2}-\frac{1}{2}(y-c_{3})^{2}/c_{4}^{2}\right)+c_{5},$
(25)

with the parameters given by
$\\{c_{0},c_{1},c_{2},c_{3},c_{4},c_{5}\\}=\\{0.0168,0,7.5\times
10^{-3},7\times 10^{-3},6\times 10^{-3},4.056\\}$. The maximum variation
$c_{0}$ from a constant base value, $c_{5}$, was chosen to reflect the
possible maximum relative variation due to martensitic phase change in $0.8\%$
carbon steel.

Figure 2: Cantilever beam geometry and coordinate system with
$l=20$\mathrm{m}\mathrm{m}$$, $h=10$\mathrm{m}\mathrm{m}$$,
$t=5$\mathrm{m}\mathrm{m}$$, $E=200$\mathrm{GPa}$$, $P=2$\mathrm{kN}$$,
$\nu=0.28$, and $I=\frac{th^{3}}{12}$.

Measurements of the form (3) were simulated for $30$ angles evenly spaced
between $0^{\circ}$ and $180^{\circ}$, with $100$ measurements per angle,
which is on the conservative side based on past experiments (Gregg et al.,
2018; Hendriks et al., 2017). The simulated measurements were corrupted with
zero-mean noise of standard deviation $\sigma_{n}=c_{5}\times 10^{-4}$ which
is equilvalent to $1\text{\times}{10}^{-4}$ standard deviation in strain
representing the typical experimental noise (Hendriks et al., 2017; Gregg et
al., 2018)

Fifty zero-traction measurements were added along the top and bottom of the
cantilever beam for both the presented method and the linear GP regression
method. Results are shown in Figure 3. These results show that the presented
method successfully reconstructs both the strain field and the $d_{0}$ field
with a relative error222$\text{relative
error}=\frac{\text{mean}\left(|\text{true}-\text{estimated}|\right)}{\text{max}\left(|\text{true}|\right)}$
of $0.0057$. By contrast ignoring the presence of a $d_{0}$ variation and
using a linear GP regression method yields a drastically degraded strain
reconstruction with a relative error of $0.3067$.

Figure 3: Simulation results for the cantilever beam strain field. The
estimated strain field using the presented method is shown as well as the
results of assuming a constant $d_{0}$ and applying standard GP regression. In
the case of the presented method, the estimated $d_{0}$ field is also shown.
Strain values are given in $\mu\text{Strain}$.

###  4.2 In-situ Loaded C-shape Sample

The method is now demonstrated on a more complex strain field given by FEA of
a mild steel C-shape sample with geometry defined in Figure 4. The sample was
subjected to a $7\text{\,}\mathrm{kN}$ compressive load distributed over
$5^{\circ}$ arcs and plane stress was assumed for the analysis. The resulting
FEA strain field is shown in Figure 5. This sample and loading conditions
correspond to the experimental setup used by Hendriks et al. (2017).

Figure 4: Geometry of the C-shape sample and in-situ loading $\mathbf{P}$. The
sample has an outer diameter of $20\text{\,}\mathrm{mm}$ and an inner diameter
of $10\text{\,}\mathrm{mm}$ with a $45^{\circ}$ segment removed. The sample
was defined to have $E=$200\text{\,}\mathrm{GPa}$$ and $\nu=0.28$.

Measurements of the form (3) were simulated through the this strain and a
synthetic, smoothly changing $d_{0}$ field is again defined by Equation (25)
with parameters given by
$\\{c_{0},c_{1},c_{2},c_{3},c_{4},c_{5}\\}=\\{0.01,$5\text{\times}{10}^{-3}$,7.5\times
10^{-3},7\times 10^{-3},6\times 10^{-3},4.056\\}$. A total of 60 strain images
were simulated with angles evenly spaced between $0^{\circ}$ and
$180^{\circ}$, and 180 measurements per image. The simulated measurements were
corrupted with zero-mean noise of standard deviation $\sigma_{n}=c_{5}\times
10^{-4}$ which is equivalent to $1\text{\times}{10}^{-4}$ standard deviation
in strain representing the typical experimental noise. A total of $131$ zero-
traction measurements were added around the boundary of the C-shape excluding
the regions within $10^{\circ}$ of the loading points. Reconstruction from the
LRT and traction measurements was performed using both the presented method
and the linear GP regression method, and the results are shown in Figure 5.

The presented method achieves a mean relative error of $0.023$ and it can be
seen that the reconstruction has achieved the correct shape. Whereas assuming
a constant $d_{0}$ value gives a mean relative error of $0.137$ and the
resulting strain fields show incorrect concentrated peaks in the strain field
and areas of tension and compression that are reversed. Despite this
improvement there is still some observable difference. In particular, the
presented method has a concentrated tensile region on the top left boundary of
the C, and does not capture the very concentrated peaks in magnitude on the
inside of the C. These peak strains on the boundary are the hardest for the
algorithm to reconstruct as they are poorly sampled by the LRT; i.e. they make
up only a very small part of each line integral. Additionally, some of this
remaining error is due to systematic error in the simulation of the
measurements. Which are generated by numerically performing a line integral
with each function evaluation being given by an interpolation of the FEA
results.

Figure 5: Simulation results for the FEA C-shape sample strain field. Strain
values are given as $\mu\text{Strain}$. The estimated strain field using the
presented method is shown as well as the results of assuming a constant
$d_{0}$ and applying standard GP regression. In the case of the presented
method, the estimated $d_{0}$ field is also shown.

##  5 Additional Remarks

###  5.1 Sensitivity to the Traction Measurement Variance

Boundary conditions given by unloaded surfaces are a natural inclusion as they
are an artefact of the physical world. This information is included in the
form of artificial measurements of zero traction, however it was found that a
small variance needed to be placed on these measurements and the rate of
convergence was impacted by the size of this variance. Conceptually, this
variance is analogous to a constraint tolerance for optimisation procedures.
Too large a variance (or not enough traction measurements) and the algorithm
may fail to converge to the correct strain field. This indicates that the
traction measurements are ensuring that the problem is observable, which is
supported by the findings of Hendriks et al. (2018b) where the inclusion of
traction measurements allowed a constant $d_{0}$ value to be found as a
hyperparameter. Conversely, too small a variance and the algorithm is unable
to take optimisation steps of significant size, resulting in a large number of
iterations to converge. Methods for optimally choosing this variance is an
avenue for future research.

Despite this, it was found that the algorithm worked well over a reasonable
range of traction variances. Typically the standard deviation of the traction
measurements could be set two orders of magnitude smaller than the measurement
standard deviation or in the range of $1\times 10^{-5}$ to $1\times 10^{-7}$.

###  5.2 Hyperparameter Optimisation

The hyperparameters $\theta=\\{l_{x},l_{y},\sigma_{f}\\}$ can be found by
performing an optimisation using $\mathcal{F}$ as the objective function.
However, the gradients of $\mathcal{F}$ with respect to the GP hyperparameters
are not trivial and so Steinberg and Bonilla (2014) suggests that gradient
free optimisation approaches could be used. In this work, both Bayesian
optimisation (Mockus, 2012) and the Nelda-Mead method (Nelder and Mead, 1965)
were found to work; with the Nelda-Mead method requiring less computation
time.

###  5.3 Computational cost

The majority of the computational burden comes from building the building the
forward model of the LRT measurements. In particular, all the combinations of
the basis function of the stress-free lattice spacing and the strain field
need to be computed for each measurement. The resulting matrix has
$nm_{d_{0}}m_{\varphi}$ elements. However, as this matrix only needs to be
computed once, it is still feasible to solve even with a large number of basis
functions.

##  6 Conclusion

This paper considers an extension of the strain tomography problem where the
stress-free lattice parameter is a known constant, to the more general case
where it is unknown and varies throughout the sample. A method for the joint
reconstruction of a strain field and a varying stress-free lattice parameter
from a set of neutron transmission strain images has been presented. This
method extends the Gaussian process based approach previously used for strain
tomography to the subsequently non-linear problem, and ensures that the
estimated strain fields satisfy equilibrium and can include knowledge of
boundary conditions. This was achieved by reformulating the problem in terms
of basis functions and unknown coefficients. Variational inference was then
employed to find estimates of the coefficients.

The method was tested on a set of simulated data, and importantly, these
results demonstrate that it is possible to perform this joint reconstruction.
Further, the results obtained by ignoring variations in $d_{0}$ and applying
the linear GP regression method are provided and show that this assumption, if
incorrect, severely degrades the accuracy of reconstruction.

Future work will involve planning an experiment to acquire a data set on which
to further evaluate the methods performance.

## Acknowledgements

This work is supported by the Australian Research Council through the
Discovery Project scheme (ARCDP170102324), as well as the Swedish Foundation
for Strategic Research (SSF) via the project _ASSEMBLE_ (contract number:
RIT15-0012), and by the Swedish Research Council via the projects _Learning
flexible models for nonlinear dynamics_ (contract number: 2017-03807) and
_NewLEADS - New Directions in Learning Dynamical Systems_ (contract number:
621-2016-06079).

##  Appendix A Basis Functions

The basis functions for the strain field were defined in (13) as

$\phi_{\epsilon,j}=\begin{bmatrix}\frac{\partial^{2}}{\partial
y^{2}}-\nu\frac{\partial^{2}}{\partial x^{2}}\\\
-(1+\nu)\frac{\partial^{2}}{\partial x\partial y}\\\
\frac{\partial^{2}}{\partial x^{2}}-\nu\frac{\partial^{2}}{\partial
y^{2}}\end{bmatrix}\phi_{\varphi,j}.$ (26)

As such, the components of $\phi_{\epsilon,j}$ can be built from

$\begin{split}\frac{\partial^{2}}{\partial
x^{2}}\phi_{\varphi,j}&=\frac{-\lambda_{\varphi x,j}^{2}}{\sqrt{L_{\varphi
x}L_{\varphi y}}}\sin(\lambda_{\varphi x,j}(x+L_{\varphi
x}))\sin(\lambda_{\varphi y,j}(y+L_{\varphi y})),\\\
\frac{\partial^{2}}{\partial y^{2}}\phi_{\varphi,j}&=\frac{-\lambda_{\varphi
y,j}^{2}}{\sqrt{L_{\varphi x}L_{\varphi y}}}\sin(\lambda_{\varphi
x,j}(x+L_{\varphi x}))\sin(\lambda_{\varphi y,j}(y+L_{\varphi y})),\\\
\frac{\partial^{2}}{\partial x\partial
y}\phi_{\varphi,j}&=\frac{\lambda_{\varphi x,j}\lambda_{\varphi
y,j}}{\sqrt{L_{\varphi x}L_{\varphi y}}}\cos(\lambda_{\varphi
x,j}(x+L_{\varphi x}))\cos(\lambda_{\varphi y,j}(y+L_{\varphi
y})).\end{split}$ (27)

A predicted LRT measurement was defined by (14) as

$y_{*}=\frac{1}{L}\int\limits_{0}^{L}\bar{\mathbf{n}}\left(\sum_{j}\sum_{k}\phi_{\epsilon,j}(\mathbf{p}+\hat{\mathbf{n}}s)w_{\varphi,j}\phi_{d_{0},k}(\mathbf{p}+\hat{n}s)w_{d_{0},k}\right)+\left(\sum_{k}\phi_{d_{0},k}(\mathbf{p}+\hat{n}s)w_{d_{0},k}\right)\,\mathrm{d}s$
(28)

where for clarity we restrict ourselves to a single measurement. Therefore, we
need the components

$\begin{split}&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\,\mathrm{d}s,\\\
&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{\partial^{2}}{\partial
x^{2}}\phi_{\varphi,j}(\mathbf{p+\hat{\mathbf{n}}}s)\right)\,\mathrm{d}s,\\\
&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{\partial^{2}}{\partial
y^{2}}\phi_{\varphi,j}(\mathbf{p+\hat{\mathbf{n}}}s)\right)\,\mathrm{d}s,\\\
&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{\partial^{2}}{\partial
x\partial
x}\phi_{\varphi,j}(\mathbf{p+\hat{\mathbf{n}}}s)\right)\,\mathrm{d}s,\\\
\end{split}$ (29)

To make the expressions briefer, we introduce the notation

$\begin{split}\alpha_{\varphi x}&=\lambda_{\varphi
x,j}(x_{0}+\hat{n}_{1}s+L_{\varphi x}),\quad\alpha_{\varphi
y}=\lambda_{\varphi y,j}(y_{0}+\hat{n}_{2}s+L_{\varphi y}),\\\
\alpha_{d_{0}x}&=\lambda_{d_{0}x,k}(x_{0}+\hat{n}_{1}s+L_{d_{0}x}),\quad\alpha_{d_{0}y}=\lambda_{d_{0}y,k}(y_{0}+\hat{n}_{2}s+L_{d_{0}y}).\\\
\end{split}$ (30)

Giving

$\begin{split}\zeta_{k}=&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\,\mathrm{d}s\\\
=&\frac{1}{\sqrt{L_{d_{0}x}L_{d_{0}y}}}\left(\frac{\sin(\alpha_{d_{0}x}-\alpha_{d_{0}y})}{2(\hat{n}_{1}\lambda_{d_{0}x,k}-\hat{n}_{2}\lambda_{d_{0}y,k})}-\frac{\sin(\alpha_{d_{0}x}+\alpha_{d_{0}y})}{2(\hat{n}_{1}\lambda_{d_{0}x,k}+\hat{n}_{2}\lambda_{d_{0}y,k})}\right)\Bigg{|}_{s=0}^{s=L},\\\
\psi_{1,kj}=&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{\partial^{2}}{\partial
x^{2}}\phi_{\varphi,j}(\mathbf{p+\hat{\mathbf{n}}}s)\right)\,\mathrm{d}s\\\
=&\frac{-\lambda_{\varphi x,j}^{2}}{\sqrt{L_{\varphi x}L_{\varphi
y}L_{d_{0}x}L_{d_{0}y}}}\left(-\Gamma_{1}-\Gamma_{2}+\Gamma_{3}+\Gamma_{4}-\Gamma_{5}+\Gamma_{6}-\Gamma_{7}+\Gamma_{8}\right)\Bigg{|}_{s=0}^{s=L},\\\
\psi_{2,kj}=&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{\partial^{2}}{\partial
y^{2}}\phi_{\varphi,j}(\mathbf{p+\hat{\mathbf{n}}}s)\right)\,\mathrm{d}s\\\
=&\frac{-\lambda_{\varphi y,j}^{2}}{\sqrt{L_{\varphi x}L_{\varphi
y}L_{d_{0}x}L_{d_{0}y}}}\left(-\Gamma_{1}-\Gamma_{2}+\Gamma_{3}+\Gamma_{4}-\Gamma_{5}+\Gamma_{6}-\Gamma_{7}+\Gamma_{8}\right)\Bigg{|}_{s=0}^{s=L},\\\
\psi_{3,kj}=&\int\limits_{0}^{L}\phi_{d_{0},k}(\mathbf{p}+\hat{\mathbf{n}}s)\left(\frac{\partial^{2}}{\partial
x\partial
y}\phi_{\varphi,j}(\mathbf{p+\hat{\mathbf{n}}}s)\right)\,\mathrm{d}s\\\
=&\frac{\lambda_{\varphi x,j}\lambda_{\varphi x,j}}{\sqrt{L_{\varphi
x}L_{\varphi
y}L_{d_{0}x}L_{d_{0}y}}}\left(-\Gamma_{1}+\Gamma_{2}-\Gamma_{3}+\Gamma_{4}-\Gamma_{5}+\Gamma_{6}+\Gamma_{7}-\Gamma_{8}\right)\Bigg{|}_{s=0}^{s=L},\\\
\end{split}$ (31a) where
$\begin{split}\Gamma_{1}&=\frac{\sin(\alpha_{d_{0}x}-\alpha_{\varphi
x}+\alpha_{d_{0}y}+\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}-\hat{n}_{1}\lambda_{\varphi
x,j}+\hat{n}_{2}\lambda_{d_{0}y,k}+\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{2}&=\frac{\sin(\alpha_{d_{0}x}-\alpha_{\varphi
x}-\alpha_{d_{0}y}-\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}-\hat{n}_{1}\lambda_{\varphi
x,j}-\hat{n}_{2}\lambda_{d_{0}y,k}-\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{3}&=\frac{\sin(\alpha_{d_{0}x}-\alpha_{\varphi
x}+\alpha_{d_{0}y}-\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}-\hat{n}_{1}\lambda_{\varphi
x,j}+\hat{n}_{2}\lambda_{d_{0}y,k}-\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{4}&=\frac{\sin(\alpha_{d_{0}x}+\alpha_{\varphi
x}-\alpha_{d_{0}y}-\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}+\hat{n}_{1}\lambda_{\varphi
x,j}-\hat{n}_{2}\lambda_{d_{0}y,k}-\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{5}&=\frac{\sin(\alpha_{d_{0}x}+\alpha_{\varphi
x}+\alpha_{d_{0}y}-\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}+\hat{n}_{1}\lambda_{\varphi
x,j}+\hat{n}_{2}\lambda_{d_{0}y,k}-\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{6}&=\frac{\sin(\alpha_{d_{0}x}-\alpha_{\varphi
x}-\alpha_{d_{0}y}+\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}-\hat{n}_{1}\lambda_{\varphi
x,j}-\hat{n}_{2}\lambda_{d_{0}y,k}+\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{7}&=\frac{\sin(\alpha_{d_{0}x}+\alpha_{\varphi
x}-\alpha_{d_{0}y}+\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}+\hat{n}_{1}\lambda_{\varphi
x,j}-\hat{n}_{2}\lambda_{d_{0}y,k}+\hat{n}_{2}\lambda_{\varphi y,j})},\\\
\Gamma_{8}&=\frac{\sin(\alpha_{d_{0}x}+\alpha_{\varphi
x}+\alpha_{d_{0}y}+\alpha_{\varphi
y})}{8(\hat{n}_{1}\lambda_{d_{0}x,k}+\hat{n}_{1}\lambda_{\varphi
x,j}+\hat{n}_{2}\lambda_{d_{0}y,k}+\hat{n}_{2}\lambda_{\varphi
y,j})}.\end{split}$ (31b) Returning to the measurement model in (14), we can
now write
$y_{*}=\frac{1}{L}\left(\sum_{k}\sum_{j}\bar{\mathbf{n}}w_{d_{0},k}w_{\varphi
j}\begin{bmatrix}\psi_{2,kj}-\nu\psi_{1,kj}\\\ -(1+\nu)\psi_{3,kj}\\\
\psi_{1,kj}-\nu\psi_{2,kj}\end{bmatrix}_{s=0}^{s=L}+\sum_{k}w_{d_{0},k}\zeta_{k}\big{|}_{s=0}^{s=L}\right).$
(31c)

##  Appendix B Measurement Model Derivatives

Here we give the derivates of the measurement model $g(\cdot)$ about the
current $\hat{\mathbf{w}}$. The measurement model is a concatenation of
equations (14) and (15), and so we require the derivatives of the predicted
LRT measurement, $\hat{y}$, and the predicted traction $\hat{y}_{t}$. For
clarity we restrict the following to a single $\hat{y}$ and $\hat{y}_{t}$. The
first derivatives are given by

$\begin{split}&\frac{\partial\hat{y}}{\partial
w_{d_{0},k}}=\frac{1}{L}\left(\sum_{j}\bar{\mathbf{n}}w_{\varphi*,j}\begin{bmatrix}\psi_{2,kj}-\nu\psi_{1,kj}\\\
-(1+\nu)\psi_{3,kj}\\\
\psi_{1,kj}-\nu\psi_{2,kj}\end{bmatrix}_{s=0}^{s=L}+\zeta_{k}\big{|}_{s=0}^{s=L}\right),\\\
&\frac{\partial\hat{y}}{\partial
w_{\varphi,k}}=\frac{1}{L}\left(\sum_{k}\bar{\mathbf{n}}w_{d_{0}*,k}\begin{bmatrix}\psi_{2,kj}-\nu\psi_{1,kj}\\\
-(1+\nu)\psi_{3,kj}\\\
\psi_{1,kj}-\nu\psi_{2,kj}\end{bmatrix}_{s=0}^{s=L}\right),\\\
&\frac{\partial\hat{y}_{t}}{\partial\mathbf{w}_{d_{0}}}=0,\\\
&\frac{\partial\hat{y}_{t}}{\partial\mathbf{w}_{\varphi}}=\mathbf{T}\boldsymbol{\phi}_{\varphi}(\mathbf{x}).\end{split}$
(32)

The second derivates are

$\begin{split}\frac{\partial^{2}\hat{\mathbf{y}}_{t}}{\partial\mathbf{w}^{2}}&=0,\\\
\frac{\partial^{2}\hat{y}}{\partial\mathbf{w}^{2}}&=\begin{bmatrix}0&\frac{\partial^{2}\hat{y}}{\partial\mathbf{w}_{d_{0}}\partial\mathbf{w}_{\varphi}}\\\
\left(\frac{\partial^{2}\hat{y}}{\partial\mathbf{w}_{d_{0}}\partial\mathbf{w}_{\varphi}}\right)^{\mathsf{T}}&0\end{bmatrix},\end{split}$
(33)

where

$\left[\frac{\partial^{2}\hat{y}}{\partial\mathbf{w}_{d_{0}}\partial\mathbf{w}_{\varphi}}\right]_{kj}=\bar{\mathbf{n}}\begin{bmatrix}\psi_{2,kj}-\nu\psi_{1,kj}\\\
-(1+\nu)\psi_{3,kj}\\\ \psi_{1,kj}-\nu\psi_{2,kj}\end{bmatrix}_{s=0}^{s=L}.$
(34)

Explicit formulation of the second derivatives allows the cost functions
curvature to be taken into account in the optimisation procedure, greatly
improving the rate of convergence.

## References

  * Abbey et al. (2009) B Abbey, SY Zhang, W Vorster, and AM Korsunsky.  Feasibility study of neutron strain tomography.  _Procedia Engineering_ , 1(1):185–188, 2009.
  * Abbey et al. (2012) B Abbey, SY Zhang, W Vorster, and AM Korsunsky.  Reconstruction of axisymmetric strain distributions via neutron strain tomography.  _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ , 270:28–35, 2012.
  * Beer et al. (2010) FP Beer, ER Johnston Jr, JT Dewolf, and DF Mazurek.  Mechanics of materials, sixth edit edition, 2010.
  * Bishop et al. (1995) CM Bishop et al.  _Neural networks for pattern recognition_.  Oxford university press, 1995.
  * Blei et al. (2017) D Blei, A Kucukelbir, and J McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American Statistical Association_ , 112(518):859–877, 2017.
  * Choi et al. (2007) WB Choi, L Li, V Luzin, R Neiser, T Gnäupel-Herold, HJ Prask, S Sampath, and A Gouldstone.  Integrated characterization of cold sprayed aluminum coatings.  _Acta Materialia_ , 55(3):857–866, 2007.
  * Gregg et al. (2017) AWT Gregg, JN Hendriks, CM Wensrich, and MH Meylan.  Tomographic reconstruction of residual strain in axisymmetric systems from bragg-edge neutron imaging.  _Mechanics Research Communications_ , 85:96–103, 2017.
  * Gregg et al. (2018) AWT Gregg, JN Hendriks, CM Wensrich, Adrian Wills, AS Tremsin, Vladimir Luzin, Takenao Shinohara, Oliver Kirstein, MH Meylan, and EH Kisi.  Tomographic reconstruction of two-dimensional residual strain fields from bragg-edge neutron imaging.  _Physical Review Applied_ , 10(6):064034, 2018\.
  * Hendriks (2019) J Hendriks.  Joint tomography of strain and the strain-free lattice spacing from neutron transmission strain images.  https://github.com/jnh277/Joint_strain_d0_tomography.git, 2019.  URL https://github.com/jnh277/Joint_strain_d0_tomography.git.
  * Hendriks et al. (2017) JN Hendriks, AWT Gregg, CM Wensrich, AS Tremsin, T Shinohara, M Meylan, EH Kisi, V Luzin, and O Kirsten.  Bragg-edge elastic strain tomography for in situ systems from energy-resolved neutron transmission imaging.  _Physical Review Materials_ , 1(5):053802, 2017\.
  * Hendriks et al. (2018a) JN Hendriks, AWT Gregg, CM Wensrich, and A Wills.  Implementation of traction constraints in bragg-edge neutron transmission strain tomography.  _arXiv preprint arXiv:1805.09760_ , 2018a.
  * Hendriks et al. (2018b) JN Hendriks, CM Wensrich, A Wills, V Luzin, and AW T Gregg.  Robust inference of two-dimensional strain fields from diffraction-based.  _arXiv preprint arXiv:1808.06282_ , 2018b.
  * Hennig and Kiefel (2013) P Hennig and M Kiefel.  Quasi-newton method: A new direction.  _Journal of Machine Learning Research_ , 14(Mar):843–865, 2013.
  * Jidling et al. (2017) C Jidling, N Wahlström, A Wills, and TB Schön.  Linearly constrained Gaussian processes.  In _Advances in Neural Information Processing Systems (NIPS)_ , Long Beach, CA, USA, December 2017.
  * Jidling et al. (2018) C Jidling, J Hendriks, N Wahlström, A Gregg, TB Schön, C Wensrich, and A Wills.  Probabilistic modelling and reconstruction of strain.  _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ , 436:141–155, 2018.
  * Jordan et al. (1999) MI Jordan, Z Ghahramani, TS Jaakkola, and LK Saul.  An introduction to variational methods for graphical models.  _Machine learning_ , 37(2):183–233, 1999.
  * Kirkwood et al. (2015) HJ Kirkwood, SY Zhang, AS Tremsin, AM Korsunsky, N Baimpas, and B Abbey.  Neutron strain tomography using the radon transform.  _Materials Today: Proceedings_ , 2:S414–S423, 2015.
  * Lionheart and Withers (2015) WRB Lionheart and PJ Withers.  Diffraction tomography of strain.  _Inverse Problems_ , 31(4):045005, 2015.
  * Luzin et al. (2011) V Luzin, K Spencer, and M-X Zhang.  Residual stress and thermo-mechanical properties of cold spray metal coatings.  _Acta Materialia_ , 59(3):1259–1270, 2011.
  * Mockus (2012) J Mockus.  _Bayesian approach to global optimization: theory and applications_ , volume 37.  Springer Science & Business Media, 2012.
  * Murray and Adams (2010) I Murray and RP Adams.  Slice sampling covariance hyperparameters of latent gaussian models.  In _Advances in neural information processing systems_ , pages 1732–1740, 2010.
  * Nelder and Mead (1965) JA Nelder and R Mead.  A simplex method for function minimization.  _The computer journal_ , 7(4):308–313, 1965.
  * Papoulis and Pillai (2002) A Papoulis and S Pillai.  _Probability, random variables, and stochastic processes_.  Tata McGraw-Hill Education, 2002.
  * Paradowska et al. (2005) A Paradowska, JWH Price, R Ibrahim, and T Finlayson.  A neutron diffraction study of residual stress due to welding.  _Journal of materials processing technology_ , 164:1099–1105, 2005.
  * Rasmussen and Williams (2006) CE Rasmussen and CKI Williams.  _Gaussian processes for machine learning_ , volume 1.  MIT press Cambridge, 2006.
  * Sadd (2009) MH Sadd.  _Elasticity: theory, applications, and numerics_.  Academic Press, 2009.
  * Santisteban et al. (2002a) JR Santisteban, L Edwards, ME Fitzpatrick, A Steuwer, PJ Withers, MR Daymond, MW Johnson, N Rhodes, and EM Schooneveld.  Strain imaging by bragg edge neutron transmission.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 481(1):765–768, 2002a.
  * Santisteban et al. (2002b) JR Santisteban, L Edwards, ME Fizpatrick, A Steuwer, and PJ Withers.  Engineering applications of Bragg-edge neutron transmission.  _Applied Physics A_ , 74(1):s1433–s1436, 2002b.
  * Solin and Särkkä (2014) A Solin and S Särkkä.  Hilbert space methods for reduced-rank gaussian process regression.  _arXiv preprint arXiv:1401.5508_ , 2014.
  * Steinberg and Bonilla (2014) DM Steinberg and EV Bonilla.  Extended and unscented gaussian processes.  In _Advances in Neural Information Processing Systems_ , pages 1251–1259, 2014.
  * Tremsin et al. (2011) AS Tremsin, JB McPhate, W Kockelmann, JV Vallerga, OHW Siegmund, and WB Feller.  High resolution Bragg edge transmission spectroscopy at pulsed neutron sources: Proof of principle experiments with a neutron counting MCP detector.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 633:S235–S238, 2011.
  * Tremsin et al. (2012) AS Tremsin, JB McPhate, A Steuwer, W Kockelmann, A M Paradowska, JF Kelleher, JV Vallerga, OHW Siegmund, and WB Feller.  High-resolution strain mapping through time-of-flight neutron transmission diffraction with a microchannel plate neutron counting detector.  _Strain_ , 48(4):296–305, 2012.
  * Wahlström (2015) N Wahlström.  _Modeling of Magnetic Fields and Extended Objects for Localization Applications_.  PhD thesis, Linköping University Electronic Press, 2015.


# Adapting Pretrained Language Models For Solving Tabular Prediction Problems
In The Electronic Health Record

Christopher McMaster  Department of Clinical Pharmacology & Therapeutics,
Austin Health, Melbourne, Australia  Department of Rheumatology, Austin
Health, Melbourne, Australia  The Centre for Digital Transformation of Health,
University of Melbourne, Melbourne, Australia  School of Computing and
Information Systems, University of Melbourne, Melbourne, Australia  David FL
Liew  Department of Clinical Pharmacology & Therapeutics, Austin Health,
Melbourne, Australia  Department of Rheumatology, Austin Health, Melbourne,
Australia  Department of Medicine, University of Melbourne, Melbourne,
Australia  Douglas EV Pires  The Centre for Digital Transformation of Health,
University of Melbourne, Melbourne, Australia  School of Computing and
Information Systems, University of Melbourne, Melbourne, Australia

###### Abstract

We propose an approach for adapting the DeBERTa model for electronic health
record (EHR) tasks using domain adaptation. We pretrain a small DeBERTa model
on a dataset consisting of MIMIC-III discharge summaries, clinical notes,
radiology reports, and PubMed abstracts. We compare this model’s performance
with a DeBERTa model pre-trained on clinical texts from our institutional EHR
(MeDeBERTa) and an XGBoost model. We evaluate performance on three benchmark
tasks for emergency department outcomes using the MIMIC-IV-ED dataset. We
preprocess the data to convert it into text format and generate four versions
of the original datasets to compare data processing and data inclusion. The
results show that our proposed approach outperforms the alternative models on
two of three tasks (p<0.001) and matches performance on the third task, with
the use of descriptive columns improving performance over the original column
names.

_K_ eywords Tabular data $\cdot$ Language models $\cdot$ Transformers $\cdot$
Electronic health records $\cdot$ Multimodal data $\cdot$ Clinical informatics

##  1 Introduction

Tabular data is ubiquitous in real world data science problems. Whilst a large
amount of data is naturally generated and collected in tables, such as
internet traffic metadata, numerical data from scientific experiments and
nursing observations in the intensive care unit, the historically poor
performance of many machine learning algorithms on unstructured data has
motivated the development of data wrangling methods to convert even non-
tabular data into tabular format for optimal algorithm performance. But data
wrangling and cleaning to convert unstructrued data into tabular data can
result in the loss of information, particularly where the data is entered as
free text. In healthcare applications, this often involves encoding drugs,
pathology results and diagnoses as factors in data tables, even when these
data are entered in free text fields. For example, when developing a
predictive model using medication prescription data, the process of data
wrangling might involve encoding the medications as factors, or categorising
them into groups, such as antihypertensives, anticoagulants, and
antiplatelets. Some details might be retained, such as dose and/or frequency
of administration, but often the full rich semantic information present in the
free text is lost, such as the route of administration, the indication,
tapering or titration schedules, and a detailed history of previous
medications for each indication.

Pretrained language models have demonstrated excellent performance on a
variety of tasks, including text classification, question answering, and
natural language inference[1, 2, 3]. These models have been trained on large
amounts of text data, and have been shown to be able to learn rich
representations of language, including the semantic meaning of words and
phrases. Recent work has shown that these models can be used to solve tabular
prediction problems by converting tabular data into strings of tokens [4, 5,
6]. However, there are known limitations with how pretrained language models
are used for numeric tokens. Wallace et al. [7] demonstrated that language
models trained to decode token embeddings of numbers (i.e. predict the
floating point number from its embedding) can only accurately recapture within
the range of numbers they were trained on. This could introduce significant
error when applied to numeric data outside of the training range.

In the setting of electronic health record (EHR) prediction problems where the
data is a 20:80 mix of tabular and unstructured (including free text)[8], we
investigate whether language models can be used to solve tabular prediction
problems without the need for extensive data wrangling. Given data privacy
concerns limiting cloud service use and limited local compute, we test whether
competitive performance can be achieved with small language models. Because
pure tabular datasets of EHR data often involve the transformation of free
text into categorical factors, we also investigate whether retaining the
original free text data improves performance. We also investigate strategies
for handling numeric data in the tokeniation and embedding process. Finally,
given the abundance of free text notes in the EHR, we evaluate performance
with and without further pretraining on free text EHR data.

##  2 Related work

Several authors have proposed using language models for tabular data. Liu et
al. [6] proposed a method for using BERT [1] to solve tabular prediction
problems, by converting tabular data into strings of tokens, performing masked
language modelling (MLM) on the resulting strings and then fine-tuning on the
target task. The method was shown to be competitive with state-of-the-art
tree-based methods on a variety of pure tabular prediction problems. LIFT [4]
and TabLLM [5] use large pretrained large language models to tackle tabular
data problems, with a focus on few-shot performance. These methods have
demonstrated not only competitive overall performance with best-in-class tree-
based methods, but also excellent few-shot performance. LIFT uses GPT-J[9] (6
billion parameters) and GPT-3[10] (off-site finetuning via OpenAI API), while
TabLLM uses T0pp[11] (11 billion parameters). The size of GPT-J and T0pp, and
the inavailability of GPT-3 weights, limit their usability in settings with
low compute, particularly where data privacy concerns present a barrier to
train and deploy models on a cloud service.

Levin et al. [12] benchmarked two different tabular transformer
architectures[13, 14] on multiple EHR tasks from the MetaMIMIC repository[15].
The authors explored a number of pretraining regimes and found that supervised
pretraining on multiple diagnostic tasks improves downstream performance on a
different diagnostic task. Whilst these models do use transformer
architectures, they are not language models and therefore cannot be used for
tasks with mixed tabular and free text data. There has been a growing interest
in using deep learning methods for working with tabular data in the context of
clinical research. Previous studies have shown that convolutional neural
networks (CNNs) can be effective for automatically extracting features from
tabular data, and have demonstrated improved performance on a variety of
tasks, such as predicting patient outcomes and identifying disease
subtypes[16, 17].

##  3 Methodology

###  3.1 Models

Our proposed approach is based on the DeBERTa model[3, 18], a variant of
BERT[1] that has been shown to perform well on tasks related to natural
language processing in the biomedical domain[19]. DeBERTa is a deep learning
model that uses a transformer architecture[20] to learn contextual
representations of text, allowing it to capture the rich semantic information
present in clinical data. As a baseline model, we used a pretrained DeBERTa
small model (version 3)111https://huggingface.co/microsoft/deberta-v3-small
(DeBERTa Small), a 44 million parameter model. We then pretrained this model
for domain adaptation on the masked language modelling task, using a dataset
consisting of 700k MIMIC-III discharge summaries, clinical notes and radiology
reports[21, 22], plus 1 million abstracts from the PubMed dataset[23]
published by the National Library of Medicine (Deberta Small [adapted]). This
allowed the model to learn a wide range of clinical and biomedical concepts
with the goal of improving performance on downstream EHR tasks. Additionally,
we compared performance with a DeBERTa base model (verison 1) that we have
previously pretrained on 1.3 million clinical texts from our institutional EHR
at Austin Health[19] (MeDeBERTa). For comparison, we trained XGBoost[24]
models as strong baselines, given the consistently good performance of XGBoost
on tabular data tasks. We also include the results from the best performing
model for each task in the original paper.

###  3.2 Data Processing

To use DeBERTa for working with tabular data, we first preprocess the data to
convert it into text format. This involves refining the column headings to
make them more descriptive and then concatenating these with the values. Where
several columns naturally belong to a list (e.g. one-hot-encoded comorbidity
columns), these can be combined into a comma-separated list of features. For
example, the following comma-separated values (CSV) file,

⬇

1cci_DM2,cci_Cancer1,age

21,1,50

30,1,73

would be converted into the following text representations:

⬇

1Age in years: 50, Comorbidities: Diabetes Mellitus with Complications, Cancer
Without Metastases

2Age in years: 73, Comorbidities: Cancer Without Metastases

To evaluate the performance of our proposed approach, we used the Medical
Information Mart for Intensive Care IV Emergency Department (MIMIC-IV-ED)
dataset[25, 26, 27] to derive three benchmark tasks for emergency department
outcomes, as described by Xie et al. [28]. These tasks are as follows:

  * •

Task 1: The hospitalization outcome predicts whether a patient will be
admitted to the hospital following an ED visit.

  * •

Task 2: The critical outcome predicts whether a patient will require urgent
care, such as transfer to an ICU within 12 hours.

  * •

Task 3: The ED reattendance outcome predicts whether a patient will return to
the ED within 72 hours after their previous discharge.

To generate the datasets for these tasks, we reuse the code from Xie et al.
[28]222https://github.com/nliulab/mimic4ed-benchmark. In their code, the
authors use snake case for column names, using abbreviations and shorthand to
encode longer descriptions (e.g. "n_ed_90d" for the number of emergency
department presentations in the previous 90 days). Additionally, the authors
used a rule-based method to extract 10 common chief complaints from the free
text description (abdominal pain, dizziness, etc.), converting this into 10
one-hot-encoded columns representing the presence or absence of each of these
complaints.

We generated 4 versions of the original datasets to compare data processing
and data inclusion:

  * •

Values only: Text strings of comma-separated values (e.g. "99.4, 76, 16")

  * •

Original columns: Text strings of original column names and values separated
by a colon, with name/value pairs joined in a semi-colon separated list (e.g.
"triage_temperature: 99.4; triage_heartrate: 76; triage_resprate: 16")

  * •

Descriptive columns: Text strings of more descriptive column names and values
separated by a colon, with name/value pairs joined in a semi-colon separated
list (e.g. "triage temperature: 99.4; triage heart rate: 76; triage
respiratory rate: 16")

  * •

Descriptive columns, text problems: The same as above, with the free text
chief complaint added as an extra data point.

The datasets were split into train and test sets, as per the original paper.

###  3.3 Task Finetuning

Models were finetuned on each of the three task separately. Finetuning was
performed by splitting the predefined training sets into train (98%) and
validation (2%) sets. Each model was trained for 20 epochs, with evaluation
performed and model checkpoints saved after every 500 steps. The checkpoint
with the best performance on the validation set was then loaded and evaluation
performed on the predefined test set defined by Xie et al. [28].

The XGBoost models were trained using the same training and test sets.
Hyperparameters were tuned on the training set using 5-fold cross validation
and Bayesian optimisation.

###  3.4 Model Comparison

Models were compared using the area under the receiver operating
characteristic curve (AUC) metric, with bootstrapped confidence intervals
calculated using the method used in the original paper. To compare the
performance of the different data processing methods, we used the two-sample
t-test.

###  3.5 Explainability

To investigate the contributions of each token at the individial level and
across samples, we calculated layer integrated gradients with respect to the
embedding layer[29, 30] using the Captum library[31]. Here, the integrated
gradient for the $i$th token is calculated by taking the path integral of the
layer gradients between two embeddings $x$ and $x^{\prime}$, where
$x^{\prime}$ is the embedding of the pad token and $x$ is the embedding of the
token of interest:

$IG_{i}=(x_{i}-x^{\prime}_{i})\times(\int_{0}^{1}\frac{\partial\mathbf{F}(x^{\prime}+\alpha(x-x^{\prime}))}{\partial
x_{i}}d\alpha)$ (1)

In practice, the integral is approximated by taking the average of the
gradients at $n$ points along the path between $x$ and $x^{\prime}$. Averaging
these gradients across the embedding dimensions gives a token-specific
attribution score.

To investigate variable-level attribuition, we classify each token by the
variable to which it belongs (e.g. in "Co morbidities : Co agul opathy ,
Dementia" the following tokens are comorbidity tokens: "Co", "agul", "opathy",
"Dementia"). We then calculate the maximum value of the attribution scores for
each token to get a variable-level attribution score. We do this for a random
sample of 1,000 samples from the test set, and then average the attribution
scores across these samples.

##  4 Results

###  4.1 Performance

In all three tasks, we were able to improve on the performance of the best
performing model from the original publication (see Table 1). The addition of
column names improved performance for every task, although improvements were
small for tasks 1 and 2. More descriptive column names had a smaller impact on
performance, however these column names produced shorter tokenized sequences
(owing to the tokenizer vocabulary). To minimize sequence length, the
descriptive column names were retained when adding the free text chief
complaint as an additional data point. This addition improved performance
across all three tasks, although the improvement was largest for task 1.

| Model | DeBERTa Small | DeBERTa Small [adapted] | MeDeBERTA | XGBoost | Original Paper[28]
---|---|---|---|---|---|---
Task | Dataset |  |  |  |  |
1 | Descriptive columns, text problem | 0.864 (0.862, 0.866) | 0.865 (0.862, 0.867) | 0.866 (0.864, 0.868) | - | -
| Descriptive columns | 0.827 (0.825, 0.83) | 0.827 (0.824, 0.829) | - | - | -
| Original columns | 0.828 (0.826, 0.831) | 0.827 (0.825, 0.83) | - | - | -
| Values only | 0.828 (0.825, 0.83) | 0.827 (0.825, 0.83) | - | 0.828 (0.828, 0.828) | 0.822 (0.821, 0.825)
2 | Descriptive columns, text problem | 0.909 (0.904, 0.913) | 0.91 (0.906, 0.914) | 0.911 (0.907, 0.915) | - | -
| Descriptive columns | 0.889 (0.884, 0.893) | 0.889 (0.884, 0.894) | - | - | -
| Original columns | 0.889 (0.884, 0.893) | 0.888 (0.883, 0.892) | - | - | -
| Values only | 0.886 (0.882, 0.891) | 0.888 (0.883, 0.892) | - | 0.891 (0.891, 0.891) | 0.883 (0.879, 0.888)
3 | Descriptive columns, text problem | 0.702 (0.690, 0.714) | 0.705 (0.693, 0.716) | 0.699 (0.687, 0.711) | - | -
| Descriptive columns | 0.689 (0.676, 0.702) | 0.690 (0.677, 0.704) | - | - | -
| Original columns | 0.694 (0.682, 0.707) | 0.694 (0.681, 0.707) | - | - | -
| Values only | 0.653 (0.639, 0.664) | 0.665 (0.651, 0.68) | - | 0.703 (0.703, 0.703) | 0.700 (0.691, 0.713)

Table 1: Area Under the Receiver Operating Characteristic (ROC-AUC) for Each
Task

Pretraining DeBERTa on a combination of MIMIC-III and PubMed resulted in small
improvements, seen most prominently in task 3. Overall, our previously
published MeDeBERTa model using the free text chief complaint achieved the
highest ROC-AUC on tasks 1 and 2 Xie et al. [28] by 0.043 (P<0.001) and 0.031
(p<0.001) respectively. On task 3, although the mean ROC-AUC was 0.007 higher
for the DeBERTa Small [adapted] model, this difference was not statistically
significant (p=0.538).

Figure 1: Receiver Operating Characteristic Curves for Each Task

###  4.2 Explainability

Tokenwise attribution scores were calculated using the embedding gradients to
interrogate the contributions of different variables and column headings.
Table 2 demonstrates a visualization of the attribution scores for 2 different
example patients, with and without the free text chief complaint field. Whilst
the model (DeBERTa Small [adapted] in this case) correctly predicted
hospitalization outcome for both patients, in the absence of a chief complaint
that fit any of the 10 predefined categories, the model attribution scores for
the free text chief complaints were high and contributed to higher model
confidence.

Table 2: Word importance for two example patients with model predictions for
task 1 (hospitalization), with and without the free text chief complaint.
Figure 2: Mean attribution scores across 5,000 randomly sampled patients in
the test set for models 1, 2 and 3.

Grouping tokens into variables headings, 5,000 examples from the test set were
sampled and the average maximum attribution score was calculated for each
variables (Figure 2). The free text chief complaints was, on average, the
variable with the largest attribution scores.

##  5 Discussion

In this study, we evaluated the effectiveness of our proposed approach, which
is based on the DeBERTa model, in predicting emergency department outcomes. We
found that the DeBERTa small [adapted] model outperformed the baseline DeBERTa
small and XGBoost models on all three benchmark tasks and showed similar
performance to a larger DeBERTa model pretrained on a different corpus of
clinical notes.

Our results showed that the inclusion of free text data and the use of the
d̈escriptive columnsv̈ersion of the dataset resulted in the best model
performance. This highlights the importance of data processing and the ability
to mix free text and tabular data in EHR datasets. The ability to mix free
text and tabular data allows for a more comprehensive representation of
contributing risk factors, which can lead to improved model performance and
decision-making.

We demonstrated that small DeBERTa models can achieve competitive performance
when fine-tuned on EHR datasets, which is important for compute-constrained
settings such as hospitals. For comparison, a 6 billion parameter model like
GPT-J[9] will require approximately 22GB of RAM using single-precision
floating point numbers for the weights alone, whereas the DeBERTa small model
used in this study requires approximately 0.17GB of RAM.

We used attribution scores to determine the importance of different input
features in our machine learning models 2. The attribution scores align with
our understanding of the factors leading to the three outcomes we predicted.
For example, the free text chief complaint was the most important variable in
predicting the likelihood of hospitalization, which is consistent with the
fact that the decision to admit a patient depends largely on the type of
presentation. The attribution scores for the models predicting the likelihood
of needing urgent care showed a reliance on patient vital signs, particularly
respiratory rate, which is consistent with the fact that urgent care is often
used to support patients with physiological decompensation. Finally,
predicting the likelihood of representation within 72 hours was strongly
attributed to the number of prior ED, hospital and ICU stays, particularly in
the prior 30 days, which is consistent with the fact that patients with a
history of frequent visits are more likely to have future visits.

Attribution scores, particularly when determined at the token level, may hold
significant clinical value in shedding light on the factors influencing
patient outcomes and pinpointing potentially modifiable risk factors. The
explainability provided by these scores may help in guiding the
appropriateness and impact of clinical interventions. For instance, within
patient subgroups where EHR notes exhibit high attribution scores related to
medications (including dosing, frequency and indication), it is plausible that
deprescribing interventions could have a disproportionately positive impact.
By identifying patients with high polypharmacy or elevated
sedative/anticholinergic burdens, healthcare professionals may leverage these
attribution scores to optimize medication regimens, reduce adverse drug
interactions, and improve overall patient safety.

Although the general approach of using small language models for mixed
structured and free text EHR data is promising, we have not presented a
comprehensive evaluation on a wide range of tasks and have not directly
compared our approach to much larger models. We are planning to conduct a more
comprehensive evaluation of our approach on a wider range of tasks and compare
performance to larger models.

We have presented an approach for adapting the DeBERTa model for EHR tasks
using domain adaptation. Our proposed approach outperforms other models on
three benchmark tasks for emergency department outcomes using the MIMIC-IV-ED
dataset. We also compared our model’s performance with a DeBERTa model pre-
trained on clinical texts from our institutional EHR and an XGBoost model. Our
results demonstrate the effectiveness of our approach and highlight the
importance of data preprocessing and inclusion. Our study provides a
foundation for future work in adapting transformer-based models for EHR tasks
and highlights the potential of these models for improving prediction models
and therefore outcomes for hospitalised patients.

## References

  * Devlin et al. [2018] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the NAACL-HLT Conference_ , pages 4171–4186, 2018\.
  * Liu et al. [2019] Yiren Liu, Myle Ott, Nikhil Goyal, Jiajun Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_ , 2019.
  * He et al. [2021a] Pengcheng He, Xiaodong Liu, Jianfeng Gao, and Weizhu Chen.  Deberta: Decoding-enhanced bert with disentangled attention.  In _International Conference on Learning Representations_ , 2021a.  URL https://openreview.net/forum?id=XPZIaotutsD.
  * Dinh et al. [2022] Tuan Dinh, Yuchen Zeng, Ruisu Zhang, Ziqian Lin, Michael Gira, Shashank Rajput, Jy-Yong Sohn, Dimitris Papailiopoulos, and Kangwook Lee.  LIFT: Language-Interfaced Fine-Tuning for Non-Language machine learning tasks.  _arXiv_ , June 2022.
  * Hegselmann et al. [2022] Stefan Hegselmann, Alejandro Buendia, Hunter Lang, Monica Agrawal, Xiaoyi Jiang, and David Sontag.  TabLLM: Few-shot classification of tabular data with large language models.  _arXiv_ , October 2022.
  * Liu et al. [2022] Guang Liu, Jie Yang, and Ledell Wu.  PTab: Using the pre-trained language model for modeling tabular data.  _arXiv [cs.LG]_ , September 2022.
  * Wallace et al. [2019] Eric Wallace, Yizhong Wang, Sujian Li, Sameer Singh, and Matt Gardner.  Do NLP models know numbers? probing numeracy in embeddings.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 5307–5315, Hong Kong, China, November 2019. Association for Computational Linguistics.  doi: 10.18653/v1/D19-1534.  URL https://aclanthology.org/D19-1534.
  * Kong [2019] Hyoun-Joong Kong.  Managing unstructured big data in healthcare system.  _Healthcare informatics research_ , 25(1):1–2, January 2019.  ISSN 2093-3681.  doi: 10.4258/hir.2019.25.1.1.
  * Wang [2021] Ben Wang.  Mesh-Transformer-JAX: Model-Parallel Implementation of Transformer Language Model with JAX.  https://github.com/kingoflolz/mesh-transformer-jax, May 2021.
  * Brown et al. [2020] Tom B Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei.  Language models are Few-Shot learners.  _arXiv_ , May 2020.
  * Sanh et al. [2021] Victor Sanh, Albert Webson, Colin Raffel, Stephen H Bach, Lintang Sutawika, Zaid Alyafeai, Antoine Chaffin, Arnaud Stiegler, Teven Le Scao, Arun Raja, Manan Dey, M Saiful Bari, Canwen Xu, Urmish Thakker, Shanya Sharma Sharma, Eliza Szczechla, Taewoon Kim, Gunjan Chhablani, Nihal Nayak, Debajyoti Datta, Jonathan Chang, Mike Tian-Jian Jiang, Han Wang, Matteo Manica, Sheng Shen, Zheng Xin Yong, Harshit Pandey, Rachel Bawden, Thomas Wang, Trishala Neeraj, Jos Rozen, Abheesht Sharma, Andrea Santilli, Thibault Fevry, Jason Alan Fries, Ryan Teehan, Stella Biderman, Leo Gao, Tali Bers, Thomas Wolf, and Alexander M Rush.  Multitask prompted training enables Zero-Shot task generalization.  _arXiv_ , October 2021.
  * Levin et al. [2022] Roman Levin, Valeriia Cherepanova, Avi Schwarzschild, Arpit Bansal, C Bayan Bruss, Tom Goldstein, Andrew Gordon Wilson, and Micah Goldblum.  Transfer learning with deep tabular models.  _arXiv [cs.LG]_ , June 2022.
  * Huang et al. [2020] Xin Huang, Ashish Khetan, Milan Cvitkovic, and Zohar Karnin.  TabTransformer: Tabular data modeling using contextual embeddings.  _arXiv_ , December 2020.
  * Gorishniy et al. [2021] Yury Gorishniy, Ivan Rubachev, Valentin Khrulkov, and Artem Babenko.  Revisiting deep learning models for tabular data.  _arXiv_ , June 2021.
  * Grzyb et al. [2021] Mateusz Grzyb, Zuzanna Trafas, Katarzyna Woznica, and Przemyslaw Biecek.  metamimic: analysis of hyperparameter transferability for tabular data using mimic-iv database.  https://github.com/ModelOriented/metaMIMIC, 2021.  Accessed: 2023-1-29.
  * Wang et al. [2019] Jie Wang, Xingyu Jiang, and Xiaodong He.  Tabular data representation learning with graph convolutional networks.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 33:7712–7719, 2019.
  * Chen et al. [2020] Yufan Chen, Zhe Shao, Zhiyuan Liu, and Xiaodan Zhu.  Graph convolutional networks for multi-omics data integration and disease subtype identification.  _Bioinformatics_ , 36(8):2820–2828, 2020.
  * He et al. [2021b] Pengcheng He, Jianfeng Gao, and Weizhu Chen.  Debertav3: Improving deberta using electra-style pre-training with gradient-disentangled embedding sharing, 2021b.
  * McMaster et al. [2023] Christopher McMaster, Julia Chan, David F L Liew, Elizabeth Su, Albert G Frauman, Wendy W Chapman, and Douglas E V Pires.  Developing a deep learning natural language processing algorithm for automated reporting of adverse drug reactions.  _Journal of biomedical informatics_ , 137:104265, January 2023.  ISSN 1532-0464, 1532-0480.  doi: 10.1016/j.jbi.2022.104265.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _arXiv_ , June 2017.
  * Johnson et al. [2020a] Alistair Johnson, Tom Pollard, and Roger Mark.  MIMIC-III clinical database, 2020a.
  * Johnson et al. [2016] Alistair E W Johnson, Tom J Pollard, Lu Shen, Li-Wei H Lehman, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark.  MIMIC-III, a freely accessible critical care database.  _Scientific data_ , 3:160035, May 2016.  ISSN 2052-4463.  doi: 10.1038/sdata.2016.35.
  * [23] PubMed annual baseline.  https://ftp.ncbi.nlm.nih.gov/pubmed/baseline/.  Accessed: 2022-12-10.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939785.  URL http://doi.acm.org/10.1145/2939672.2939785.
  * Goldberger et al. [2000] A L Goldberger, L A Amaral, L Glass, J M Hausdorff, P C Ivanov, R G Mark, J E Mietus, G B Moody, C K Peng, and H E Stanley.  PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals.  _Circulation_ , 101(23):E215–20, June 2000.  ISSN 0009-7322, 1524-4539.  doi: 10.1161/01.cir.101.23.e215.
  * Johnson et al. [2020b] Alistair Johnson, Lucas Bulgarelli, Tom Pollard, Steven Horng, Leo Anthony Celi, and Roger Mark.  MIMIC-IV, 2020b.
  * Johnson et al. [2022] Alistair Johnson, Lucas Bulgarelli, Tom Pollard, Leo Anthony Celi, Roger Mark, and Steven Horng.  MIMIC-IV-ED, 2022.
  * Xie et al. [2022] Feng Xie, Jun Zhou, Jin Wee Lee, Mingrui Tan, Siqi Li, Logasan S/o Rajnthern, Marcel Lucas Chee, Bibhas Chakraborty, An-Kwok Ian Wong, Alon Dagan, Marcus Eng Hock Ong, Fei Gao, and Nan Liu.  Benchmarking emergency department prediction models with machine learning and public electronic health records.  _Scientific data_ , 9(1):658, October 2022.  ISSN 2052-4463.  doi: 10.1038/s41597-022-01782-9.
  * Sundararajan et al. [2017] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  _arXiv_ , March 2017.
  * Mudrakarta et al. [2018] Pramod Kaushik Mudrakarta, Ankur Taly, Mukund Sundararajan, and Kedar Dhamdhere.  Did the model understand the question?  _arXiv_ , May 2018.
  * Kokhlikyan et al. [2020] Narine Kokhlikyan, Vivek Miglani, Miguel Martin, Edward Wang, Bilal Alsallakh, Jonathan Reynolds, Alexander Melnikov, Natalia Kliushkina, Carlos Araya, Siqi Yan, and Orion Reblitz-Richardson.  Captum: A unified and generic model interpretability library for pytorch, 2020.