# Time-Resolved Reconstruction of Motion, Force, and Stiffness using Spectro-
Dynamic MRI

Max H. C. van Riel, Tristan van Leeuwen, Cornelis A. T. van den Berg,
Alessandro Sbrizzi This work was supported by a grant by the Nederlandse
Organisatie voor Wetenschappelijk Onderzoek (NWO) under grant number 18897.
(Corresponding author: Max H. C. van Riel)M. H. C. van Riel, C. A. T. van den
Berg and A. Sbrizzi are with the Department of Radiotherapy, Computational
Imaging Group for MR Diagnostics and Therapy, UMC Utrecht, Utrecht, The
Netherlands (e-mail: m.h.c.vanriel-3@umcutrecht.nl).T. van Leeuwen is with the
Mathematical Institute, Utrecht University, Utrecht, The Netherlands.© 2023
IEEE. Personal use of this material is permitted. Permission from IEEE must be
obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

###### Abstract

Measuring the dynamics and mechanical properties of muscles and joints is
important to understand the (patho)physiology of muscles. However, acquiring
dynamic time-resolved MRI data is challenging. We have previously developed
Spectro-Dynamic MRI which allows the characterization of dynamical systems at
a high spatial and temporal resolution directly from k-space data. This work
presents an extended Spectro-Dynamic MRI framework that reconstructs 1) time-
resolved MR images, 2) time-resolved motion fields, 3) dynamical parameters,
and 4) an activation force, at a temporal resolution of 11 ms. An iterative
algorithm solves a minimization problem containing four terms: a motion model
relating the motion to the fully-sampled k-space data, a dynamical model
describing the expected type of dynamics, a data consistency term describing
the undersampling pattern, and finally a regularization term for the
activation force. We acquired MRI data using a dynamic motion phantom
programmed to move like an actively driven linear elastic system, from which
all dynamic variables could be accurately reconstructed, regardless of the
sampling pattern. The proposed method performed better than a two-step
approach, where time-resolved images were first reconstructed from the
undersampled data without any information about the motion, followed by a
motion estimation step.

###### Index Terms:

Dynamical systems, magnetic resonance imaging, real-time imaging, spectro-
dynamic MRI, time-resolved imaging

##  I Introduction

Dynamic measurements of joints and muscles are important to study the
(patho)physiology of the musculoskeletal system [1, 2, 3, 4, 5]. Previous
studies have shown that dynamic scans result in different kinematics compared
to static scans [6]. Mechanical properties of biological tissues can give
additional quantitative information, but these values are often obtained ex
vivo or (quasi-)statically [7, 8]. MRI is a promising imaging modality for
measuring these dynamic properties in vivo, as it provides excellent soft-
tissue contrast. However, acquiring dynamic MRI data at a high spatial and
temporal resolution remains a challenge.

Several techniques exist for dynamic musculoskeletal imaging using MRI. Some
methods use gating to bin the data into different motion states [9, 10]. In
contrast, real-time MRI does not require periodic motion or synchronization of
the acquisition to the motion [11, 12]. Sparsity along the temporal dimension
can be exploited to reconstruct dynamic musculoskeletal images from
undersampled data [13, 14, 15]. Real-time methods are preferred for
musculoskeletal applications, since accurately repeating the same motion is
difficult to achieve, and sometimes not feasible for patients who experience
pain. Therefore, the development of time-resolved 3D MRI techniques with a
high temporal resolution is an active area of research [16].

All methods described so far extract dynamic information from a time series of
images. However, to achieve a sufficient temporal resolution, high
undersampling factors are required. Any undersampling artifacts remaining in
the images will have a detrimental effect on the estimated dynamics. It has
been shown that dynamic information can be extracted directly from k-space
data, even in the case of very high undersampling factors, without an
intermediate image reconstruction step [17].

In our previous work [18], we have proposed the Spectro-Dynamic MRI framework
in an effort to characterize dynamical systems at a high temporal resolution.
This novel acquisition paradigm uses a spectral motion model, which relates
the raw MRI measurements in the spectral domain (k-space) to the displacement
field. Thus, it can reconstruct the motion directly from k-space data, without
requiring images at a high temporal resolution as an intermediate step. In
addition, it uses a dynamical model which adds prior knowledge about the
structure of the motion. This dynamical model also allows for the estimation
of dynamical parameters from extremely undersampled data. We demonstrated
Spectro-Dynamic MRI using two coupled spherical pendula, from which the
motion, the length of the pendula, and the spring stiffness could be estimated
accurately.

In this work, we address some of the limitations of the Spectro-Dynamic MRI
reconstruction as presented in [18]. Firstly, the reconstruction method
required rigid motion fields, or at most motion fields that are spatially
varying along the readout direction, due to the high undersampling factor.
This prevented the reconstruction of more general motion fields. Secondly, the
measurement noise was not explicitly taken into account in the reconstruction.
This could lead to biases in the estimated motion.

Here, we propose an extended iterative reconstruction framework for Spectro-
Dynamic MRI to overcome these limitations. Besides the motion field and
dynamical parameters, the full time-resolved k-space is added as an additional
variable in the reconstruction. By jointly reconstructing the missing k-space
data and the motion fields in a single optimization problem, a very high
undersampling factor per time instance could be achieved. In addition, we have
added an external activation force, which allows for the dynamical model to
display the output of an arbitrary forcing function. This activation
represents an external load acting on a muscle, or neuronal activation of a
muscle. Using a dynamic motion phantom, we show that we can reconstruct 1)
time-resolved images, 2) time-resolved motion fields, 3) dynamical parameters,
and 4) the activation force, at a temporal resolution of 11 ms and a spatial
resolution of 5.0 mm $\times$ 5.0 mm. These results bring us closer to our
goal of real-time in vivo dynamical system characterization.

##  II Theory

TABLE I: Table of variables. Symbol | Description
---|---
bold | Vector quantity
$\hat{}$ | Spectral quantity
$\mathbf{r},\mathbf{k}$ | Spatial/Spectral coordinate
$\Omega_{M},\Omega_{S}$ | Moving/Stationary compartment
$d$ | Number of dimensions
$P$ | Number of basis functions
$N$ | Number of k-space samples per time instance
$M$ | Number of measured samples per time instance
$T$ | Number of time instances
$m,\hat{m}$ | Time-resolved image/k-space data
$\mathbf{v},\mathbf{u}$ | Velocity/Displacement field
$\mathcal{F}$ | Fourier transform operator
$\mathbf{d}$ | Undersampled k-space data
$\boldsymbol{\phi}$ | Basis function of the displacement field
$\mathbf{m}$ | Time-resolved k-space
$\mathbf{q}$ | Motion field coefficients
$\kappa$ | Elastic stiffness
$c$ | Damping coefficient
$\mathbf{f}$ | Activation force
$G$ | Motion model
$F$ | Dynamical model
$H$ | Data consistency model
$R$ | Regularization
$\lambda_{F},\lambda_{H},\lambda_{R}$ | Regularization parameters
$\alpha$ | Readout rotation angle

The extended Spectro-Dynamic MRI reconstruction framework contains four
components (Fig. 1): a motion model ($G$), a dynamical model ($F$), a data
consistency term ($H$), and a regularization term ($R$). These four components
are combined as penalty terms in a minimization problem. We use penalty terms
instead of hard equality constraints to capture model inaccuracies and noise
in all four components [19]. Each component is discussed in more detail in the
following subsections. The minimization problem is solved with an iterative
reconstruction algorithm, which will be described in Section II-E. Some
symbols used throughout this work are listed in Table I.

###  II-A Motion Model

The motion model $G$ provides the relation between the MRI data and the motion
field. Its derivation starts by assuming that all signal is conserved during
motion [18]. This assumption holds if the magnetization is in steady state,
the receive and excitation fields are homogeneous, and the readouts are short.

Let $\Omega_{r}\subseteq\mathbb{R}^{d}$ be a spatial domain in $d$ dimensions,
and let $\Omega_{t}\subseteq\mathbb{R}$ be a time interval. Then
$\Omega_{rt}=\Omega_{r}\times\Omega_{t}$ is the combined spatio-temporal
domain. In this work, only 2D experiments are shown, but all methods are valid
in 3D as well. The extension to 3D will be the topic of future research.

Let $m(\mathbf{r},t):\Omega_{rt}\to\mathbb{C}$ be a temporal series of
complex-valued images, and let
$\mathbf{v}(\mathbf{r},t):\Omega_{rt}\to\mathbb{R}^{d}$ be the $d$-dimensional
velocity field. The conservation of magnetization results in the following
partial differential equation (PDE), known as the continuity equation:

$\frac{\partial m}{\partial t}+\nabla\cdot\left(m\mathbf{v}\right)=0.$ (1)

Note that combined with the incompressibility constraint
($\nabla\cdot\mathbf{v}=0$), the motion model becomes identical to the optical
flow model [20] as used in image registration.

MRI measurements are acquired in the spectral or spatial frequency domain,
called k-space. Therefore, it is convenient to transform the continuity
equation (1) to the spectral domain. Let $\mathcal{F}$ be the $d$-dimensional
spatial Fourier transform operator. It transforms $m$ into its k-space
representation $\hat{m}(\mathbf{k},t):\Omega_{kt}\to\mathbb{C}$, such that
$\hat{m}=\mathcal{F}m$ and $m=\mathcal{F}^{-1}\hat{m}$. Here,
$\Omega_{k}\subseteq\mathbb{R}^{d}$ is the spectral domain in k-space, and
${\Omega_{kt}=\Omega_{k}\times\Omega_{t}}$. Similarly, the spectral velocity
field $\hat{\mathbf{v}}(\mathbf{k},t):\Omega_{kt}\to\mathbb{C}^{d}$ is
obtained by applying $\mathcal{F}$ to each component of $\mathbf{v}$, such
that $\hat{v}_{j}=\mathcal{F}v_{j}$ for $j=1,\cdots,d$. Note that an
undersampled and noisy version of $\hat{m}$ is measured during an MRI scan.

Figure 1: Overview of the different elements used in the iterative
reconstruction algorithm. The four different model components are indicated in
the middle column. The measured undersampled data ($\mathbf{d}$) and the
chosen basis functions ($\boldsymbol{\phi}$) are the inputs of the
reconstruction, and remain constant during the reconstruction. The
optimization variables that are reconstructed are given on the right. The
connections indicate which inputs and variables are used in each model. Note
that every variable is shared between at least two models.

Using the properties of the Fourier transform, (1) can be converted to the
spectral domain:

$\frac{\partial\hat{m}}{\partial t}+2\pi
i\sum_{j=1}^{d}k_{j}(\hat{m}*\hat{v}_{j})=0.$ (2)

Here, $i=\sqrt{-1}$ is the imaginary unit, $k_{j}$ is the spectral coordinate
in dimension $j$, and the symbol $*$ is a convolution over the $d$ spectral
dimensions.

The equation above is written in terms of the velocity field, while the
dynamical model (Section II-B) is described in terms of the displacement
field. Let $\mathbf{u}(\mathbf{r},t):\Omega_{rt}\to\mathbb{R}^{d}$ be the
$d$-dimensional displacement vector field. We will parameterize the
displacement field using $P$ basis functions
$\boldsymbol{\phi}_{p}(\mathbf{r}):\Omega_{r}\to\mathbb{R}^{d}$ and
corresponding generalized coordinates $q_{p}(t):\Omega_{t}\to\mathbb{R}$ for
$p=1,\cdots,P$:

$\mathbf{u}(\mathbf{r},t)=\sum_{p=1}^{P}\boldsymbol{\phi}_{p}(\mathbf{r})q_{p}(t).$
(3)

By using the basis functions $\boldsymbol{\phi}_{p}$, prior information about
the motion field can be introduced, thereby reducing the number of unknown
variables. For example, if the motion field is smooth, splines can be used as
basis functions. Alternatively, the basis functions can be piecewise constant
if the domain $\Omega_{r}$ can be separated into different compartments. The
latter approach will be used in this work (Section III-D and Fig. 5).

The velocity field and the displacement field are related through the material
derivative
(${\mathbf{v}=\frac{D\mathbf{u}}{Dt}=\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{u}}$). In this work we will limit ourselves to
piecewise constant basis functions. Therefore, $\nabla u_{j}=\mathbf{0}$ for
all $j$, and thus $\mathbf{v}=\frac{\partial\mathbf{u}}{\partial t}$.
Combining (2) and (3) gives:

$\frac{\partial\hat{m}}{\partial t}+2\pi
i\sum_{j=1}^{d}\sum_{p=1}^{P}k_{j}(\hat{m}*\hat{\phi}_{p,j)}\frac{dq_{p}}{dt}=0,$
(4)

where $\hat{\phi}_{p,j}=\mathcal{F}\phi_{p,j}$.

During an MRI scan, the continuous time-resolved k-space $\hat{m}$ is sampled
at a discrete number of sample points. The k-space domain $\Omega_{k}$ is
discretized to $N$ points, and $\Omega_{t}$ is discretized to $T$ time
instances. When all discrete points of $\hat{m}$ are put into one vector, we
get $\mathbf{m}\in\mathbb{C}^{NT}$. Similarly, we can discretize the
generalized coordinates over time to get $\mathbf{q}\in\mathbb{R}^{PT}$. Note
that $\mathbf{m}$ cannot be measured directly at a high spatio-temporal
resolution. Instead, a highly undersampled and noisy version of $\mathbf{m}$
is acquired (see Section II-C).

After discretization, the motion model term $G$ is the squared $L^{2}$-norm of
the residual of (4), which is a function of $\mathbf{m}$ and $\mathbf{q}$:

$G(\mathbf{m},\mathbf{q})=\frac{1}{2}\bigl{\|}D_{t}\mathbf{m}+2\pi
i\sum_{j=1}^{d}K_{j}C(\mathbf{m},D_{t}\Phi_{j}\mathbf{q})\bigr{\|}_{2}^{2}.$
(5)

Here, $D_{t}$ is the first-order temporal finite difference operator, $K_{j}$
is a diagonal matrix that performs multiplication with the k-space
coefficients $k_{j}$, $C$ is a bilinear map that performs the spatial
convolution, and $\Phi_{j}$ is a block diagonal matrix whose columns contain
the spectral basis functions $\hat{\phi}_{p,j}$.

Note that the term inside the norm in (5) is bilinear in $\mathbf{m}$ and
$\mathbf{q}$, making $G(\mathbf{m},\mathbf{q})$ a biconvex function.

###  II-B Dynamical Model

If the measured data were fully sampled and free of noise, the coefficients of
the motion field could be retrieved directly from the measured data. However,
a fully-sampled k-space at a high spatial and temporal resolution is not
feasible.

To overcome this issue, the dynamical model is added to constrain the solution
of the motion field. This model provides additional information about the
expected type of motion. In general, the dynamical model is a PDE constructed
from the balance of momentum (Newton’s second law) combined with a material’s
constitutive relation, like Hooke’s law for linear elasticity. This
constitutive model contains the dynamical parameters, which can either be
provided a priori or can be estimated from the data to provide additional
information. Examples of dynamical parameters are stiffness and viscosity.

For now, we will limit ourselves to dynamical models that are only time-
dependent, and can thus be written as an ordinary differential equation (ODE).
All spatial dependencies are captured by the basis functions
$\boldsymbol{\phi}_{p}$. In this work, we use a linear elastic model with an
external activation for each generalized coordinate $q_{p}$:

$\frac{d^{2}q_{p}}{dt^{2}}+c\frac{dq_{p}}{dt}+\kappa q_{p}=f_{p}\qquad\forall
p,$ (6)

with $c\in\mathbb{R}_{+}$ the damping coefficient in Ns/m,
$\kappa\in\mathbb{R}_{+}$ the elastic stiffness in N/m, and
${f_{p}(t):\Omega_{t}\to\mathbb{R}}$ the activation force in N for the $p$-th
degree of freedom.

In our experiments, we set $c$ to a particular value, which is fixed during
the reconstruction. Also note the absence of the mass in the first term of
(6). Since the equation can be scaled by an arbitrary constant, we regard all
quantities as mass-normalized. Thus, $\kappa$ remains as the only unknown
dynamical parameter.

Despite its simplicity, this model can be used to describe muscle dynamics in
certain tasks [21], or it can act as a linear approximation to more complex
dynamical systems. In this case, $f_{p}$ can be interpreted as the internal
muscular activation caused by a neuronal input, or as an external loading,
depending on the application.

The dynamical model term $F$ takes the residual of (6) and takes its squared
$L^{2}$-norm summed over the degrees of freedom. After discretization, this
can be written as a function of $\mathbf{q}$, $\kappa$, and $\mathbf{f}$:

$F(\mathbf{q},\kappa,\mathbf{f})=\frac{1}{2}\bigl{\|}\left(D_{tt}+cD_{t}+\kappa
I\right)\mathbf{q}-\mathbf{f}\bigr{\|}_{2}^{2}.$ (7)

$D_{t}$ and $D_{tt}$ are the first- and second-order temporal finite
difference operators, $I$ is the $PT\times PT$ identity operator, and
$\mathbf{f}\in\mathbb{R}^{PT}$ is the discretized activation force.

The term inside the norm of (7) is linear in $\mathbf{q}$ and $\mathbf{f}$ for
a fixed value of $\kappa$, and it is linear in $\kappa$ and $\mathbf{f}$ for a
fixed value of $\mathbf{q}$. This makes $F(\mathbf{q},\kappa,\mathbf{f})$,
like $G(\mathbf{m},\mathbf{q})$ in (5), a biconvex function.

###  II-C Data Consistency

Sampling k-space during an MRI scan is not instantaneous. If all $N$ samples
are to be acquired at each time instance, the time step $\Delta t$ between two
time instances becomes too large to capture motion at a sub-second timescale.
Therefore, the measured k-space data $\mathbf{d}\in\mathbb{C}^{MT}$ will be an
undersampled version of $\mathbf{m}$, with $M$ measured samples per time
instance, and $M$ much smaller than the number of samples in the fully-sampled
data ($M\ll N$). In practice, $M$ can be the number of samples during a single
repetition time (TR), or multiple readouts can be grouped together if the
repetition time is small enough. The sampling pattern used during the
experiment is given by $E\in\\{0,1\\}^{MT\times NT}$, such that
$E\mathbf{m}=\mathbf{d}$ in the absence of noise.

In our previous work [18], we used the measurements $\mathbf{d}$ directly in
the motion model. However, this created two problems. First of all, the motion
model contains a convolution between the k-space data and the velocity field.
This convolution can only be evaluated when the k-space is fully sampled at
each time instance. In practice, the k-space of each time instance will be
heavily undersampled if a high temporal resolution is required. The data
points that are not sampled leave gaps in k-space, while these data points are
required to calculate the convolution. In [18], this issue was circumvented by
only considering motion fields that vary spatially in the readout direction
and are constant in all other directions, resulting in a 1D convolution in the
fully-sampled readout direction. However, for more general motion fields or
more complex acquisition trajectories, this approach is no longer possible.

Furthermore, $\mathbf{d}$ will invariably contain measurement noise. By
plugging $\mathbf{d}$ directly into the motion model, this noise gets scaled
by the different terms in (5). For example, the high spatial frequency noise
is amplified through the operator $K_{j}$. This will result in a biased
reconstruction of the motion field.

The first limitation is addressed by including the full time-resolved k-space,
including the data points that were not sampled, as an additional variable in
the reconstruction. This increases the number of unknown variables, but gives
the reconstruction more flexibility. By using a specific interleaved sampling
pattern, the information coming from the sampled data at a certain time
instance is shared through the temporal derivatives in the motion model to the
time instances where those data points are not sampled. That way, the gaps in
the time-resolved k-space can be filled in, allowing the convolution to be
evaluated for more general motion patterns.

Instead of using a hard data consistency constraint in the form of
$E\mathbf{m}-\mathbf{d}=\mathbf{0}$, we will use a relaxation term in the
reconstruction that minimizes the error of this equation:

$H(\mathbf{m})=\frac{1}{2}\|E\mathbf{m}-\mathbf{d}\|_{2}^{2}.$ (8)

This relaxation term allows the reconstructed data $\mathbf{m}$ to deviate
slightly from the measured data $\mathbf{d}$. We expect that this deviation
can capture the measurement noise and other measurement imperfections, thereby
reducing the bias in the motion field caused by these differences.

Figure 2: Simulated inputs for the continuous (left) and discontinuous (right)
experiments. On the top row, the activation forces are shown for both cases.
The middle and bottom rows show the displacements and velocities respectively,
as obtained by numerically solving (6). Only the two small tubes moved along a
single direction, while the water compartment and the large cylinder remained
stationary.

Figure 3: (a) The Quasar MRI4D Motion Phantom. A large water compartment
surrounds a stationary gel-filled tube and a cylinder that is actuated by an
MR-compatible piezoelectric motor. Inside the moving cylinder, two smaller
gel-filled tubes were placed. (b) Stationary MR images of the phantom setup,
acquired with the same scan parameters as those used during the dynamic
experiments. The red and blue arrows indicate the readout and phase encode
directions respectively for the three different readout orientations
($\alpha=$ 0°, 45°, and 90°). These directions correspond to the $x$\- and
$y$-directions respectively in the reconstructed motion. (c) A stationary MR
image, where the two yellow lines indicate the division between the moving
compartment ($\Omega_{M}$) and the stationary compartment ($\Omega_{S}$).

###  II-D Regularization

By treating $\mathbf{q}$, $\kappa$ and $\mathbf{f}$ as unknown variables, we
have created an ill-posed reconstruction problem. Any error in the dynamical
model given by (6) can be compensated by changing its right-hand side,
resulting in unfeasible solutions for $\mathbf{f}$.

Therefore, we added a fourth term $R$ to the reconstruction. This term acted
as a regularization of the activation force and depends on the prior knowledge
that is available about the shape of the activation force. For example, if
$f_{p}$ is continuous, a smoothness regularization can be used. Alternatively,
if $f_{p}$ is known to be piecewise constant, a regularization function based
on the temporal Total Variation (TV) is more suitable. Both regularization
terms can be discretized using the second- and first-order finite difference
operators $D_{tt}$ and $D_{t}$ respectively, resulting in:

$R_{\text{S}}(\mathbf{f})=\frac{1}{2}\|D_{tt}\mathbf{f}\|_{2}^{2},$ (9)

$R_{\text{TV}}(\mathbf{f})=\|D_{t}\mathbf{f}\|_{1}.$ (10)

###  II-E Iterative Reconstruction

Combining all four components, the extended Spectro-Dynamic MRI framework aims
to solve the following optimization problem:

$\min_{\mathbf{m},\mathbf{q},\kappa,\mathbf{f}}G(\mathbf{m},\mathbf{q})+\lambda_{F}F(\mathbf{q},\kappa,\mathbf{f})+\lambda_{H}H(\mathbf{m})+\lambda_{R}R(\mathbf{f}).$
(11)

The parameters $\lambda_{F}$, $\lambda_{H}$, and $\lambda_{R}$ determine the
trade-off between the residuals of the different models.

The minimization in (11) must be solved for four different parameters:
$\mathbf{m}$, $\mathbf{q}$, $\kappa$, and $\mathbf{f}$, for a total of
${NT+2PT+1}$ unknowns. It is a multiconvex optimization problem, for which
several optimization algorithms exist. In this work, we solve (11) with an
iterative scheme called Block Coordinate Descent (BCD) [22] because it is easy
to implement and interpret the results. In every step of BCD, one subset of
the unknown variables is fixed, while problem is solved for the remaining
variables using convex optimization. This is done most efficiently by fixing a
minimal number of variables [23]. In the case of (11), these minimal fixed
sets are $\\{\mathbf{m},\kappa\\}$ and $\\{\mathbf{u}\\}$. The resulting
iterative optimization is given by Algorithm 1.

Figure 4: Sampling pattern as used during the dynamic experiments. The readout
direction is $k_{x}$, while $k_{y}$ is the phase encode direction. Only eight
of the 64 phase encode lines are visualized. Each time instance contains two
readout lines, for an effective temporal resolution of two times the
repetition time (TR). By distributing the readout lines evenly over the
k-space, information about different spatial frequencies is acquired during
each time instance. Algorithm 1 Spectro-Dynamic MRI Reconstruction.

1: $\mathbf{m}^{0},\mathbf{q}^{0},\kappa^{0},\mathbf{f}^{0}\leftarrow 0$

2: for $k=1$ to $K$ do

3:
$\displaystyle\mathbf{m}^{k}\leftarrow\operatorname*{arg\,min}_{\mathbf{m}}G(\mathbf{m},\mathbf{q}^{k-1})+\lambda_{H}H(\mathbf{m})$

4: $\begin{aligned}
\displaystyle\mathbf{q}^{k},\bar{\mathbf{f}}^{k}\leftarrow\operatorname*{arg\,min}_{\mathbf{q},\mathbf{f}}\Bigl{[}&G(\mathbf{m}^{k},\mathbf{q})+\lambda_{F}F(\mathbf{q},\kappa^{k-1},\mathbf{f})\\\
&+\lambda_{R}R(\mathbf{f})\Bigr{]}\end{aligned}$

5:
$\displaystyle\kappa^{k},\mathbf{f}^{k}\leftarrow\operatorname*{arg\,min}_{\kappa,\mathbf{f}}\lambda_{F}F(\mathbf{q}^{k},\kappa,\mathbf{f})+\lambda_{R}R(\mathbf{f})$

6: end for

Note that line 4 in Algorithm 1 is a convex optimization with $\mathbf{m}$ and
$\kappa$ fixed, while lines 3 and 5 together form a convex optimization with
$\mathbf{q}$ fixed. This last minimization problem over $\mathbf{m}$, $\kappa$
and $\mathbf{f}$ is separable into two independent minimization problems: one
over $\mathbf{m}$ (line 3), and one over $\kappa$ and $\mathbf{f}$ (line 5).
The ordering of the subproblems is chosen such that the influence of the
initial values is minimal. The activation force $\mathbf{f}$ is optimized
twice in each iteration, therefore $\bar{\mathbf{f}}^{k}$ is a nuisance
variable only used to improve the solution of $\mathbf{u}$.

When the smooth regularization function (9) is used, all three subproblems
(lines 3-5 in Algorithm 1) can be solved using linear least squares solvers.
In the case of the Total Variation regularizer (10), lines 4 and 5 are still
convex, but no longer smooth. In this case, these are solved using the
Alternating Direction Method of Multipliers (ADMM) [24, pp. 43–44].

Figure 5: The four basis functions used during the reconstruction for
$\alpha=0\degree{}$. Two basis functions are defined in the moving compartment
($\Omega_{M}$, top row), the other two in the stationary compartment
($\Omega_{S}$, bottom row). For each compartment, one basis function points in
the vertical direction ($x$, left column), the other in the horizontal
direction ($y$, right column). The coefficients of all four basis functions
are estimated during the reconstruction, but only those in the moving
compartment should be nonzero.

##  III Methods

###  III-A Activation Functions

We used two different kinds of activations to observe different kinds of
dynamics in our system. First, a smooth continuous activation function was
constructed by adding two sinusoidal functions with frequencies of 0.15 Hz and
0.33 Hz (see Fig. 2, top left). For this activation, the reconstruction was
regularized by $R_{\text{S}}$. The next experiment used a discontinuous
activation, where the activation switched between an ’off’ and an ’on’ state
(Fig. 2, top right). For this experiment, $R_{\text{TV}}$ was used as a
regularization function.

The differential equation of motion (6) was solved using $\kappa=\textrm{30
N/m}$ and the initial condition ${q(0)=0}$, where $f(t)$ is given by one of
the activations in Fig. 2. We used a damping coefficient of $c=\textrm{0
Ns/m}$ for the continuous activation, and $c=\textrm{1 Ns/m}$ for the
discontinuous activation. This value was fixed during the reconstruction.

###  III-B Phantom Setup

We used the Quasar MRI4D Motion Phantom (Modus Medical Devices Inc., London,
ON, Canada) to acquire dynamic experimental data. This phantom (Fig. 3(a))
contains an MR-compatible piezoelectric motor, which drives a cylindrical
moving compartment inside a stationary compartment filled with water. Inside
the moving cylinder, two gel-filled tubes (TO5, Eurospin II test system,
Scotland) were placed. A second, static cylinder filled with gel was inserted
next to the moving cylinder. The motion fields were piecewise constant, which
means that the motion was locally rigid. As such, they can be described by
four degrees of freedom (two for the moving compartment and two for the
stationary compartment). The simulated motion profile (Fig. 2, middle row) was
used as position setpoint in the phantom’s motion control software.

###  III-C Data Acquisition

Figure 6: Images reconstructed using (a) the two-step reconstruction, (b) the
result of the Spectro-Dynamic MRI reconstruction, and (c) stationary data. The
readout direction was parallel to the direction of motion
($\alpha=0\degree{}$). The bottom row shows one line in the image (indicated
by the dashed yellow line) over time. Note the removal of the motion blurring
in the Spectro-Dynamic MRI reconstruction. A video showing all frames of the
time-resolved Spectro-Dynamic MRI reconstruction is available online as
Supplementary Video S1.

The motion phantom was placed in a 1.5T MRI scanner (Ingenia, Philips
Healthcare, Best, The Netherlands). Both stationary and dynamic data were
acquired with a spoiled gradient echo sequence (TR $=$ 5.5 ms, TE $=$ 2.2 ms,
Flip angle $=$ 9°). The field of view was a 320 mm $\times$ 320 mm coronal
slice with a thickness of 5 mm. The phantom was leveled to prevent through-
slice motion. To achieve a homogeneous receive sensitivity, the body coil was
used for data acquisition. Data sampling was performed on a 64 $\times$ 64
matrix using an interleaved pattern (Fig. 4). Thus, the acquired information
at each time instance was distributed over k-space. The sampling pattern was
repeated 40 times, for a total of 2560 readouts, and a total acquisition time
of 14 seconds.

To show that the extended Spectro-Dynamic MRI framework can estimate motion
along an arbitrary direction, the experiment was repeated three times for both
activations. Each time, the readout direction ($k_{x}$) was rotated
respectively by 0°, 45°, and 90° with respect to the direction of motion (Fig.
3(b)).

###  III-D Reconstruction

Four degrees of freedom ($P=\textrm{4}$) were used during reconstruction, for
two compartments (the moving compartment $\Omega_{M}$ and the stationary
compartment $\Omega_{S}$, Fig. 3(c)) and two directions ($x$ and $y$) per
compartment. Thus, four piecewise constant basis functions were defined to
parameterize the displacement vector field according to (3), corresponding
with four generalized coordinates (Fig. 5):

$\displaystyle\boldsymbol{\phi}_{1}(\mathbf{r})$
$\displaystyle=\begin{cases}(1,0)^{T}\qquad\forall\mathbf{r}\in\Omega_{M}\\\
(0,0)^{T}\qquad\forall\mathbf{r}\in\Omega_{S}\end{cases}$ (12)
$\displaystyle\boldsymbol{\phi}_{2}(\mathbf{r})$
$\displaystyle=\begin{cases}(0,1)^{T}\qquad\forall\mathbf{r}\in\Omega_{M}\\\
(0,0)^{T}\qquad\forall\mathbf{r}\in\Omega_{S}\end{cases}$
$\displaystyle\boldsymbol{\phi}_{3}(\mathbf{r})$
$\displaystyle=\begin{cases}(0,0)^{T}\qquad\forall\mathbf{r}\in\Omega_{M}\\\
(1,0)^{T}\qquad\forall\mathbf{r}\in\Omega_{S}\end{cases}$
$\displaystyle\boldsymbol{\phi}_{4}(\mathbf{r})$
$\displaystyle=\begin{cases}(0,0)^{T}\qquad\forall\mathbf{r}\in\Omega_{M}\\\
(0,1)^{T}\qquad\forall\mathbf{r}\in\Omega_{S}\end{cases}$

Since the moving compartment moves along a straight line, the reconstructed
displacements of only one degree of freedom should be nonzero for $\alpha=$ 0°
and 90°. For $\alpha=$ 45°, the displacements of two degrees of freedom should
be nonzero and identical.

Multiple successive readouts were grouped together during reconstruction. This
was done to increase the amount of available data per time instance while
guaranteeing a high temporal resolution. Including more lines in a single time
instance reduces the undersampling of k-space, at the cost of a lower temporal
resolution. We grouped together every two successive readouts (Fig. 4). This
resulted in an effective temporal resolution of $\Delta
t=\textrm{2}\cdot\textrm{TR}=\textrm{11 ms}$ and an undersampling factor of 32
for each time instance. This temporal resolution enables future applications
in real-time dynamic cardiac or speech imaging, which require similar frame
rates [12].

The reconstruction, as described by Algorithm 1, was implemented in Matlab
(The MathWorks Inc., Natick, MA, USA). The large-scale and sparse subproblem
in line 3 of Algorithm 1 was solved using the iterative linear least-squares
solver lsqr. The other two subproblems were solved using Matlab’s mldivide
(\\) or an ADMM implementation, depending on $R(\mathbf{f})$. The
regularization parameters
(${\lambda_{F}=\textrm{5.0}\cdot\textrm{10\textsuperscript{6}}}$,
${\lambda_{H}=\textrm{1.0}\cdot\textrm{10\textsuperscript{3}}}$, and
${\lambda_{R}=\textrm{1.0}\cdot\textrm{10\textsuperscript{3} /
2.0}\cdot\textrm{10\textsuperscript{4}}}$ for $R_{\text{S}}\textrm{ /
}R_{\text{TV}}$) were determined empirically and fixed for all
reconstructions. These values resulted in comparable contributions of all four
model terms to the objective function. After 15 iterations, the updates to the
estimated values were negligible (less than 1%), and the algorithm was
terminated. Each reconstruction took approximately 2 hours on a workstation
with a 2.90 GHz Intel Xeon W-2102 CPU and 32GB RAM.

The estimated fully-sampled k-space $\mathbf{m}$ was transformed into the
image domain with a Fast Fourier Transform (FFT) followed by a geometry
correction step to correct for gradient nonlinearities. The velocity was
calculated using central finite differences from the estimated displacements.
Together, the activation force and the velocity are important for
understanding muscle physiology [5]. For the estimated displacements,
velocities, and activation forces, the root-mean-square error (RMSE) compared
to the ground truth values (Fig. 2) was calculated as a measure for the
reconstruction accuracy.

###  III-E Two-step Reconstruction

We compared the proposed Spectro-Dynamic MRI reconstruction to a two-step
approach. In the first step, the time-resolved k-space $\mathbf{m}$ was
reconstructed by only using a temporal smoothness regularization and no
information about the motion fields:

$\min_{\mathbf{m}}H(\mathbf{m})+\lambda\|D_{t}F^{H}\mathbf{m}\|_{1},$ (13)

with $F^{H}$ the inverse Fourier transform operator and
$\lambda=\textrm{0.05}$. The minimization in (13) was solved using the
nonlinear conjugate gradient algorithm [25]. Next, the remaining variables
$\mathbf{q}$, $\kappa$, and $\mathbf{f}$ were estimated by iterating over
lines 4 and 5 of Algorithm 1 while keeping $\mathbf{m}$ fixed. This approach
is referred to as the ”two-step” reconstruction.

##  IV Results

TABLE II: Reconstruction errors and estimated dynamical parameter. Experiment |  | Two-step |  | Spectro-Dynamic
---|---|---|---|---
Activation | $\alpha$ |  | RMSE u | RMSE v | RMSE f |  $\kappa$* |  | RMSE u | RMSE v | RMSE f |  $\kappa$*
| (°) |  | (mm) | (mm/s) | (N) | (N/m) |  | (mm) | (mm/s) | (N) | (N/m)
Continuous | 0 |  | 1.38 | 3.55 | 27$\cdot$10-3 | 53.8 |  | 0.24 | 1.00 | 9.5$\cdot$10-3 | 31.8
Continuous | 45 |  | 1.48 | 3.62 | 14$\cdot$10-3 | 58.4 |  | 0.11 | 0.87 | 7.0$\cdot$10-3 | 30.1
Continuous | 90 |  | 1.54 | 3.72 | 20$\cdot$10-3 | 67.0 |  | 0.20 | 0.96 | 8.6$\cdot$10-3 | 30.4
Discontinuous | 0 |  | 1.36 | 4.95 | 38$\cdot$10-3 | 28.4 |  | 0.23 | 1.42 | 17$\cdot$10-3 | 28.0
Discontinuous | 45 |  | 1.45 | 5.11 | 54$\cdot$10-3 | 18.3 |  | 0.11 | 0.85 | 15$\cdot$10-3 | 26.8
Discontinuous | 90 |  | 1.38 | 5.07 | 47$\cdot$10-3 | 22.3 |  | 0.19 | 1.36 | 16$\cdot$10-3 | 27.4
* Ground truth value: 30.0 N/m |  |  |  |  |  |

In the dynamic image series estimated in the two-step reconstruction, residual
motion blurring was visible (Fig. 6(a)). The images calculated from the
reconstructed spatio-temporal k-space data $\mathbf{m}$ did not show this
blurring (Fig. 6(b)). For validation purposes, we also reconstructed an image
from fully-sampled k-space data of a stationary phantom (Fig. 6(c)). Note that
Fig. 6(b) and (c) are very similar, although the phantom moved while acquiring
data for the Spectro-Dynamic MRI reconstruction (Fig. 6(b)). The reconstructed
image series had a temporal resolution of 11 ms (Supplementary Video S1).

The estimated displacements, velocities, and activation forces were close to
their respective ground truth values (Fig. 7 and Table II). Additionally, the
estimated dynamical parameter ($\kappa$) was very close to its ground truth
value of 30 N/m (Table II).

When rotating the readout direction with respect to the direction of motion,
the reconstruction accuracy remained high (Table II and Figs. S1 and S2 in the
Supplementary Material). In addition, when changing the activation force to a
discontinuous ’on/off’ function (Fig. 2, right column), the reconstructed
motion remained accurate (Table II and Fig. S3 in the Supplementary Material),
although the accuracy of the estimated activation force and dynamical
parameter was lower than for the reconstructions with the continuous
activation. Again, the accuracy remained similar when the readout direction
was rotated (Table II and Figs. S4 and S5 in the Supplementary Material).

##  V Discussion

Figure 7: Estimated (a) displacements, (b) velocities, and (c) activation
forces in both directions resulting from the two-step reconstruction (top) and
Spectro-Dynamic reconstruction (bottom) for the experiment with a continuous
activation function and $\alpha=0\degree{}$. The velocities were calculated
from the displacements using finite differences. The motion of both the moving
and the stationary compartment was estimated, but only the moving compartment
moved along one direction ($x$). The small insets show zoomed versions of the
graphs around the point where the motion starts.

In this work, we presented an extended reconstruction framework for Spectro-
Dynamic MRI that estimates time-resolved images, time-resolved motion fields,
dynamical parameters, and the activation force. We formulated an optimization
problem in terms of four components: a motion model that relates the
displacements to the fully-sampled k-space data, a dynamical model that
describes the expected type of motion, a data consistency term that describes
the sampling and measurement process, and a regularization term on the
activation force. We solved this optimization problem by using an iterative
reconstruction algorithm with convex subproblems. Phantom experiments with
different kinds of activations showed that the dynamics can be reconstructed
accurately at a temporal resolution of 11 ms. The reconstruction did not
assume any periodicity in the motion, and the accuracy was independent of the
orientation of the motion with respect to sampling pattern.

The strength of our method is the joint reconstruction of the missing k-space
data and the motion fields. We have shown the increased accuracy of our
proposed reconstruction method compared to a two-step approach, where time-
resolved images are reconstructed in the first step, after which the motion is
estimated from this data. Since no motion information is available during the
first step, a temporal smoothness regularization is used to resolve the high
undersampling, which results in motion blurring in the images. Estimating the
motion from this blurred data results in a large underestimation of the
displacements. In contrast, the iterative Spectro-Dynamic MRI reconstruction
leverages the motion information during the estimation of the time-resolved
data and vice versa. Thus, the measured information can be shared more
efficiently across time, and accurate reconstructions with high acceleration
factors per time instance (32 in this work) can be achieved.

The ordering of the subproblems in Algorithm 1 is important. In the first
iteration, the optimization for $\mathbf{m}$ (with
$\mathbf{q}^{0}=\mathbf{0}$) results in a smooth temporal interpolation of the
missing k-space data. This rough estimate for $\mathbf{m}$ is subsequently
used to generate a first estimate of the motion fields and activation force.
Finally, the dynamical parameters are estimated, while the activation force is
refined. Thus, the result after one iteration already resembles the eventual
solution. Subsequent iterations refine the solution, bringing it closer to the
ground truth values. Choosing a different ordering can result in convergence
problems, especially in the first iteration, and increases the sensitivity of
the solution to the chosen initial values.

The activation force $\mathbf{f}$ is optimized twice during each iteration of
Algorithm 1. Omitting the estimation of $\mathbf{f}$ in either line 4 or 5
results in a much slower convergence of $\mathbf{q}$ or $\kappa$,
respectively. Even though the temporary variable $\bar{\mathbf{f}}^{k}$ is not
used explicitly, it gives more freedom during the convex optimization of line
4, resulting in a better estimate for $\mathbf{q}$.

The damping coefficient $c$ was fixed during the optimization, since
estimating both $c$ and $\mathbf{f}$ would lead to complexities in the
reconstruction. For more complex and more realistic in vivo applications, some
mechanical parameters could be fixed to literature values. Alternatively, the
time-dependent activation function could be measured using e.g. ECG or EMG,
depending on the application.

Due to the limitations of the motion phantom we used, only one of the four
degrees of freedom of the motion was actually nonzero. Since the motion
estimation remained accurate when the readout direction was rotated, we expect
that motion in any of the four degrees of freedom can be reconstructed
accurately.

By using the data consistency term $H$ as a penalty term instead of an
equality constraint, the noise in the measurements can be captured by the
residual of this term. The idea is that this relaxation strategy reduces the
noise in $\mathbf{m}$ which is used in the motion model $G$. This in turn
results in a reduced bias of the estimated displacement field. This bias was
visible for example in the estimated displacements in the phase encode
direction in our previous work [18], which were less accurate than the
displacements in the readout direction. In our current results, no such
difference was present between the different directions. The addition of the
data consistency term in the reconstruction also removed the need for the low-
pass filter which was previously applied to the k-space data.

The motion model leverages the reconstructed displacement field to fill the
missing k-space data. Other techniques also aim to reconstruct images from
undersampled k-space data. Parallel imaging [26, 27] uses information about
the coil sensitivities to estimate the missing data, thereby removing the
undersampling artifacts. Compressed sensing [25] combines an incoherent
undersampling pattern with a sparse regularization to reconstruct an image. In
addition, multiple methods have been developed to accelerate dynamic MRI data
using temporal correlations in the data, without using the motion fields
explicitly [28]. These techniques could complement the motion model in
Spectro-Dynamic MRI, although this is not done in this work.

Our method requires several extensions before it can be used for in vivo
applications. Firstly, through-slice motion violates the assumption of signal
conservation used in the derivation of the motion model. Since pure in-plane
in vivo motion is rare, we will next investigate a 3D implementation for
Spectro-Dynamic MRI. Acquiring 3D data at a high temporal resolution will
require undersampling in both phase encode directions, resulting in an even
higher undersampling factor than in the current 2D implementation. We expect
that carefully choosing the sampling pattern will be important to keep a good
trade-off between spatial and temporal resolution. More than two readouts may
have to be included in a single time instance. Again, techniques such as
parallel imaging or compressed sensing may be used to deal with the increased
undersampling in 3D. Random undersampling patterns or non-Cartesian
trajectories might also increase the robustness of the reconstructions,
especially when using compressed sensing reconstruction approaches. Prior
information about the estimated images could also be added to the
reconstruction in the form of a smoothness or total variation regularizer.

Secondly, a partial differential equation (PDE) has to replace the ODE used as
dynamical model in this work. This PDE allows for the reconstruction of
deformable motion fields. For example, the Navier-Cauchy equation can be used
to model the behavior of isotropic linear elastic materials. In this case, the
non-rigid motion fields become a function of space as well as time. The
mechanical parameters might also become spatially dependent when different
tissue types or inhomogeneous tissues are considered. We will investigate
deformable motion estimation and the use of PDEs as dynamical model in the
Spectro-Dynamic MRI framework in the future.

A non-rigid motion field also calls for different basis functions
$\boldsymbol{\phi}(\mathbf{r})$. B-splines [29] could be used, since motion
fields are smooth almost everywhere. Furthermore, a low-rank description of
the motion field has been shown to be very effective [17]. Choosing adequate
basis functions to describe the motion field will be an important step in
extending Spectro-Dynamic MRI towards 3D in vivo applications.

Finally, determining constitutive models to describe the dynamic behavior of
tissues in vivo is not trivial, and many different (nonlinear) models can be
found in the literature [30]. Alternatively, data-driven discovery of
dynamical systems can be used to learn constitutive relations from the
dynamical data itself [31, 32]. Discovering dynamical models using Spectro-
Dynamic MRI is currently being investigated by the authors [33].

##  VI Conclusion

The extended Spectro-Dynamic MRI reconstruction framework presented in this
work allows for time-resolved dynamic MRI reconstructions at a high temporal
resolution. In addition, it gives insights into the underlying dynamics of the
motion through the dynamical model. The extended framework enables the
reconstruction of more general motion fields, regardless of the sampling
pattern adopted during the MRI scan.

## Acknowledgment

The authors would like to thank David Heesterbeek for fruitful discussions.

## References

  * [1] F. G. Shellock, “Functional Assessment of the Joints Using Kinematic Magnetic Resonance Imaging,” _Seminars in Musculoskeletal Radiology_ , vol. 7, no. 4, pp. 249–276, 2003.
  * [2] G. E. Gold, “Dynamic and Functional Musculoskeletal Imaging,” _Seminars in Musculoskeletal Radiology_ , vol. 7, no. 4, pp. 245–248, 2003\.
  * [3] L. M. Shapiro and G. E. Gold, “MRI of weight bearing and movement,” _Osteoarthritis and Cartilage_ , vol. 20, no. 2, pp. 69–78, 2012.
  * [4] B. Borotikar, M. Lempereur, M. Lelievre, V. Burdin, D. Ben Salem, and S. Brochard, “Dynamic MRI to quantify musculoskeletal motion: A systematic review of concurrent validity and reliability, and perspectives for evaluation of musculoskeletal disorders,” _PLoS ONE_ , vol. 12, no. 12, 2017\.
  * [5] J. Alcazar, R. Csapo, I. Ara, and L. M. Alegre, “On the Shape of the Force-Velocity Relationship in Skeletal Muscles: The Linear, the Hyperbolic, and the Double-Hyperbolic,” _Frontiers in Physiology_ , vol. 10, 2019.
  * [6] A. G. d’Entremont, J. A. Nordmeyer-Massner, C. Bos, D. R. Wilson, and K. P. Pruessmann, “Do dynamic-based MR knee kinematics methods produce the same results as static methods?” _Magnetic Resonance in Medicine_ , vol. 69, no. 6, pp. 1634–1644, 2013.
  * [7] L. E. Bilston and K. Tan, “Measurement of Passive Skeletal Muscle Mechanical Properties In Vivo: Recent Progress, Clinical Applications, and Remaining Challenges,” _Annals of Biomedical Engineering_ , vol. 43, no. 2, pp. 261–273, 2 2015.
  * [8] C. F. Guimarães, L. Gasperini, A. P. Marques, and R. L. Reis, “The stiffness of living tissues and its implications for tissue engineering,” _Nature Reviews Materials_ , vol. 5, pp. 351–370, 5 2020.
  * [9] F. T. Sheehan, F. E. Zajac, and J. E. Drace, “Using cine phase contrast magnetic resonance imaging to non-invasively study in vivo knee dynamics,” _Journal of Biomechanics_ , vol. 31, no. 1, pp. 21–26, 1998.
  * [10] J. Kaiser, R. Bradford, K. Johnson, O. Wieben, and D. G. Thelen, “Measurement of Tibiofemoral Kinematics Using Highly Accelerated 3D Radial Sampling,” _Magnetic Resonance in Medicine_ , vol. 69, no. 5, pp. 1310–1316, 2013.
  * [11] D. S. Asakawa, K. S. Nayak, S. S. Blemker, S. L. Delp, J. M. Pauly, D. G. Nishimura, and G. E. Gold, “Real-Time Imaging of Skeletal Muscle Velocity,” _Journal of Magnetic Resonance Imaging_ , vol. 18, no. 6, pp. 734–739, 12 2003.
  * [12] K. S. Nayak, Y. Lim, A. E. Campbell-Washburn, and J. Steeden, “Real-Time Magnetic Resonance Imaging,” _Journal of Magnetic Resonance Imaging_ , vol. 55, no. 1, pp. 81–99, 1 2022.
  * [13] V. Mazzoli, J. Schoormans, M. Froeling, A. M. Sprengers, B. F. Coolen, N. Verdonschot, G. J. Strijkers, and A. J. Nederveen, “Accelerated 4D self-gated MRI of tibiofemoral kinematics,” _NMR in Biomedicine_ , vol. 30, no. 11, 11 2017.
  * [14] C. B. Shaw, B. H. Foster, M. Borgese, R. D. Boutin, C. Bateni, P. Boonsri, C. O. Bayne, R. M. Szabo, K. S. Nayak, and A. J. Chaudhari, “Real-time three-dimensional MRI for the assessment of dynamic carpal instability,” _PLoS ONE_ , vol. 14, no. 9, 9 2019.
  * [15] R. G. Menon, M. V. Zibetti, and R. R. Regatte, “In Vivo Tibiofemoral Cartilage Strain Mapping Under Static Mechanical Loading Using Continuous GRASP-MRI,” _Journal of Magnetic Resonance Imaging_ , vol. 51, no. 2, pp. 426–434, 2 2020.
  * [16] L. Feng, Q. Wen, C. Huang, A. Tong, F. Liu, and H. Chandarana, “GRASP-Pro: imProving GRASP DCE-MRI through self-calibrating subspace-modeling and contrast phase automation,” _Magnetic Resonance in Medicine_ , vol. 83, no. 1, pp. 94–108, 1 2020.
  * [17] N. R. F. Huttinga, T. Bruijnen, C. A. T. van den Berg, and A. Sbrizzi, “Non-rigid 3D motion estimation at high temporal resolution from prospectively undersampled k-space data using low-rank MR-MOTUS,” _Magnetic Resonance in Medicine_ , pp. 1–18, 2020.
  * [18] M. H. C. Van Riel, N. R. F. Huttinga, and A. Sbrizzi, “Spectro-Dynamic MRI: Characterizing Mechanical Systems on a Millisecond Scale,” _IEEE Access_ , vol. 10, pp. 271–285, 2022.
  * [19] T. van Leeuwen and F. J. Herrmann, “A penalty method for PDE-constrained optimization in inverse problems,” _Inverse Problems_ , vol. 32, no. 1, pp. 1–26, 2016.
  * [20] B. K. Horn and B. G. Schunck, “Determining optical flow,” _Artificial Intelligence_ , vol. 17, no. 1-3, pp. 185–203, 1981.
  * [21] J. M. Winters and L. Stark, “Muscle Models: What Is Gained and What Is Lost by Varying Model Complexity,” _Biological Cybernetics_ , vol. 55, no. 6, pp. 403–420, 1987.
  * [22] Y. Xu and W. Yin, “A block Coordinate Descent Method for Regularized Multiconvex Optimization with Applications to Nonnegative Tensor Factorization and Completion,” _SIAM Journal on Imaging Sciences_ , vol. 6, no. 3, pp. 1758–1789, 2013.
  * [23] X. Shen, S. Diamond, M. Udell, Y. Gu, and S. Boyd, “Disciplined Multi-Convex Programming,” in _Proceedings of the 29th Chinese Control And Decision Conference_ , 2017, pp. 895–900.
  * [24] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,” _Foundations and Trends in Machine Learning_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [25] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” _Magnetic Resonance in Medicine_ , vol. 58, no. 6, pp. 1182–1195, 2007.
  * [26] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, “SENSE: Sensitivity encoding for fast MRI,” _Magnetic Resonance in Medicine_ , vol. 42, no. 5, pp. 952–962, 1999.
  * [27] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase, “Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA),” _Magnetic Resonance in Medicine_ , vol. 47, no. 6, pp. 1202–1210, 2002.
  * [28] J. Tsao and S. Kozerke, “MRI Temporal Acceleration Techniques,” _Journal of Magnetic Resonance Imaging_ , vol. 36, no. 3, pp. 543–560, 9 2012\.
  * [29] T. Rohlfing, C. R. Maurer, D. A. Bluemke, and M. A. Jacobs, “Volume-Preserving Nonrigid Registration of MR Breast Images Using Free-Form Deformation With an Incompressibility Constraint,” _IEEE Transactions on Medical Imaging_ , vol. 22, no. 6, pp. 730–741, 6 2003.
  * [30] G. Chagnon, M. Rebouah, and D. Favier, “Hyperelastic Energy Densities for Soft Biological Tissues: A Review,” _Journal of Elasticity_ , vol. 120, no. 2, pp. 129–160, 12 2015.
  * [31] S. L. Brunton, J. L. Proctor, J. N. Kutz, and W. Bialek, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [32] M. Flaschel, S. Kumar, and L. De Lorenzis, “Unsupervised discovery of interpretable hyperelastic constitutive laws,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 381, 8 2021.
  * [33] D. G. J. Heesterbeek, M. H. C. van Riel, C. A. T. van den Berg, and A. Sbrizzi, “Data-driven Discovery of mechanical models directly from k-space data: Initial results with an actively forced linear elastic system,” in _Proceedings of the 31st Annual Meeting of ISMRM_ , no. 2552, 2023.

Max H. C. van Riel received the B.Sc. degree in Biomedical Engineering from
Eindhoven University of Technology, Eindhoven, The Netherlands, in 2018, and
the M.Sc. degree in Medical Imaging from Eindhoven University of Technology,
Eindhoven, The Netherlands in 2020. In 2019, he was a Student Intern at the
New York University School of Medicine for six months. He is currently
pursuing a Ph.D. degree at the UMC Utrecht, Utrecht, The Netherlands. In his
research, he focuses on characterizing dynamical systems at a high temporal
resolution using MRI.
---
Tristan van Leeuwen received the M.Sc. degree in Computational Science from
Utrecht University, Utrecht, The Netherlands, in 2006, and the Ph.D. degree in
Geophysics from Delft University of Technology, Delft, The Netherlands, in
2010. He was a Postdoctoral Fellow at the University of British Columbia from
2010-2013, and at the Centrum voor Wiskunde en Informatica (CWI) from
2013-2014. He was an Assistant Professor at the Mathematical Institute of
Utrecht University, Utrecht, The Netherlands. Currently, he leads the
Computational Imaging Group at CWI in Amsterdam, The Netherlands. His research
interests include: inverse problems, computational imaging, tomography and
numerical optimisation.
---
Cornelis A. T. van den Berg received the M.Sc. degree in Applied Physics from
University of Twente, Enschede, The Netherlands, in 2000, and the Ph.D. degree
from Utrecht University, Utrecht, The Netherlands, in 2006. At present, he is
a Full Professor at the UMC Utrecht, Utrecht, The Netherlands, where he leads
the Computational Imaging Group. Here, he is (co)leading various research
projects on MRI methods & hardware for MRI guided radiotherapy, new MRI
diagnostics using MR-STAT technology, automated imaging workflows for
radiotherapy using deep learning and his first love on electromagnetics of
MRI.
---
Alessandro Sbrizzi received the B.Sc. and M.Sc. degrees in Mathematics in 2009
and the Ph.D. degree in 2013, all from Utrecht University, Utrecht, The
Netherlands. He is currently an Associate Professor at the Computational
Imaging group of the UMC Utrecht, Utrecht, The Netherlands. His research
focuses on fast multi-parametric MRI (in particular MR-STAT), real-time
motion-estimation (MR-MOTUS technique), radiofrequency pulse design and the
application of machine learning in MRI.
---

Time-Resolved Reconstruction of Motion, Force, and Stiffness using Spectro-
Dynamic MRI

Supplementary material

Max H. C. van Riel1, Tristan van Leeuwen2, Cornelis A. T. van den Berg1, and
Alessandro Sbrizzi1

1Department of Radiotherapy, Computational Imaging Group for MR Diagnostics
and Therapy, UMC Utrecht, Utrecht, The Netherlands

2Mathematical Institute, Utrecht University, Utrecht, The Netherlands

Figure S1: Estimated (a) displacements, (b) velocities, and (c) activation
forces in both directions resulting from the two-step reconstruction (top) and
Spectro-Dynamic reconstruction (bottom) for the experiment with a continuous
activation function and $\alpha=45\degree{}$. The velocities were calculated
from the displacements using finite differences. The motion of both the moving
and the stationary compartment was estimated, but only the moving compartment
moved along one direction. Since $\alpha=45\degree{}$, this motion was
separated over the $x$\- and $y$-directions.

Figure S2: Estimated (a) displacements, (b) velocities, and (c) activation
forces in both directions resulting from the two-step reconstruction (top) and
Spectro-Dynamic reconstruction (bottom) for the experiment with a continuous
activation function and $\alpha=90\degree{}$. The velocities were calculated
from the displacements using finite differences. The motion of both the moving
and the stationary compartment was estimated, but only the moving compartment
moved along one direction ($y$).

Figure S3: Estimated (a) displacements, (b) velocities, and (c) activation
forces in both directions resulting from the two-step reconstruction (top) and
Spectro-Dynamic reconstruction (bottom) for the experiment with a
discontinuous activation function and $\alpha=0\degree{}$. The velocities were
calculated from the displacements using finite differences. The motion of both
the moving and the stationary compartment was estimated, but only the moving
compartment moved along one direction ($x$).

Figure S4: Estimated (a) displacements, (b) velocities, and (c) activation
forces in both directions resulting from the two-step reconstruction (top) and
Spectro-Dynamic reconstruction (bottom) for the experiment with a
discontinuous activation function and $\alpha=45\degree{}$. The velocities
were calculated from the displacements using finite differences. The motion of
both the moving and the stationary compartment was estimated, but only the
moving compartment moved along one direction. Since $\alpha=45\degree{}$, this
motion was separated over the $x$\- and $y$-directions.

Figure S5: Estimated (a) displacements, (b) velocities, and (c) activation
forces in both directions resulting from the two-step reconstruction (top) and
Spectro-Dynamic reconstruction (bottom) for the experiment with a
discontinuous activation function and $\alpha=90\degree{}$. The velocities
were calculated from the displacements using finite differences. The motion of
both the moving and the stationary compartment was estimated, but only the
moving compartment moved along one direction ($y$).

Video S1: Video showing the time-resolved image data. Each frame is obtained
by performing a Fast Fourier Transform (FFT) on one time instance of the
estimated time-resolved k-space data $\mathbf{m}$, followed by a geometry
correction to account for gradient nonlinearities. The video can be viewed
online, this document only shows a single frame.


# Position Bias Estimation with Item Embedding for Sparse Dataset

Shion Ishikawa<EMAIL_ADDRESS>Rakuten Group, Inc.TokyoJapan ,
Yun Ching Liu<EMAIL_ADDRESS>Rakuten Group, Inc.TokyoJapan ,
Young-joo Chung<EMAIL_ADDRESS>Rakuten Group, Inc.CaliforniaUSA
and  Yu Hirate<EMAIL_ADDRESS>Rakuten Group, Inc.TokyoJapan

(2024)

###### Abstract.

Estimating position bias is a well-known challenge in Learning to Rank (L2R).
Click data in e-commerce applications, such as targeted advertisements and
search engines, provides implicit but abundant feedback to improve
personalized rankings. However, click data inherently includes various biases
like position bias. Based on the position-based click model, Result
Randomization and Regression Expectation-Maximization algorithm (REM) have
been proposed to estimate position bias, but they require various paired
observations of (item, position). In real-world scenarios of advertising,
marketers frequently display advertisements in a fixed pre-determined order,
which creates difficulties in estimation due to the limited availability of
various pairs in the training data, resulting in a sparse dataset. We propose
a variant of the REM that utilizes item embeddings to alleviate the sparsity
of (item, position). Using a public dataset and internal carousel
advertisement click dataset, we empirically show that item embedding with
Latent Semantic Indexing (LSI) and Variational Auto-Encoder (VAE) improves the
accuracy of position bias estimation and the estimated position bias enhances
Learning to Rank performance. We also show that LSI is more effective as an
embedding creation method for position bias estimation.

EM algorithm; Position bias estimation; Item embedding

††journalyear: 2024††copyright: acmlicensed††conference: Companion Proceedings
of the ACM Web Conference 2024; May 13–17, 2024; Singapore,
Singapore.††booktitle: Companion Proceedings of the ACM Web Conference 2024
(WWW ’24 Companion), May 13–17, 2024, Singapore, Singapore††isbn:
979-8-4007-0172-6/24/05††doi: 10.1145/3589335.3651546††ccs: Information
systems Display advertising

##  1\. Introduction

Learning to Rank (L2R) (Cao et al., 2007) is a valuable approach to improving
ranking performance in search engines, recommender systems (Yuan et al., 2017)
(Rendle et al., 2009), and targeted advertisement (Aharon et al., 2019). Click
data from rankings applications is a promising source for enhancing
personalized rankings. Unlike explicit feedback, such as user ratings, click
data reflects users’ interaction with applications, making it a valuable form
of _implicit feedback_. However, due to its implicit nature, click data
contains various biases, including position bias (Joachims et al., 2017).
Position bias implies that the location of content influences its likelihood
of being selected by users. To extract reliable signals from users’ implicit
feedback, previous research has focused on removing bias, in other words,
denoising click data.

Click-Modeling is one approach to addressing the aforementioned issues. For
instance, the position-based click model (PBM) (Richardson et al., 2007)
represents the probability of a user clicking on an item as the product of two
probabilities: the probability of a user examining a specific position
(position bias), and the probability of a user clicking on an item after
examining it (item-user relevance).

In this study, we focus on estimating position bias under the PBM approach.
While the cascade model (Craswell et al., 2008) has been extensively studied,
it assumes a linear order of examining items and may not be applicable in all
industrial scenarios where the order of item examination is not linear or
follows a specific pattern, such as in the case of carousel ads with
navigation buttons. For example, in the case of carousel ads with navigation
buttons, when users click the left arrow button, the displayed ads shift to
the left, causing users to examine the ads in a right-to-left order.
Conversely, if they click the right arrow button, the order of examination
will be reversed, with users examining the ads from left to right.

The accuracy of position bias estimation is critical for the performance of
recommender systems, because an inaccurate estimation will directly lead to a
biased evaluation of item-user relevance. We have implemented the regression
EM algorithm (Wang et al., 2018) in the production carousel advertisement to
estimate position bias from logged click data, without affecting real user
traffic. Unfortunately, we discovered that this model didn’t work well for
this application due to _the training data sparsity in the position bias
estimation_.

The regression EM algorithm takes the tuple (item, position) as input data.
Ideally, the items should be evenly distributed across all positions with a
sufficient number of samples. _However, in real-world scenarios, marketers
often display advertisements in a fixed pre-determined order, resulting in
most items being placed in only a few positions._ Moreover, in many cases, we
observe one item is consistently placed in a single position. When we have 15
ads in our e-commerce platform, It means that we observe only
$6.7\%(=15/15^{2})$ of total combinations in (ad, position) tuples. The
position bias estimation suffers from it.

In this study, we propose a method to address this sparsity issue. Our
approach utilizes item similarities to estimate missing values in the (item,
position) matrix. Specifically, when the (item $X$, position 1) entry is
missing, we estimate its position using a similar item $X^{\prime}$. To
achieve this, we employ embedding techniques to create a representation that
preserves similarity between items. In our investigation, we discovered that
item embeddings for estimating position bias have not been tested extensively
with well-known tabular embedding baseline models like Latent Semantic
Indexing (LSI) (Rosario, 2000) and Variational Autoencoder (VAE) (Kingma et
al., 2019). Therefore, we aim to leverage these two models as an initial
proposal in this study.

In summary, our study makes the following key contributions:

  * •

We introduce the novel problem of data sparsity in (item, position) pairs for
position bias estimation in a real-world scenario.

  * •

We propose a novel method111We plan to publish code on
https://github.com/rakutentech for estimating position bias using Latent
Semantic Indexing and Variational auto-encoder.

  * •

We conduct extensive experiments on public and internal datasets, which
demonstrates that our method, REM with item embedding, estimates position
biases more accurately and improves learning-to-rank performance.

  * •

We find that LSI outperforms VAE in terms of accuracy of position bias
estimation and personalized ranking.

##  2\. Problem Setting

###  2.1. Position-based click model

Let $i\in\mathcal{I}$ be an item to be recommended (e.g., an advertisement or
a product). Let $C\in\mathbb{R}_{\geqslant 0}$ denote a reward variable (e.g.,
whether the displayed item was clicked or not). Let $\mathcal{U}$ denote a set
of contexts, such as a user’s demographic features (e.g., age and occupation).
Let $k$ be a position at which an item is displayed. The position-based click
model assumes that the probability of clicks $P(C=1|i,u,k)$ can be represented
as a product of two latent probabilities, as follows:

(1) $P(C=1|i,u,k)=P(R=1|i,u)P(E=1|k)$

where $P(E=1|k)$ represents the probability that position $k$ is examined by a
user and $P(R=1|i,u)$ represents the probability that the item is related to
the user. In short, the position-based click model assumes that if a user
examines a particular position and the item in the position is relevant, they
will click on the item. This model also assumes the examination process solely
depends on a position and does not depend on other items. For simplicity, we
denote relevance as $\mu(i,u)=P(R=1|i,u)$ and position bias as
$\theta_{k}=P(E=1|k)$.

###  2.2. Data sparsity in the position-based click model

In the context of the position-based click model presented in Eq 1, let
$\mathcal{D}:=\\{(u_{j},i_{j},r_{j},k_{j})\\}^{n}_{j=1}$ be a logged dataset
of $n$ observations. Let $\mathcal{A}$ = $\\{(i,k)\\}$ be a set of actions
containing all possible pairs of $(i,k)$. We define a function
$\pi:\mathcal{U}\to\Delta(\mathcal{A})$ as a policy, which maps each context
$u\in\mathcal{U}$ to a distribution over the actions, where $\pi(i,k)$
represents the probability of assigning item $i$ to position $k$.
$\Delta(\mathcal{A})$ is a simplex which satisfies
$\Delta(\mathcal{A})=\\{a\in\mathbb{R}^{K}|\sum_{j=0}^{K}a_{j}=1,a_{j}\geq
0\\}$ where $K$ is a number of actions. The policy that creates the logging
dataset $D$ is referred to as the logging policy.

In most cases, the logging policy is often determined by marketing experts
based on their experience and expertise. As a result, it is mostly
deterministic and static. As a result, the logged dataset contains a limited
variation of $(i,k)$ pairs, which can lead to data sparsity that affects the
estimation precision of position bias.

##  3\. Method

###  3.1. Item embeddings and modified representation of position-based click
model

As we discussed in the above sections, the sparsity problem exists for the
$(i,k)$ pairs in the dataset $\mathcal{D}$. Hence, the estimation of position
bias could be inaccurate due to spare and skewed data distribution. To
mitigate this issue, we represent an item with latent embedding vectors
$e\in\mathcal{E}$. Formally, we create a function that maps the $n\times l$
(item, feature) matrix into (item, embedding vector) $n\times m$ matrix where
$m<<l$. With the obtained embeddings, we then transform a sparse tuple $(i,k)$
into a dense tuple $(e,k)$. We first introduce the probability of an embedding
vector given an item $P(e|i)$. Then we represent the probability of item-user
relevance $P(R=1|i,u)$ with embedding vectors, as follows.

(2) $P(R=1|i,u)=\sum_{e\in\mathcal{E}}P(e|i)P(R=1|e,u)$

where $\sum_{e\in\mathcal{E}}P(e|i)=1$. Based on $P(e|i)$ and Eq 2, we sample
partial reward $w$ from click $C$ based on the following equation.

(3) $P(w=1|e,u,k)=P(e|i)P(C=1|i,u,k)$

Then we can rewrite dataset
$\mathcal{D}=\\{(u_{j},i_{j},c_{j},k_{j})\\}^{n}_{j=1}$ as
$\mathcal{D}_{e}:=\\{(u_{j},e_{j},w_{j},k_{j})\\}^{n|\mathcal{E}|}_{j=1}$.
Based on $\mathcal{D}_{e}$, we can rewrite the position-based click model in
Eq. 1 with embedding vectors.

(4) $P(w=1|e,u,k)=P(R=1|e,u)P(E=1|k)$

Table 1. A toy example illustrating the benefit of item embedding against a
sparse dataset. Left table: logging policy (the distribution of item and their
position assignment) which induces sparsity. It can be seen that item $i_{0}$
is only shown in position $k_{0}$. Right table: Assignment of embedding vector
to each item. Bottom table: the resampled dataset with a new logging policy
$\pi(e_{l},k_{l})=\sum_{j}p(e_{l}|i_{j})\pi(i_{j},k_{l})$, showcasing reduced
sparsity through item embeddings.

| $\pi_{b}(i,k_{0})$ | $\pi_{b}(i,k_{1})$ | $\pi_{b}(i,k_{2})$
---|---|---|---
$i_{0}$ | 1/3 | 0 | 0
$i_{1}$ | 0 | 1/3 | 0
$i_{2}$ | 0 | 0 | 1/3

| $p(e_{0}|i)$ | $p(e_{1}|i)$
---|---|---
$i_{0}$ | 1/2 | 1/2
$i_{1}$ | 1/3 | 2/3
$i_{2}$ | 1/4 | 3/4

| $\pi(e,k_{0})$ | $\pi(e,k_{1})$ | $\pi(e,k_{2})$
---|---|---|---
$e_{0}$ | 1/6 | 1/9 | 1/12
$e_{1}$ | 1/6 | 2/9 | 3/12

Table 1 shows a toy example in which the sparsity of the original dataset is
mitigated through item embedding.

###  3.2. Calculating the probability of embedding vector given item

To use embeddings for position-based click models, we need to convert
$\mathcal{E}$ into the probability space $[0,1]$ to represent $p(e|i)$. We
utilize the soft-max function to represent the conditional probability of
embeddings with a given item, as introduced by Saito et al. (Saito and
Joachims, 2022). However, it is important to note that while Saito’s work
focused on building an off-policy evaluation with action embeddings for a
large number of actions, our work addresses the issue of accurately estimating
position bias in a sparse dataset. We define $p(e|i)$ as shown below:

(5)
$p(e|i)=\prod_{j=1}^{m}\frac{\exp(e_{i,j})}{\sum_{j^{\prime}=1}^{m}\exp(e_{i,j^{\prime}})}$

where $e_{i,j}\in\mathcal{E}$ refers to an element obtained from the embedding
matrix (item $\times$ embedding). Here, we assume that when two items (e.g.,
item $i$ and $i^{\prime}$) are similar, their corresponding representations in
embedding space ($e_{i,j}$ and $e_{i^{\prime},j}$) will be similar, and which
leads to the similar conditional probabilities of embedding with a given item
($p(e|i)$ and $p(e|i^{\prime})$).

###  3.3. Regression EM algorithm with embedding vectors

Our goal is to obtain an accurate position bias $\theta_{k}$. Under the
position-based click model in Eq. 1, we apply the Regression EM algorithm
(Wang et al., 2018) to optimize $P(R=1|i,u)$ and $P(E=1|k)$ from the observed
$P(C=1|i,u,k)$.

We summarize our Regression EM algorithm with embedding in Algorithm 1. At
iteration $t+1$, the Expectation step will update the distribution of hidden
variable $E$ and $R$ from $\theta_{k}^{(t)}$ and $\mu^{(t)}(i,u)$. We use the
Gradient Boosted Decision Tree (GBDT) method to learn the function $\mu(u,e)$.

$\displaystyle P(E=1,R=1|w=1,e,k)$ $\displaystyle=1$ $\displaystyle
P(E=1,R=0|w=0,e,k)$
$\displaystyle=\frac{\theta_{k}^{(t)}(1-\mu^{(t)}(e,u))}{1-\theta_{k}^{(t)}\mu^{(t)}(e,u)}$
$\displaystyle P(E=0,R=1|w=0,e,k)$
$\displaystyle=\frac{(1-\theta_{k}^{(t)})\mu^{(t)}(e,u)}{1-\theta_{k}^{(t)}\mu^{(t)}(e,u)}$
(6) $\displaystyle P(E=0,R=0|w=0,e,k)$
$\displaystyle=\frac{(1-\theta_{k}^{(t)})(1-\mu^{(t)}(e,u))}{1-\theta_{k}^{(t)}\mu^{(t)}(e,u)}$

The Maximization step calculates $\theta_{k}^{(t+1)}$ and $\mu^{(t+1)}(e,u)$
using the probabilities from the Expectation step.

$\displaystyle\theta_{k}^{(t+1)}=\frac{\sum_{(u,e,w,k^{\prime})\in\mathcal{D}_{e}}\mathcal{I}_{k^{\prime}=k}(w+(1-w)P(E=1|u,e,w,k))}{\sum_{(u,e,w,k^{\prime})\in\mathcal{D}_{e}}\mathcal{I}_{k^{\prime}=k}}$
(7)
$\displaystyle\mu^{(t+1)}(e,u)=\frac{\sum_{(u,e^{\prime},w,k)\in\mathcal{D}_{e}}\mathcal{I}_{e^{\prime}=e}(w+(1-w)P(R=1|u,e,w,k))}{\sum_{(u,e^{\prime},w,k)\in\mathcal{D}_{e}}\mathcal{I}_{e^{\prime}=e}}$

where $\mathcal{I}$ shows the indicator function.

Algorithm 1 Regresssion EM algorithm with embedding

logging dataset $\mathcal{D}=\\{($user’s context $u,$ item $i,$ click $c,$
position $k)\\}$, position bias: $\\{\theta_{k}\\}$, relevance between
embedding $e$ and user’s context $u$: $\mu(u,e)$, assignment of embedding to
item: $P(e|i)$

for all $(u,i,c,k)\in\mathcal{D}$ do

Sample $w\in\\{0,1\\}$ from $c$ with $P(e|i)$ based on Eq. 3.

end for

Prepare $\mathcal{D}_{e}=\\{(u,e,w,k)\\}$ from sampled $\mathbf{w}$

repeat

Let $S=\\{\\}$

for all $(u,e,w,k)\in\mathcal{D}_{e}$ do

Sample $r\in\\{0,1\\}$ from $P(R=1|u,e,w,k)$ based on Eq. 3.3.

$S=S\cup(u,e,r)$

end for

Update $\mu(u,e)=GBDT(\mu(u,e),S)$

Update $\\{\theta_{k}\\}$ based on Eq. 3.3.

until Convergence

return $\\{\theta_{k}\\}$, $\mu(e,u)$

##  4\. Experiments

Table 3. Impact on ranking metrics with improved position bias | Open Bandit Dataset | Rakuten Ichiba dataset
---|---|---
| MRR | MAP@5 | MAP@10 | MRR | MAP@5 | MAP@10
Regression EM (REM) | 8.773e-3 | 1.709e-3 | 8.773e-4 | 0.7179 | 0.1573 | 7.394e-2
VAE + REM | 8.777e-3 | 1.712e-3 | 8.777e-4 | 0.7599 | 0.1583 | 7.673e-2
LSI + REM | 9.049-3 | 1.708e-3 | 9.049-4 | 0.7716 | 0.1613 | 7.795e-2

###  4.1. Experiments setup

In this section, we conduct extensive experiments to answer the following
questions.

  * •

(RQ1): How accurately does the proposed method estimate the position bias from
a sparse dataset?

  * •

(RQ2): How does the proposed method perform in a real-world application?

####  4.1.1. Dataset

We evaluate the proposed method with following two datasets:

  * •

Open Bandit Dataset (Saito et al., 2020): It is a large public dataset, which
contains 1,374,327 instances with 88 features and 80 items. The dataset has
two logging policies: a uniform random and Thompson Sampling. It has only
three positions which users see simultaneously and the difference among
trained position biases was less than 0.1%. Therefore, we instead introduce 10
virtual positions and resample clicks. We follow the procedure in (Ren et al.,
2022) to sample clicks. In the PBM model (Eq 1), we replace $P(E=1|k)$ with
position biases from eBay’s experiment (Aslanyan and Porwal, 2019), and we
resample clicks as $click_{i,u,k}\thicksim U(P(C=1|i,u,k))$.

  * •

Rakuten Ichiba Dataset: Rakuten is one of the largest e-commerce companies in
Japan. This platform has a widget of carousel ads containing 15 ads, with
three ads displayed simultaneously. From the four-day access log of this
widget, we obtained 494,329 instances with 156 features. As an A/B test, one
logging policy was determined heuristically by marketers, leading to sparse
pairs of (item, position). The other policy was a uniform random as control
group.

###  4.2. Effect on position bias estimation on sparse datasets (RQ1)

To evaluate the position bias estimation against a sparse datasets, we
constructed a sparse dataset by restricting an item to being assigned to a
single position from a dataset with a uniform random logging policy. We
applied the proposed method to this sparse dataset and assessed its ability to
recall the position bias of the original dataset in terms of RMSE. But in Open
Bandit Dataset, because we know the true values of the virtual position
biases, we compare our trained values with them.

Figure 1. The comparison of position bias estimation in Open Bandit Dataset
(top) and Rakuten Ichiba Dataset (bottom).

Lines show position biases estimated from different conditions. (Blue):
original dataset with Regression EM. (Yellow): sparse dataset with Regression
EM. (Green): sparse dataset with Latent Semantic Indexing. (Red): sparse
dataset with Variational Auto-Encoder. (Violet): Virtual position biases from
the experiment (Aslanyan and Porwal, 2019). We regard them as true position
biases.

Table 2. RMSE of position bias estimation for a sparse dataset | Open Bandit Dataset | Rakuten Ichiba dataset
---|---|---
REM | 0.271 $\pm$ 0.0277 | 0.292 $\pm$ 0.1138
VAE + REM | 0.257 $\pm$ 0.1147 | 0.265 $\pm$ 0.1203
LSI + REM | 0.219 $\pm$ 0.0989 | 0.260 $\pm$ 0.1137

Fig.1 and Table 2 show the performance when we applied the regression EM
algorithm (REM) with Latent Semantic Index, REM with Variational Auto-Encoder,
and vanilla REM to the sparse dataset. Table 2 shows the mean and standard
deviation of RMSE with 5 trials. It can be seen that VAE improved RMSE
relatively by 5.1% and LSI improved RMSE relatively by 19.2%. This result
indicates the item embedding leads to the accurate estimation of position
bias, and LSI outperforms VAE in terms of RMSE.

###  4.3. Effects on Ranking recommendation (RQ2)

We evaluate the performance of learned position bias estimation in the
recommendation model. We first fixed position biases with the learned values
from the previous experiment, and then we trained GDBT estimators under the
PBM model. Table 3 shows the performance of mean reciprocal rank (MMR) and
mean average precision (MAP) with the top 5 positions and top 10 positions.
Values are the mean of 5 trials. This result indicates that item embedding
improved ranking metrics, and except for MAP@5 in the Open Bandit Dataset, LSI
performs better than VAE.

##  5\. Conclusion

In this paper, we first introduced the data sparsity issue in the position
bias estimation under the Position-based click model. To overcome this issue,
we proposed the Regression EM algorithm with embedding (with LSI and VAE). By
utilizing the assignment probability $p(e|i)$ obtained from the item
embedding, we sampled various pairs of (embedding vector, position) and
applied the Regression EM algorithm to the sampled dataset. With the two real-
world datasets, we validated that our method outperforms the vanilla REM in
terms of both position bias estimation and recommendation.

## References

  * (1)
  * Aharon et al. (2019) Michal Aharon, Oren Somekh, Avi Shahar, Assaf Singer, Baruch Trayvas, Hadas Vogel, and Dobri Dobrev. 2019.  Carousel ads optimization in yahoo gemini native. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 1993–2001.
  * Aslanyan and Porwal (2019) Grigor Aslanyan and Utkarsh Porwal. 2019.  _Position Bias Estimation for Unbiased Learning-to-Rank in eCommerce Search_.  Springer International Publishing, 47–64.  https://doi.org/10.1007/978-3-030-32686-9_4
  * Cao et al. (2007) Zhe Cao, Tao Qin, Tie-Yan Liu, Ming-Feng Tsai, and Hang Li. 2007.  Learning to rank: from pairwise approach to listwise approach. In _Proceedings of the 24th international conference on Machine learning_. 129–136.
  * Craswell et al. (2008) Nick Craswell, Onno Zoeter, Michael Taylor, and Bill Ramsey. 2008.  An experimental comparison of click position-bias models. In _Proceedings of the 2008 international conference on web search and data mining_. 87–94.
  * Joachims et al. (2017) Thorsten Joachims, Laura Granka, Bing Pan, Helene Hembrooke, and Geri Gay. 2017.  Accurately interpreting clickthrough data as implicit feedback. In _Acm Sigir Forum_ , Vol. 51. Acm New York, NY, USA, 4–11.
  * Kingma et al. (2019) Diederik P Kingma, Max Welling, et al. 2019\.  An introduction to variational autoencoders.  _Foundations and Trends® in Machine Learning_ 12, 4 (2019), 307–392.
  * Ren et al. (2022) Yi Ren, Hongyan Tang, and Siwen Zhu. 2022.  Unbiased Learning to Rank with Biased Continuous Feedback. In _Proceedings of the 31st ACM International Conference on Information & Knowledge Management_ _(CIKM ’22)_. ACM.  https://doi.org/10.1145/3511808.3557483
  * Rendle et al. (2009) Steffen Rendle, Leandro Balby Marinho, Alexandros Nanopoulos, and Lars Schmidt-Thieme. 2009.  Learning optimal ranking with tensor factorization for tag recommendation. In _Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining_. 727–736.
  * Richardson et al. (2007) Matthew Richardson, Ewa Dominowska, and Robert Ragno. 2007\.  Predicting clicks: estimating the click-through rate for new ads. In _Proceedings of the 16th international conference on World Wide Web_. 521–530.
  * Rosario (2000) Barbara Rosario. 2000\.  Latent semantic indexing: An overview.  _Techn. rep. INFOSYS_ 240 (2000), 1–16.
  * Saito and Joachims (2022) Yuta Saito and Thorsten Joachims. 2022.  Off-Policy Evaluation for Large Action Spaces via Embeddings.  _arXiv preprint arXiv:2202.06317_ (2022).
  * Saito et al. (2020) Yuta Saito, Aihara Shunsuke, Matsutani Megumi, and Narita Yusuke. 2020.  Open Bandit Dataset and Pipeline: Towards Realistic and Reproducible Off-Policy Evaluation.  _arXiv preprint arXiv:2008.07146_ (2020).
  * Wang et al. (2018) Xuanhui Wang, Nadav Golbandi, Michael Bendersky, Donald Metzler, and Marc Najork. 2018\.  Position bias estimation for unbiased learning to rank in personal search. In _Proceedings of the Eleventh ACM International Conference on Web Search and Data Mining_. 610–618.
  * Yuan et al. (2017) Fajie Yuan, Guibing Guo, Joemon M Jose, Long Chen, Haitao Yu, and Weinan Zhang. 2017\.  Boostfm: Boosted factorization machines for top-n feature-based recommendation. In _Proceedings of the 22nd International Conference on Intelligent User Interfaces_. 45–54.