Weierstrass elliptic and related functions have been recently shown to enable analytical explicit solutions to classical problems in astrodynamics. These include the constant radial acceleration problem, the Stark problem and the two-fixed center (or Euler's) problem. In this paper we review the basic technique that allows for these results and we discuss the limits and merits of the approach. Applications to interplanetary trajectory design are then discussed including low-thrust planetary fly-bys and the motion of an artificial satellite under the influence of an oblate primary including $J_2$ and $J_3$ harmonics.

§ NOMENCLATURE

$f$  Generic function, polynomial

$P$  Generic 3rd or 4th order polynomial

$z, w$  Complex variables

$\Im$  Imaginary part of a complex number

$\Re$  Real part of a complex number

$\wp$  Weierstrass elliptic function

$\wp^{-1}$  Inverse of the Weierstrass elliptic function

$\zeta$  Weierstrass Zeta function

$\sigma$  Weierstrass Sigma function

$g_2, g_3$  Lattice invariants

$\omega_1, \omega_2, \omega_3$  Half-periods of the Weierstrass elliptic and related functions

$e_1, e_2, e_3$  Lattice roots

$\tilde e_1, \tilde e_2, \tilde e_3 (, \tilde e_4)$  Roots of the third (or fourth) order polynomial

$\tau$  Sundmann transformed pseudo-time (or anomaly)

$\tau_m$  Pseudo-time of root passage

$A, B$  Constants

$\mathcal E$  Specific energy

$\mathcal T$  Kinetic energy

$\mathcal V$  Potential

$h$  angular momentum

$a$  semi-major axis

$e$  eccentricity

$p$  orbital parameter

$\gamma$  flight path angle

$\delta$  asymptote deflection angle

$\mu$  gravitational parameter

$\alpha$  spacecraft acceleration

$r, v$  radius and velocity magnitudes

$\xi, \eta$  parabolic coordinates

$\theta$  true anomaly

Special notation

$\mathbf r, \mathbf v$  vectors

$f'$  derivative

$\dot f$  time derivative

$m$  referred to a polynomial root (root passage, pericenter)

$K$  Keplerian

$\infty$  At infinity

§ INTRODUCTION

Weierstrass elliptic and related functions appear in the solution of many problems in physics. In general relativity, for example, they are an established tool to tackle complex issues [1, 2, 3, 4]. In astrodynamics, only recently, they have been used to find explicit solutions to three fundamental problems: the constant radial acceleration problem[5], the Stark problem[6] and the two fixed center problem[7], also known as Euler's three body problem. The constant radial acceleration problem consists in describing the motion of a point mass particle subject to a central gravity field and to an additional constant radial acceleration. The Stark problem consists in describing the motion of a point mass particle subject to a central gravity field and to an additional acceleration constant in the inertial reference frame. The Euler's three body problem consists in describing the motion of a point mass particle subject to the gravity field of two masses fixed in the inertial frame. In all cases, the resulting dynamical system is integrable in the Liouville sense and in all cases the resulting dynamics, extensively studied both from a theoretical and an applicative perspective, admits an explicit analytical solution via the use of Weierstrass functions and the introduction of anomalies (or pseudo-times) introduced via Sundmann transformations. The actual solution in the real time is recovered by solving a transcendental function of such anomalies or pseudo-times (analogues to Kepler's equation). It is interesting to both note and further study the close analogy to the solution to Kepler's problem. Indeed the procedures and expressions are, at least formally, analogues of the Keplerian ones: the main difference being in their use of Weierstrass elliptic and related functions rather than of circular functions.

In this paper we review the generic solution procedure that allowed to obtain these results and we discuss the computer implementation of the new resulting procedures. We start with a basic introduction to Weierstrass elliptic and related functions and their relation to applications in astrodynamics. We then discuss the computer implementation of these functions showing how their evaluation cost is, essentially, comparable to that of circular functions when we restrict the evaluation to the real axis and assume the lattice properties as known. In the following section the constant radial acceleration problem is considered and a procedure to solve the initial value problem is detailed where evaluations of the Weierstrass functions are kept in the real domain in most of the cases. Then, the case of a radially powered planetary fly-by is studied in detail and an expression returning the asymptote deflection angle is developed. In the following section we turn our attention to the 2D Stark problem, first deriving the full solution, and then studying a second case of powered fly-by, here called the Stark fly-by. We find analytical expression that allow to design such a fly-by with ease thus allowing to study the exploitation of Oberth effect in a low-thrust trajectory. In the final section, we briefly show how the problem of artificial satellite motion under a gravity field including the $J_2$ and $J_3$ perturbations also has an analytical, explicit solution in terms of the Weierstrass functions.

§ WEIERSTRASS ELLIPTIC AND RELATED FUNCTIONS

Weierstrass elliptic and related functions are a group of special functions that were studied and introduced by Karl Weierstrass at the end of the 19th century as an improvement over the Jacobian elliptic functions $\mbox{sn}$, $\mbox{cn}$ and $\mbox{dn}$. Today it is accepted [8] that they constitute a superior tool to construct a generic theory of elliptic functions and that they are often advantageous to solve integrals in the form:

\int f\left(x, \sqrt{P(x)}\right) dx

where $f$ is a function of $x$ and of the square root of a third or fourth order polynomial $P$ not resolved into factors (Byrd[9] pag. 1889) (if the polynomial is resolved into known fixed factors, the Jacobian approach often offers a valid alternative). We report briefly the definition of these functions and a few theorems that establish their fundamental relation to fundamentals problems in astrodynamics. We follow and use the conventions and developments discussed extensively in the on-line version of the NIST Handbook of Mathematical Functions[10].

§.§ Definition

(Origin)   at (0,0);
(XAxisMin) at (-3,0);
(XAxisMax) at (5,0);
(YAxisMin) at (0,-2);
(YAxisMax) at (0,5);
[thin, gray,-latex] (XAxisMin) – (XAxisMax);
[thin, gray,-latex] (YAxisMin) – (YAxisMax);

(-5,-5) rectangle (4cm,4cm);

(Bone) at (0,1);
(Btwo) at (1,-1);
[style=help lines,dashed] (-14,-14) grid[step=2cm] (14,14);

in -7,-6,...,7
in -7,-6,...,7
[draw=gray,circle,inner sep=1pt,fill=gray] at (2*,2*) ;

[thick,-latex,red] (Origin)
– (Bone) node [above left] $\omega_3$;
[thick,-latex,red] (Origin)
– ($(Bone)+(Btwo)$) node [below right] $\omega_1$;
[thick,-latex,red] (Origin)
– (-1,-1) node [left] $\omega_2$;

[fill=gray, fill opacity=0.3, draw=black] (Origin)
rectangle ($2*(Bone)+(Btwo)$);
[fill=gray, fill opacity=0.3, draw=black] (Origin)
rectangle ($4*(Bone)+2*(Btwo)$);

[] at (1,-2.5) $\Delta > 0$;

(Origin)   at (0,0);
(XAxisMin) at (-3,0);
(XAxisMax) at (5,0);
(YAxisMin) at (0,-2);
(YAxisMax) at (0,5);
[thin, gray,-latex] (XAxisMin) – (XAxisMax);
[thin, gray,-latex] (YAxisMin) – (YAxisMax);

(-5,-5) rectangle (4cm,4cm);

(Bone) at (0,1);
(Btwo) at (1,-1);
[style=help lines,dashed] (-14,-14) grid[step=2cm] (14,14);

in -7,-6,...,7
in -7,-6,...,7
[draw=gray,circle,inner sep=1pt,fill=gray] at (2*,2*) ;

[thick,-latex,red] (Origin)
– (Bone) node [above left] $\omega_3$;
[thick,-latex,red] (Origin)
– (-1,-2) node [below right] $\omega_2$;
[thick,-latex,red] (Origin)
– (1,1) node [midway, below] $\omega_1$;

[fill=gray, fill opacity=0.3, draw=black] (Origin)
rectangle ($2*(Bone)+(Btwo)$);
[fill=gray, fill opacity=0.3, draw=black] (Origin)
rectangle ($4*(Bone)+2*(Btwo)$);

[] at (1,-3) $\Delta < 0$;

The two possible lattices in $\mathbb{C}$ having real invariants $g_2$ and $g_3$

Consider any pair of complex numbers $\omega_1, \omega_3 \in \mathbb C$ such that $\Im(\omega_3 / \omega_1) > 0$. This last request makes sure that one can rotate counterclockwise $\omega_1$ until an overlap to $\omega_3$ spanning less than $180$ degrees. The set of points $2n\omega_1 + 2m\omega_3$ with $n,m\in \mathbb{Z}$ define a lattice $\mathbb{L}$ in the complex plane. The quantities $2\omega_1$ and $2\omega_3$ are called lattice generators and are not unique. If, for example,  $\omega_1 + \omega_2 + \omega_3 = 0$, then $2\omega_3, 2\omega_2$ and $2\omega_2, 2\omega_1$ are also generator of the same lattice $\mathbb{L}$. Weierstrass defined the following function:

\wp(z|\mathbb L) = \frac{1}{z^2} + \sum_{w\in\mathbb L \backslash \left\{0\right\}}\left( \frac{1}{(z-w)^2} - \frac{1}{w^2} \right)

where the series is uniformly and absolutely convergent, and thus the exact order of its terms is irrelevant. When not needed, the underlying lattice $\mathbb L$ is omitted from the notation. From the above definition it follows that $\wp(z+2\omega_i) = \wp(z)$, hence $\wp$ is a doubly periodic function in $\mathbb{C}$, that is $\wp$ is an elliptic function. In a similar way, Weierstrass introduced two more functions, $\sigma(z|\mathbb{L})$ and $\zeta(z|\mathbb{L})$, quasi-periodic, thus not elliptic, and having the following differential relations to $\wp$:

\begin{array}{l}
\wp(z) = \zeta'(z) \\
\zeta(z) = \sigma'(z) / \sigma(z)
\end{array}

The following quantities are the so-called lattice invariants:

\begin{array}{l}
g_2 = 60 \sum_{w\in\mathbb{L}\backslash \left\{0\right\}} w^{-4}\\
g_3 = 140 \sum_{w\in\mathbb{L}\backslash \left\{0\right\}} w^{-6}
\end{array}

they are constants defined as a sum over all the lattice points except the origin and are thus determined solely by the lattice itself. Each couple of lattice generators (for example $\omega_1$, $\omega_3$) thus determine univocally the lattice invariants. Conversely, given a couple of complex numbers, there is only one lattice $\mathbb{L}$ having them as invariants. The polynomial $g(w) = 4w^3 - g_2 w -g_3$, factorized as $g(w) = (w-e_1)(w-e_2)(w-e_3)$, defines the lattice roots $e_i$. Given any pair of generators $w_1, w_3$ and $w_2 = -w_1-w_3$, the lattice roots can be ordered and identified through the relation:

\wp(\omega_i) = e_i

which we will assume valid in all of the following developments. In most of the applications that we are interested in, the lattice invariants are real numbers. From the lattice invariants definition it is immediate to see that when a lattice is symmetric with respect to the real axis then necessarily $g_2, g_3 \in \mathbb{R}$. It is possible to show that this is also a sufficient condition so that only two types of lattices will be possible resulting in real lattice invariants and are visualized in Figure <ref>. The following relations derive from the identity $4w^3 - g_2 w -g_3 = (w-e_1)(w-e_2)(w-e_3)$ and link the lattice roots to its invariants :

\begin{array}{l}
e_1+e_2+e_3 = 0 \\
g_2 = -4(e_1e_2+e_1e_3+e_2e_3) \\
g_3 = 4e_1e_2e_3
\end{array}

The discriminant $\Delta = g_2^3-27g_3^2$ determines whether the roots will be all real and distinct ($\Delta>0$) or one real and two complex conjugates ($\Delta < 0$), as well as the lattice type (see Figure <ref>). Note that under the selected convention $\omega_1$ is always real and $\omega_3$ is either a pure imaginary number ($\Delta > 0$, the lattice roots are all real) or a complex quantity with positive imaginary part ($\Delta < 0$, only one lattice root is real). The fundamental result revealing the importance of lattice invariants is the differential identity:

\wp'^2 = 4\wp^3 - g_2\wp-g_3

which also implies the important integral definition for $\wp$,

z = \int_{\wp(z)}^\infty (4s^3 - g_2s - g_3)^{- \frac 12}ds

Plot of the Weierstrass functions $\wp(x), \zeta(x)$ and $\sigma(x)$ for $g_3=0$ (lemniscatic case), $g_2=[0.1,0.2,0.5,0.8]$ and a real $x$.

§.§ Relevance to Astrodynamics

The importance of Weierstrass elliptic and related functions to astrodynamics can be best appreciated considering the following integral definition of a function $\tau$:

\begin{equation}
\tau(x) - \tau_m = \pm \int_{x_m}^{x} \frac{1}{\sqrt{a_0s^4 + 4a_1s^3 + 6a_2 s^2 + 4a_3 s + a_4}} ds
\label{eq:weierstrass1}
\end{equation}

which appears in many fundamentals problems of astrodynamics as the relation between a pseudo-time $\tau$ (or an anomaly) and a state variable $x$. A generic procedure to solve the above integral and find $x(\tau)$ rather than $\tau(x)$ is well described in the classic book from Whittaker and Watson[11] (see 20.6), and results, when $a_0$ and $a_1$ are not both null and the polynomial $f = a_0s^4 + 4a_1s^3 + 6a_2 s^2 + 4a_3 s + a_4$ has no repeated factors, in the use of the following inversion formula:

\begin{equation}
x(\tau) = x_m + \frac{A}{\wp(\tau - \tau_m, g_2, g_3)-B}
\label{eq:weierstrass2}
\end{equation}

where $x_m$ is a root of $f$, $A = 1/4 f'(x_0)$,  $B = 1 / 24 f''(x_0)$ and the two lattice invariants are:

\begin{equation}
\begin{array}{l}
g_2 = a_0a_4 - 4a_1a_3+3a_2^2 \\
g_3 = a_0a_2a_4+2a_1a_2a_3-a_2^3-a_0a_3^2-a_1^2a_4
\end{array}
\label{eq:lattice_invariants}
\end{equation}

The quantity $\tau_m$ is what we call time of root passage since $x(\tau_m) = x_m$, and can be computed writing Eq.(<ref>) at $\tau=0$ and finding $\wp$ from it:

\begin{equation}
\wp(\tau_m) = B + \frac{A}{x_0 - x_m}
\label{eq:timeofrootpassage}
\end{equation}

where $x_0$ is the initial value $x(0)$. The inversion of the Weierstrass function will return two valid values for $\tau_m$ which reflect the original ambiguity in the integral sign. Such an ambiguity is solved forcing the initial condition. The derivative of Eq.(<ref>) with respect to $\tau$ is:

\begin{equation}
x'(\tau) = -\frac{(x(\tau) - x_m)^2}{A}\wp'(\tau - \tau_m)
\label{eq:weierstrass3}
\end{equation}

which, for $\tau=0$, holds:

\wp'(\tau_m) = \frac{Ax'_0}{(x_0-x_m)^2}

which can be used to univocally determine $\tau_m$ from the initial condition $x_0'$.
It is noteworthy that any root $x_m$ of the polynomial $f(s)$ can be used with the above formulae be it a real, pure imaginary or complex root. Choosing a real root, when possible, has two main advantages: it allows to define the origin of the $\tau$ pseudo-time variable so that $\tau_m=0$, and it keeps all the computations in the real domain (i.e. no complex quantities involved) resulting in a significant increase in the efficiency of evaluating the expressions on a computer as shown in the next section.

§.§ Notes on the computer implementation

Function     type     evaluations     ns / evaluation

real $\sin$    double                       20000000             62

real $\wp$    double                        10000000            178

real $\wp'$    double                       5000000            215

real $\zeta$    double                      10000000            166

real $\sigma$    double                     10000000            175

complex $\sin$    double                    10000000            175

complex $\wp$    double                     2000000           532

complex $\wp'$    double                    2000000           787

complex $\wp^{-1}$    double                1000000           1643

complex $\zeta$    double                   2000000           531

complex $\sigma$    double                  2000000           659

Performances of Weierstrass elliptic and related functions as implemented in the project w_elliptic. Results for the standard sin function (std::sin) are also shown for comparison.. The table has been obtained running the project test suite compiled using gcc using the -ffast-math flag on a Intel(R) Core(TM) i7-3610QM CPU.

The computer implementation of Weierstrass elliptic and related functions received little or no attention from the computer science community in the last decades. As a consequence,  not many languages nor tools are offering the possibility to compute these functions and it is very difficult to assess their use in terms of computational efficiency, for example with respect to the much more studied and popular Jacobi's elliptic functions. An attempt was recently made, limited to the Stark problem, by Hatten and Russel[12] who, not having access to efficient implementations of the Weierstrass functions nor of the expressions using them, were forced to conclude that Weierstrass functions carry a computational penalty with respect to other methods. Such a conclusion was also later used to justify the methodology adopted by Beth et al. in their study on planetary exospheres[13, 14]. In reality, while it is true that the construction of the lattice from the invariants $g_2, g_3$ is a necessary step whose cost must be paid, once this step is performed the actual evaluation of the Weierstrass elliptic and related function is extremely fast as shown in Table <ref> where their speed is compared to the computational speed of the simplest trigonometric function (std::sin). In order to build such a table we programmed a C++ (and python) open source project called w_elliptic (<https://github.com/bluescarni/w_elliptic>) that implements efficient versions of these functions distinguishing with respect to the argument type being a complex or a real (double) number. It is not in the scope of this paper to discuss the implementation details of w_elliptic, which is still undergoing further optimization and improvements and is the subject of a dedicated paper under preparation[15]. It is, though, very clear that these functions can be computed extremely efficiently, especially in the real domain. We will see in the rest of this paper how indeed most of the expressions involved in the solution of fundamental astrodynamical problems can be written as to keep the argument of the Weierstrass functions in the real domain. Any use of explicit solutions in terms of Weierstrass functions using their complex implementation is bound to be much slower as shown in the Table.

§ SOLUTION TO THE CONSTANT RADIAL ACCELERATION INITIAL VALUE PROBLEM

Previous work [5] reported the fundamental theoretical developments that lead to solve explicitly the radial acceleration problem. We here use those developments to establish a procedure to find the spacecraft position and velocity  $\mathbf r, \mathbf v$ at any time $t$.

Define the motion invariants, i.e. the angular momentum and the specific energy:

\begin{equation}
h = r_0^2 \dot\theta = \mathbf r_0 \times \mathbf v_0
\label{eq:cons_h}
\end{equation}

\begin{equation}
\mathcal E = \frac{v_0^2}2 - \frac \mu{r_0} - \alpha r_0
\label{eq:vis_viva}
\end{equation}

Note that $\alpha > 0$ corresponds to an outward pointing acceleration.
Compute the two invariants of the Weierstrass elliptic and related functions in the constant radial acceleration case:

\begin{equation}
\begin{array}{l}
g_2 = \frac{\mathcal E^2}{3} - \alpha \mu \\
g_3 = \frac{\alpha^2}{4}(h^2 + \frac{2 \mathcal E \mu}{3\alpha} - \frac{4\mathcal E^3}{27\alpha^2})
\end{array}
\end{equation}

Define the following third order polynomial and its derivatives:

f(r) = 2\alpha r^3 +2 \mathcal E r^2 + 2 \mu r - h^2

f'(r) = 6\alpha r^2 +4 \mathcal E r + 2\mu

f''(r) = 12\alpha r +4 \mathcal E

Compute the three roots $\tilde e_1, \tilde e_2, \tilde e_3$ of $f$. Compute the radius of pericenter passage $r_m$ as the real root closest to $r_0$ and such that $r_m \le r_0$. Using the energy conservation equation define the velocity of pericenter passage $v_m^2 = 2(\mathcal E + \frac \mu {r_m} + \alpha r_m)$. Compute the constants $B = \frac 1{24} f''(r_m)$, $A = \frac 14 f'(r_m)$ and the quantity $\xi$ defined as $\wp(\xi) = (B - \frac{A}{r_m})$, $\wp'(\xi) = \frac{Ah}{r_m^2} \imath $. Note that $\xi$ will be a real number if $A < 0$.

Introducing the radial anomaly $\tau(t)$ (via a Sundmann transformation $\frac{dt}{d\tau} = r$), the final explicit solution is given by:

\begin{equation}
\label{eq:radial1}
\left\{
\begin{array}{l}
r(\tau) = r_m + \frac{A}{\wp(\tau) - B} \\
\sin \theta(\tau) = x(\tau) \sin(v_m\tau) - y(\tau) \cos(v_m \tau) \\
\cos \theta(\tau) = y(\tau ) \sin(v_m\tau) + x(\tau) \cos(v_m \tau)
\end{array}
\right.
\end{equation}

where we defined $x(\tau) + iy(\tau) = \frac{\sigma(\dario-\tau)}{\sigma(\tau+\dario)} \exp{2\tau\zeta(\dario)}$. We also have:

\begin{equation}
\label{eq:radial2}
\left\{
\begin{array}{l}
\dot r(\tau) = - \frac{(r(\tau) - r_m)^2}{Ar(\tau)} \wp'(\tau) \\
\dot \theta(\tau) = \frac{h}{r^2(\tau)}
\end{array}
\right.
\end{equation}

The relation between the time and the radial anomaly is described via the radial Kepler's equation which admits the following equivalent forms:

\begin{equation}
t(\tau) = r_m \tau  - \frac{ 4 B A}{g_3+ 8 B^3}\left[ B\tau +
\zeta(\tau) + \frac 12 \frac{\wp'(\tau)}{\wp(\tau) - B} \right]
\label{eq:radial_kepler}
\end{equation}

\begin{equation}
t(\tau) = r_m \tau  - \frac{ 2 B A}{g_3+ 8 B^3}\left[ 2 B\tau +
\zeta(\tau-\omega_k) + \zeta(\tau + \omega_k) \right]
\label{eq:radial_kepler1}
\end{equation}

where $\wp(\omega_k) = B$. Note that in case the motion is bounded, necessarily $\wp(\tau) > B, \forall \tau \in R$, hence $\omega_k$ will be a complex quantity, while in case of unbounded motion $\omega_k$ will necessarily be a real quantity and thus the second expression, as the first one, would also involves only real quantities while is not undefined at $\tau = 0$. Note that both radial and true anomaly $\tau$, $\theta$ are zero at a pericenter passage, i.e. $t=0 \rightarrow \tau=0$, $\theta = 0$, $r = r_m$ . The radial anomaly at the initial conditions is found computing the expression:

\begin{equation}
\wp(\tau_0) = B - \frac{A}{r_m-r_0}
\label{eq:tau0}
\end{equation}

while the true anomaly at the initial conditions is found computing the expression:

\begin{equation}
\begin{array}{l}
\sin \theta_0 = x(\tau_0) \sin(v_m\tau_0) - y(\tau_0) \cos(v_m \tau_0) \\
\cos \theta_0 = y(\tau_0) \sin(v_m\tau_0) + x(\tau_0) \cos(v_m \tau_0)
\end{array}
\label{eq:theta0}
\end{equation}

The explicit Cartesian coordinates of the satellite can be finally obtained, trivially, using the following expression:

\begin{equation}
\begin{array}{l}
\mathbf r(t) = r(t)\cos\theta(t) \mathbf i_{r_m} + r(t)\sin\theta(t)\mathbf i_{p_m} \\
\mathbf v(t) = (\dot r(t)\cos\theta(t) - \frac h{r(t)} \sin\theta(t)) \mathbf i_{r_m} + (\dot r(t)\sin\theta(t) +\frac h{r(t)} \cos\theta(t))\mathbf i_{p_m}
\end{array}
\end{equation}

which involves the solution to the Kepler's radial equation in order to get $t$ from the radial anomaly $\tau$.

§ RADIAL FLY-BY

Planetary fly-bys are often modelled as Keplerian hyperbolas. In the preliminary mission design phases it is common to consider the effect of a fly-by as that of an instantaneous rotation of the relative velocity vector by the angle $\delta_K$ function of the closest passage distance $r_m$ and of the hyperbolic trajectory plane orientation. The outgoing relative velocity vector is then summed to the planet velocity vector to obtain the new spacecraft state in the interplanetary medium. The spacecraft is generally considered to be not thrusting during this phase. Here we study a simple powered fly-by manoeuvre where the spacecraft, along its planetocentric hyperbola, turns on its propulsion system when its distance $r$ from the planet is in the interval $[r_o, r_i]$ to produce a constant acceleration of magnitude $\alpha$. Since the resulting trajectory is perfectly symetric, also the effect of such a fly-by can be considered as a rotation of the relative velocity vector by an angle $\delta$. Using the analytical explicit solution in terms of the Weierstrass elliptic and related functions allow to derive simple equations to find $\delta$ as a function of $r_m, V_\infty, r_i, r_o$ and $\alpha$. In Figure <ref> the basic geometry of a radial fly-by is shown. Note that the trajectory is not an hyperbola and is obtained by patching hyperbolic arcs (outside $[r_i, r_o]$) with the constant radial acceleration solution (in $[r_i, r_o]$).

[draw=none,fill=gray, fill opacity = 0.1] (-*,0) circle (0.7);
[draw=none,fill=white] (-*,0) circle (0.4);

[->] (-*-,0,0) – (0.25-*-,0,0) node[right] $\mathbf i_{r_m}$;
[->] (-*-,0,0) – (-*-,0.25,0) node[left] $\mathbf i_{p_m}$;

[] (-*-,0) – (-*cosh(-1),*̱sinh(-1)) node[midway, above, fill=white] $r_o$;

[] (-*-,0) – (-*cosh(-0.725),*̱sinh(-0.725)) node[midway, right, fill=white] $r_i$;

[arrows=<->](-*-,-0.02)–(0.08-*-,-0.02) node[midway, below] $r_m$;

[red] plot[domain=-1:-0.725] (-*cosh(),*̱sinh());
[red] plot[domain=0.725:1] (-*cosh(),*̱sinh());

[->] plot[domain=-1.3:-1] (-*cosh(),*̱sinh());
plot[domain=1:1.3] (-*cosh(),*̱sinh());

[->] plot[domain=-0.725:0.725] (-*cosh(),*̱sinh());

[dashed] plot[domain=-1.2:0] (1.65*,0.85*+0.22);
[dashed] plot[domain=-1.2:0] (1.65*,-0.85*-0.22);

[thick,->] (1.65*1,0.85*1+0.22) – (1.65*1,0.85*1+0.22) node[above] $\mathbf v^-_\infty$;
[thick,->] (1.65*2,-0.85*2-0.22) – (1.65*2,-0.85*2-0.22) node[above] $\mathbf v^+_\infty$;

[->,domain=30:150] plot (-1.65*0.22/0.85+0.1*cos(), 0.1*sin()) node[right] $2\delta$;

The radial fly-by geometry. In the range $r\in[r_i, r_o]$, the spacecraft maintains an additional radial acceleration $\alpha$ and is thus not flying along a  hyperbola

§.§ The unpowered flyby case $\alpha=0$

First, as a benchmark, consider the case of a purely ballistic fly-by modelled as a Keplerian hyperbola. Assume as incoming conditions $v_{\infty}$, and write the specific energy at the incoming conditions at infinite:

\mathcal E_K = \frac{v_\infty^2}{2}

and at the closest distance:

v_m^2 = \frac{2\mu}{r_m} + v_\infty^2

where the subscript $K$ indicates that we are in the purely Keplerian case. We compute the angular momentum as $h_K = r_m v_m$:

h_K = \sqrt{r_m(2\mu+v_\infty^2 r_m)}

and the relative velocity rotation half angle as:

\begin{equation}
\sin \delta_K = \frac 1e
\label{eq:unpowered}
\end{equation}

where, accounting that $a(1-e^2) = \frac{h^2}{\mu}$, we have:

e^2 = 1 + \frac{2h_K^2 \mathcal E_K}{\mu^2}

so that from $r_m, v_\infty$ we can compute $\delta_K$.

§.§ The powered flyby case

We now study how the relative velocity rotation angle $\delta_K$ is modified when we assume a constant acceleration $\alpha$ acting on the spacecraft at $r \in [r_i, r_o]$ while keeping $r_m$ unchanged. From Figure <ref> it can be seen how $2 \delta = 2 (\delta_1 + \delta_2 + \delta_3)$ is the sum of three contributions. The first and the last one, indicated with $2\delta_1$ and $2\delta_3$ are due to the hyperbolic motion outside the $[r_i, r_o]$ interval and can be computed using the Keplerian solution. The second contribution, indicated with $2 \delta_2$ is due to the radial accelerated motion and must thus be computed using Weierstrass functions. Let us preliminary determine the various motion invariants for the different arcs. Starting with the out-most hyperbolic arcs we have the specific energy:

\mathcal E_{K_o} = \frac{v_\infty^2}{2}

and by its conversation:

v_o^2 =2 \left(\mathcal E_{K_o} + \frac \mu {r_o}\right)

hence the specific energy along the following propelled arc will be:

\mathcal E = \mathcal E_{K_o} - \alpha r_o

and by its conversation:

v^2_i = 2\left(\mathcal E_{K_o}+\frac\mu{r_i} +\alpha(r_i-r_o)\right)

hence the energy along the inner hyperbolic arc will be:

\mathcal E_{K_i} = \mathcal E_{K_o}+\alpha(r_i-r_o)

and the velocity in the point of closest passage:

v_m^2 = 2\left(\mathcal E_{K_o}+\frac\mu{r_m} +\alpha(r_i-r_o)\right)

which allow to compute the angular momentum for all the arcs as $h(v_\infty, r_m, r_i, r_o) = r_m v_m$

§.§.§ Computing $\delta_1$ and $\delta_3$

To compute $\delta_1$ and $\delta_3$ we develop a generic expression for the angle $\delta$ between the velocities acquired along a Keplerian arc at two positions $r_1$ and $r_2$. Along a Keplerian hyperbola we have (see Battin[16] (3.6)):

\frac{h_K}{\mu} \mathbf v =  - \sin \theta_K \mathbf{\hat i}_e + (e + \cos\theta_K) \mathbf{\hat i}_p

where $\theta_K$ is the true anomaly in the Keplerian motion, i.e counted with respect to $ \mathbf{\hat i}_e$. We will drop the subscript $K$ for the true anomaly in the following expressions to avoid cluttering our notation. We thus have:

\frac{p}{\mu}(\mathbf v_1 \cdot \mathbf v_2) = \frac{p}{\mu}v_1 v_2 \cos\delta = \sin\theta_1\sin\theta_2 + (e+\cos\theta_1)(e+\cos\theta_2)

where the orbit polar $r = p / (1+e\cos\theta)$ can be used to compute:

\begin{array}{ll}
\cos\theta = \frac 1e (\frac p{r} - 1) &
\sin\theta = -\frac 1e\sqrt{e^2-(\frac p{r} - 1)^2}
\end{array}

where the sign minus is chosen on the sine as the hyperbolic arc considered is incoming. Applying the above formulae to compute $\delta_1$ and $\delta_3$ we get:
We finally have the the final expression for $\delta_1$:

\frac{p}{\mu}v_\infty v_o \cos\delta_1 = \frac 1{e^2}\sqrt{(e^2-1)\left(e^2-\left(\frac p{r_o}-1\right)^2\right)}+\left(e-\frac 1e\right)\left(e + \frac 1e \left(\frac p{r_o} - 1\right)\right)

and for $\delta_3$:

\frac{p}{\mu}v_i v_m \cos\delta_3 = \left(e + \frac 1e \left(\frac p{r_i}-1\right)\right)(e+1)

Note that in inverting the above expressions using $\arccos$ we have $\delta_1, \delta_3 \in [0, \pi]$ which is correct whenever no passage through the pericenter happens.

§.§.§ Computing $\delta_2$

While $r_i < r < r_o$ the propulsion system of the spacecraft turns on putting the spacecraft on a radially accelerated trajectory having $\mathcal E, h$ as motion invariants. We then may also compute $r_m$ and the lattice invariants $g_2$ and $g_3$ and hence all the relevant Weierstrass elliptic and related functions.
We may then compute:

\delta_2 = (\gamma_o - \theta_o) - (\gamma_i - \theta_i)

where the true anomalies $\theta < 0$ are now referred to the non Keplerian arc and are computed by Eq.(<ref>), while $\gamma \in [-\pi,0]$ is the flight path angle which can also be computed along the Keplerian arcs. We use:

\tan\gamma = - \frac {r}{p}\sqrt{e^2-\left(\frac p{r} - 1\right)^2}

The final effect of a radial fly-by will then be to rotate the relative velocity vector by an angle $2\delta = 2(\delta_1+\delta_2+\delta_3)$ by providing a cumulative $\Delta V_{lt} = 2 \alpha \cdot \mbox{tof}$ where $\mbox{tof} = t(\tau_o) - t(\tau_i)$ is the duration of the first propelled arc as computed applying Eq.(<ref>) in correspondence of the two radial anomalies at $r_i$ and $r_o$.

§.§.§ A numerical example

Effect of a radial flyby at the Moon for $r_o=40 R_M$, $r_m = 1.1R_M$ and $m \alpha = [8, 7, 6, 5, 4, 3, 2]$ [N] with $m=2000$ [kg]. The net $\Delta V$ gained by the manoeuvre (left) is plotted aside the $\Delta V_{lt}$ necessary to perform the manoeuvre (right).

Consider a Moon fly-by, where the spacecraft approaches the sphere of influence with a relative velocity $v_\infty = 1000$ [m/s] and performs a fly-by with closest approach distance $r_m = 1.1 R_M$, where $R_M = 4905$ [km] is the Moon radius. Under these conditions the effect of an unpowered fly-by, as computed from Eq.(<ref>), is that of rotating the relative velocity vector by an angle $2 \delta_K = 1.209$ [rad.]. We study the possibility to increase such an angle (the new value indicated by $2 \delta$) by performing a powered fly-by during which the spacecraft maintains a constant outward acceleration of magnitude $\alpha$ when its distance from the atracting body is $r\in[r_i, r_o]$. We consider to start the propelled phase at $r_o = 40 R_M$ and consider $r_i\in[1.1,40]R_M$. We compare the $\Delta V_{lt}$ used by the low-thrust propulsion system to the instantaneous $\Delta V_{K}$ that would be needed at the end of the outgoing asymptote to change the relative velocity direction by the same amount: $\Delta V_{K} = 2 v_\infty \sin{(\delta-\delta_K)}$. The angle $\delta$ as well as the $\Delta V_{lt}$ are computed using the formulae developed above. The use of the newly developed formulae enables to make this study very efficiently avoiding numerical propagation altogether. In Figure <ref> we plot the results in the selected case. We show the net gain of $\Delta V$ computed as the difference between $\Delta V_K$ and $\Delta V_{lt}$ as well as the value of $\Delta V_{lt}$. Note the area where a $\Delta V$ amplification effect is present. This is related to the decrease in spacecraft velocity which allows for the planet gravity to bend the asymptotes with greater efficiency. The analytical formulae derived express the deflection angle $\delta$ and the velocity increment $\Delta V_{lt}$ as an explicit function of $r_i, r_o, \alpha, v_\infty$. They are suitable to be used in a larger interplanetary trajectory optimization scheme, as well as in the preliminary assessment of some planetary encounter.

§ SOLUTION TO THE 2D STARK INITIAL VALUE PROBLEM

Previous work [6] reported the explicit solution to the Stark problem in the full three dimensional case. That solution can be simplified if one restricts the problem to be purely two-dimensional. We present the derivation of the Stark problem solution in terms of Weierstrass functions specific for this simpler case. Note that an explicit solution using Jacobian elliptic functions is known for the 2D Stark case [17]. The reader is encouraged to compare the Weierstrass form of such a solution derived here to the solutions needed to cover all possible cases using Jacobi elliptic integrals. As recognized by Byrd[9], the Weierstrass approach has a clear advantage when the polynomial roots and their order are unknown, which is the case here as initial conditions will determine such roots hierarchy.

Consider the planar motion of a spacecraft subject to an inertially fixed acceleration of magnitude $\alpha > 0$ directed along the $x$ axis. The system specific energy, that is conserved, can be written as:

\begin{equation}
\mathcal E = \frac {v^2}2 - \frac \mu r - \alpha x
\label{eq:energystark}
\end{equation}

In a similar way as done in the full three dimensional case, we introduce the coordinates $\xi \in [-\infty, \infty]$, $\eta \in [-\infty, \infty]$ via the following transformations:

\begin{array}{lr}
x = \frac{\xi^2 - \eta^2}2 & \dot x = \xi\dot\xi - \eta\dot\eta \\
y = \xi\eta & \dot y = \eta\dot\xi + \xi\dot\eta
\end{array}

Note that with respect to the classical parabolic coordinates we do not restrict the domain of either $\eta$ or $\xi$ to the real axis. As a consequence, the transformation here used is not unique when reversed:

\begin{equation}
\begin{array}{lr}
\xi = \pm \sqrt{r + x} & \dot\xi = \pm \frac{\dot r + \dot x}{2\sqrt{r+x}} \\
\eta = \pm \sqrt{r - x}   & \dot\eta = \pm \frac{\dot r - \dot x}{2\sqrt{r-x}}
\end{array}
\label{eq:xietadef}
\end{equation}

where $r = \sqrt{x^2 + y^2} = \frac{\xi^2 + \eta^2}2$ and $\dot r = \frac{x\dot x + y\dot y}{r} = \dot\xi \xi + \dot\eta \eta$. This is not a problem here, rather an advantage, as we will only be interested in having unique values for $x$ and $y$. We may then choose any of the signs above when, for example, computing the initial conditions $\xi_0$ and $\eta_0$.

Introduce now two further motion invariants $h_\xi$ and $h_\eta$[6]:

\begin{array}{l}
h_\xi = -\alpha \frac {\xi^4}{2} - \mathcal E \xi^2 + 2 r^2\dot\xi^2 \\
h_\eta = \alpha \frac {\eta^4}{2} - \mathcal E \eta^2 + 2 r^2\dot\eta^2
\end{array}

linked by the relation $h_\xi+h_\eta = 2\mu$. We may then write the fundamental differential equations that allow to solve the 2-D Stark problem:

\begin{equation}
\begin{array}{l}
\frac{d\xi}{d\tau} = \pm \sqrt{\alpha \xi^4 + 2\mathcal E\xi^2 + 2h_\xi} = \pm \sqrt{f_\xi}\\
\frac{d\eta}{d\tau} = \pm \sqrt{-\alpha \eta^4 + 2\mathcal E\eta^2 + 2h_\eta}  = \pm \sqrt{f_\eta}
\end{array}
\label{eq:differential_equations}
\end{equation}

where the pseudo-time $2r d\tau = dt$ is used.

Plot of the Stark 2D solution as obtained from Eq.(<ref>) and Eq.(<ref>) for an unbounded case (above) and a bounded case (below). Cartesian coordinates are shown together with one of the two possible choices for the parabolic ones. Initial conditions are $\mathbf r_0 = [-1.,1e-3]$, $\mathbf v_0=[1e-3,1.5]$, $\alpha = 0.1$ for the unbounded case, and $\mathbf r_0 = [1e-3,1]$, $\mathbf v_0=[1, 1e-3]$, $\alpha = 0.003$. Non dimensional units are used so that $\mu=1$.

§.§ Roots of $f_\xi$

It is straight forward to compute the roots of the polynomial $f_\xi$ applying the quadratic equation root formula on the fictitious variable $s = \xi^2$:

s_{1,2} = \frac{-\mathcal E \pm \sqrt{\mathcal E^2-2\alpha h_\xi}}{\alpha}

which translates immediately to the four solutions in $\xi$:

\begin{equation}
\tilde e^\xi_{1,2,3,4} = \pm \sqrt{\frac{-\mathcal E \pm \sqrt{\mathcal E^2-2\alpha h_\xi}}{\alpha}}
\label{eq:rootsxi}
\end{equation}

These roots can be, depending from the initial conditions, all complex, all real or two real and two pure imaginary.

§.§ Roots of $f_\eta$

Likewise, we compute the roots of the polynomial $f_\eta$ applying the quadratic equation root formula on the fictitious variable $s = \eta^2$:

s_{1,2} = \frac{\mathcal E \pm \sqrt{\mathcal E^2+2\alpha h_\eta}}{\alpha}

which translates immediately to the four solutions in $\eta$:

\begin{equation}
\tilde e^\eta_{1,2,3,4} = \pm \sqrt{\frac{\mathcal E \pm \sqrt{\mathcal E^2+2\alpha h_\eta}}{\alpha}}
\label{eq:rootseta}
\end{equation}

Note that $\lim_{\eta\to\pm\infty} f(\eta) = -\infty$ and $f(\eta_0) > 0$ since $\dot\eta_0$ must be real. As a consequence at least one of the above roots (and hence two) must be real.

§.§ Solution for the $\xi$ coordinate

Consider the first of Eqq.(<ref>). Rewrite it as follow:

\tau - \tau_{m,\xi} = \pm \int_{\xi_m}^\xi \frac {ds}{\sqrt{\alpha s^4 + 2\mathcal Es^2 + 2h_\xi}}

where $\xi_m$ is any one of the four roots in Eq.(<ref>).
We recognize the expression has the form  Eq.(<ref>) and we thus
introduce the lattice invariants, as defined in Eq.(<ref>):

\begin{equation}
\begin{array}{l}
g_{2, \xi} = 2\alpha h_\xi + \frac{\mathcal E^2}{3}\\
g_{3, \xi} = \frac {\mathcal E}{3}\left( 2\alpha h_\xi - \frac {\mathcal E^2}9\right)
\end{array}
\label{eq:xiinv}
\end{equation}

and find the lattice roots as:

e_{1,\xi} = - \frac {\mathcal E} 3, \quad
e_{2,\xi} = \frac {\mathcal E} 6 + \frac{\sqrt 2}{2}\sqrt{\alpha h_\xi}, \quad
e_{3, \xi} = \frac {\mathcal E} 6 - \frac{\sqrt 2}{2}\sqrt{\alpha h_\xi}

Apply now Eq.(<ref>) to write:

\begin{equation}
\xi(\tau) = \xi_m + \frac{A_\xi}{\wp_\xi(\tau-\tau_{m,\xi})-B_\xi}
\label{eq:xi}
\end{equation}

where $A_\xi = 1 / 4 f'(\xi_m) = \xi_m(\alpha \xi_m^2 + \mathcal E)$, $B_\xi = 1 / 24 f''(\xi_m) = \frac{1}{2}(\alpha\xi_m^2 + \frac 13 \mathcal E)$ and $\wp_\xi(\tau_{m,\xi}) = B_\xi + \frac{A_\xi}{\xi_0 - \xi_m}$, $\wp_\xi'(\tau_{m,\xi}) = \frac{2 r_0 A_\xi \dot\xi_0 }{(\xi_0 - \xi_m)^2}$. We also have:

\begin{equation}
\dot\xi(\tau) = - \frac{(\xi(\tau) - \xi_m)^2}{2rA_\xi}\wp_\xi'(\tau - \tau_{m,\xi})
\label{eq:dxi}
\end{equation}

The expressions derived are valid regardless of the choice of the root $\xi_m$, but as explained previously, choosing a real root is desirable. The conditions to have at least one real root for the polynomial $f$, as easily verified from Eq.(<ref>), can be written as:

\begin{equation}
\begin{array}{l}
a) \quad \mathcal E>0, \alpha h_\xi <0 \\
b) \quad \mathcal E<0, \alpha h_\xi<\frac{\mathcal E^2}{2}
\end{array}
\end{equation}

in both cases $\tilde e_1 = \sqrt{\frac{- \mathcal E + \sqrt{\mathcal E^2-2\alpha h_\eta}}{\alpha}}$ will result to be the expression for one of the real roots and we thus choose it in all cases. In essence, $h_\xi$ determines whether at least one real root exist. Noting that $\lim_{\alpha\to\infty} h_\xi = -\frac{\alpha}{2}y^2$ we can be certain that high level of thrust put us in this condition. Looking into the opposite direction, we note how low thrust levels also guarantee that, for $\mathcal E < 0$, a real root exists.

§.§ Solution for the $\eta$ coordinate

Consider the second of Eq.(<ref>). Rewrite it as follows:

\tau - \tau_{m,\eta} = \pm \int_{\eta_m}^\eta \frac {ds}{\sqrt{-\alpha s^4 + 2\mathcal Es^2 + 2h_\eta}}

where $\eta_m$ is any one of the four roots in Eq.(<ref>).
We recognize the expression has the form  Eq.(<ref>) and we thus
introduce the lattice invariants, as defined in Eq.(<ref>):

\begin{equation}
\begin{array}{l}
g_{2, \eta} = -2\alpha h_\eta + \frac{\mathcal E^2}{3}\\
g_{3, \eta} = \frac {\mathcal E}{3}\left( -2\alpha h_\eta - \frac {\mathcal E^2}9\right)
\end{array}
\label{eq:etainv}
\end{equation}

and find the lattice roots as:

e_{1,\eta} = - \frac {\mathcal E} 3, \quad
e_{2,\eta} = \frac {\mathcal E} 6 + \frac{\sqrt 2}{2}\sqrt{-\alpha h_\eta}, \quad
e_{3, \eta} = \frac {\mathcal E} 6 - \frac{\sqrt 2}{2}\sqrt{-\alpha h_\eta}

Apply now Eq.(<ref>) to write:

\begin{equation}
\eta(\tau) = \eta_m + \frac{A_\eta}{\wp_\eta(\tau-\tau_{m,\eta})-B_\eta}
\label{eq:eta}
\end{equation}

where $A_\eta = 1 / 4 f'(\eta_m) = \eta_m(-\alpha \eta_m^2 + \mathcal E)$, $B_\eta = 1 / 24 f''(\eta_m) = \frac{1}{2}(-\alpha\eta_m^2 + \frac 13 \mathcal E)$ and $\wp_\eta(\tau_{m,\eta}) = B_\eta + \frac{A_\eta}{\eta_0 - \eta_m}$, $\wp_\eta'(\tau_{m,\eta}) = \frac{2 r_0 A_\eta \dot\eta_0 }{(\eta_0 - \eta_m)^2}$.
We also have:

\begin{equation}
\dot\eta(\tau) = - \frac{(\eta(\tau) - \eta_m)^2}{2rA_\eta}\wp_\eta'(\tau - \tau_{m,\eta})
\label{eq:deta}
\end{equation}

It is possible to show that $\tilde e^\eta_1$ is always real and we will thus choose it as our $\eta_m$

§.§ The time equation

The pseudo time $\tau$ is defined via the differential equation $dt = (\xi^2 + \eta^2) d\tau$ which we may now write explicitly using Eq.(<ref>) and Eq.(<ref>):

\frac{dt}{d\tau} = \xi_m + \eta_m + \frac{A_\xi^2}{[\wp(\tau-\tau_{m,\xi})-B_\xi]^2}  + \frac{A_\eta^2}{[\wp(\tau-\tau_{m,\eta})-B_\eta]^2} + \frac{2A_\xi}{\wp(\tau-\tau_{m,\xi})-B_\xi}  + \frac{2A_\eta}{\wp(\tau-\tau_{m,\eta})-B_\eta}

In order to integrate the above equation, we employ two formulae from Tannery & Molk [18]  [Chapter CXII] (see also Gradshteĭn & Ryzhik[19] [5.141]):

\begin{align}
\int\frac{du}{\wp\left(u\right) - \wp\left(v\right)} & = \frac{1}{\wp^\prime\left(v\right)}\left[\ln\frac{\sigma\left(u - v\right)}{\sigma\left(u + v\right)}+2u\zeta\left(v\right)\right],\label{eq:molk_00}\\
\int\frac{du}{\left[\wp\left(u\right) - \wp\left(v\right)\right]^2} & =  - \frac{1}{\wp^{\prime 2}\left(v\right)} \left[
\zeta\left(u - v\right)+\zeta\left(u + v\right)
\vphantom{\int\frac{du}{\wp\left(u\right) - \wp\left(v\right)}}
\right.\notag\\
&\quad \left. + 2u\wp\left(v\right)+\wp^{\prime\prime}\left(v\right) \int\frac{du}{\wp\left(u\right) - \wp\left(v\right)}
\right].\label{eq:molk_01}
\end{align}

where it is assumed that $\wp'(v) \ne 0$, that is $v$ is not a root of the Weierstrass polynomial $g(s) = 4s^3-g_2s-g_3$.
We introduce the shorthand notation:

\begin{align}
\mathcal{J}_1\left(u,v\right) & = \int\frac{du}{\wp\left(u\right) - \wp\left(v\right)},\\
\mathcal{J}_2\left(u,v\right) & = \int\frac{du}{\left[\wp\left(u\right) - \wp\left(v\right)\right]^2}
\end{align}

(with the understanding that we will add a $\xi$ or $\eta$ subscript depending on the subscript of the Weierstrassian functions appearing in the integrals). We may then derive the following time equation:

\begin{align}
t(\tau) = &(\xi_m + \eta_m) \tau + \notag\\
&A^2_\xi \left(\mathcal{J}_{2,\xi}\left(\tau-\tau_{m,\xi},b_\xi\right)+ \mathcal{J}_{2,\xi}\left(\tau_{m,\xi},b_\xi\right) \right) + \notag\\
&A^2_\eta \left(\mathcal{J}_{2,\eta}\left(\tau-\tau_{m,\eta},b_\eta\right)+ \mathcal{J}_{2,\eta}\left(\tau_{m,\eta},b_\eta\right) \right) + \notag\\
&2A_\xi \left(\mathcal{J}_{1,\xi}\left(\tau-\tau_{m,\xi},b_\xi\right)+ \mathcal{J}_{1,\xi}\left(\tau_{m,\xi},b_\xi\right) \right) + \notag\\
&2A_\eta \left(\mathcal{J}_{1,\eta}\left(\tau-\tau_{m,\eta},b_\eta\right)+ \mathcal{J}_{1,\eta}\left(\tau_{m,\eta},b_\eta\right) \right)
\end{align}

where we have defined $b_\xi = \wp_\xi^{-1}(B_\xi)$, $b_\eta = \wp_\eta^{-1}(B_\eta)$ (note that we can take any of the values returned by the inverse $\wp$ as $\mathcal J_1(u,v) = \mathcal J_1(u,-v)$, $\mathcal J_2(u,v) = \mathcal J_2(u,-v)$).

§ STARK FLY-BY

å(*) /

[draw=none,fill=gray, fill opacity = 0.05] (0,0) circle (0.7);
[draw=none,fill=gray, fill opacity = 0.05] (0,0) circle (0.9);

[->] (0,0) – (0.3,0) node[right] $\mathbf i_{x}$;
[->] (0,0) – (0,0.3) node[left] $\mathbf i_{y}$;

[red] plot[domain=-0.94:0] (*̱sinh(), *cosh() - *);
[red, ->] plot[domain=0:0.93] (*sinh(), å*cosh() + (-å) - *);

[] plot[domain=-1.5:-0.94] (*̱sinh(), *cosh() - *);
[] plot[domain=0.93:1.5] (*sinh(), å*cosh() + (-å) - *);

[draw=none,fill=black, fill opacity = 0.5] (*̱sinh(-0.94), *cosh(-0.94) - *) circle (0.02);

[dashed] plot[domain=-1.4:0] (,*+);

[dashed] plot[domain=-0.2:1.7] (,*+);

- (-) / (-)

[thick,->] (,) – (+0.1,+*0.1) node[left] $\mathbf v^-_\infty$;
[thick,->] (,) – (+0.2,+*0.2) node[right] $\mathbf v^+_\infty$;

[<->] (-0.615,0.4) – (0,0.4) node[midway, above] $-x_0$ ;
[<->] (-0.68,0.32) – (-0.68,0) node[midway, left] $y_0$ ;
[<->] (0,-0.01) – (0,-0.14) node[midway, left] $r_m$ ;

The Stark fly-by geometry.

Consider a second case of powered flyby (which we will refer to as to a Stark fly-by) where the spacecraft, incoming along a purely ballistic trajectory (hyperbola) ignites his propulsion system at $\mathbf r_0=[x_0, y_0]$ to keep a constant acceleration $\alpha$ inertially fixed along the direction $\mathbf i_x$. The spacecraft propelled arc targets a closest planetary distance of $r_m$ and the spacecraft cuts off its propulsion system after a time $T$. It is known that a single velocity increment delivered at the pericenter along the direction of the velocity vector is an efficient way to increase the spacecraft velocity as the final velocity increment obtained at the end of the outbound hyperbola arc results to be much higher than the one delivered. To leverage this effect, often  known as the Oberth effect, we require the fixed constant acceleration direction to be aligned with the spacecraft velocity at the closest passage $r_m$. Such an alignment will be lost as the constant acceleration direction is kept fixed while the spacecraft velocity will be bended by the planet gravity. The fly-by trajectory geometry for this case is shown in Figure <ref>. The use of this manoeuvre as part of an interplanetary trajectory can be studied once an efficient and simple procedure to design it is laid down. Assume $v_\infty$ and $r_m$ as known, and consider the problem of designing a Stark fly-by such that the velocity at the pericenter is $v_m$. Choosing the variable $v_m$ to parametrize the Stark fly-by allows to compute the lattice invariants only once and thus computations that make use of Weierstrass functions are extremely efficient as discussed in Section <ref>.

Start computing the velocity at the pericenter $v_{mK}$ along a Keplerian hyperbola defined by the same entry condition $v_\infty$ and pericenter distance $r_m$:

\frac{v_{mK}^2}{2} = \frac{v_\infty^2}{2} + \frac{\mu}{r_m}

The energy along the propelled arc is obtained from Eq.(<ref>) applied at the $y$ axis crossing (i.e. $x=0$):

\begin{equation}
\mathcal E = \frac{v_m^2}{2} - \frac{\mu}{r_m}
\label{eq:energystarkx0}
\end{equation}

which allows to compute $x_0$ from:

\frac{v_\infty^2}{2}-\alpha x_0 = \mathcal E

Note that $x_0$ will be negative if $v_m > v_{mK}$ as the spacecraft will need to accelerate, while
$x_0$ will be positive if $v_m < v_{mK}$ as the spacecraft will need to decelerate. Compute now the motion invariants $h_\xi$ and $h_\eta$ at $x=0$. In this point $\xi=-\eta=\sqrt{r_m}$ and the motion invariants $h_\xi$ and $h_\eta$ can be written as:

\begin{array}{ll}
h_\xi = -\alpha \frac {r_m^2}{2} - \mathcal E r_m + 2 r_m^2\dot\xi_m^2 \\
h_\eta = \alpha \frac {r_m^2}{2} - \mathcal E r_m + 2 r_m^2\dot\eta_m^2
\end{array}

Consider now the $x=0$ point as the pericenter of the propelled arc (in order to maximize the Oberth effect). From Eq.(<ref>) it is easy to derive $\dot \xi^2 = \frac{v_m^2}{4r_m}$, $\dot \eta^2 = \frac{v_m^2}{4r_m}$ since $\dot x = v_m$ and $\dot r= 0$ in this case. Hence we get the following expressions for the motion invariants :

\begin{equation}
\begin{array}{ll}
h_\xi = -\alpha\frac{r_m^2}{2} + \mu \\
h_\eta = \alpha\frac{r_m^2}{2} + \mu
\end{array}
\label{eq:motinvfb}
\end{equation}

Stark fly-bys with motion invariants $\mathcal E = 5.125\cdot 10^{-2}$, $h_\xi=0.9975$, $h_\eta=1.0025$ and $v_\infty = 0.3$ ($v_{m} = 1.45 > v_{mK}$). The $x_0$ coordinate at the beginning of the manoeuvre (left) is determined, while the $y_0$ coordinate is free to vary.

We may now compute the lattice invariants $g_{2,\xi}$, $g_{3,\xi}$ from Eq.(<ref>) and $g_{2,\eta}$, $g_{3,\eta}$ from Eq.(<ref>). The flyby geometry is, though, still not fully determined as there are infinitely many propelled arcs, parametrized by the other initial condition $y_0$, having the computed motion invariants and satisfying the entry condition on $v_\infty$. Figure <ref> visualized the different arcs parametrized by different $y_0$. The initial value $y_0$ can be selected by forcing at $x=0$ the distance from the origin to be equal to the requested $r_m$. Since at $x=0$, $\xi^2 = r$ this requires solving the equation $\xi^2(\tau^*) = r_m$. The value $\tau^*$ is thus found solving the equation $\xi(\tau^*)=\eta(\tau^*) \to x=0$. This is done using Eq.(<ref>) and Eq.(<ref>) to compute $\xi(\tau^*)$ and $\eta(\tau^*)$ and a simple Newton iteration to get $\tau^*$ (as initial guess $\tau^*=0$ reveals to be good in most cases).

§.§ Pseudo-algorithm to design a Stark fly-by

Assume to know the values $\mu, \alpha, v_\infty, r_m, v_m$. Compute $\mathcal E$ from Eq.(<ref>) and $h_\xi$, $h_\eta$ from Eq.(<ref>). Compute the lattice invariants from Eq.(<ref>) and Eq.(<ref>). Assume a value for $y_0$. At the start of the propelled arc we have:

\begin{array}{l}
x_0 = \frac 1{\alpha}\left( \frac{v_\infty^2}{2} - \mathcal E\right) \\
r_0 = \sqrt{x_0^2 + y_0^2}
\end{array}

Compute the parabolic coordinates at the starting point from:

\begin{array}{ll}
\xi_0 = \sqrt{r_0+x_0} & \dot\xi_0 = - \frac{\sign{(y_0)}}{2r_0} \sqrt{\alpha \xi_0^4 + 2\mathcal E\xi_0^2 + 2h_\xi}\\
\eta_0 = \sign{(y_0)} \sqrt{r_0-x_0} & \dot\eta_0 = - \frac{\sign{(y_0)}}{2r_0} \sqrt{-\alpha \eta_0^4 + 2\mathcal E\eta_0^2 + 2h_\eta}
\end{array}

where we have used Eq.(<ref>) and we have chosen the signs so that the concavity of the represented shape is positive. Find $\tau^*$ using Newton iterations to solve the equation:

\xi(\tau^*)=\eta(\tau^*)

where $\xi$ and $\eta$ are given by Eq.(<ref>) and Eq.(<ref>) respectively. Note that $\tau_{m,\xi}$ and $\tau_{m,\eta}$ need to be computed for each assumed $y_0 \to \eta_0, \xi_0$. Iterate using a Newton method on the assumed $y_0$ so that the following relation is satisfied:

\xi^2(\tau^*) = r_m

Having determined $x_0$ and $y_0$ all the remaining relevant quantities can be easily found solving an initial value problem.

§.§ A numerical example

The amplification factor during a Stark fly-by at Jupiter at different $\Delta V$ levels (left). The net gain is also shown (right). A Rosetta type of spacecraft is assumed with $m = 3000$ kg and $T=10$ [N] resulting in $\alpha=0.00\overline 3$. Incoming conditions of $v_\infty = 3000$ [m/s] are assumed and a closest passage radius of $r_m = 10 R_J$.

A numerical example where the formulas derived above reveal their use is shown. We look into a Jupiter fly-by, where the spacecraft approaches Jupiter sphere of influence with a relative velocity $v_\infty = 3000$ [m/s] and performs a Stark fly-by with closest approach distance $r_m = 10 R_J$, where $R_J = 71492$ [km] is the Jupiter radius. The spacecraft, inspired by the Rosetta spacecraft, has a starting mass $m = 3000$ kg and a maximum thrust capability of $T_{max}=10$ [N] resulting in the possibility to apply a constant $\alpha=0.00\overline 3$. We assume a $\Delta V_{lt} = \alpha \Delta T$ available and we compute the velocity increment $v_\infty^+ - v_\infty^-$ at infinite  resulting from a Stark fly-by. The amplification factor, defined as $(v_\infty^+ - v_\infty^-) / \Delta V_{lt}$ as well as the net gain defined as $(v_\infty^+ - v_\infty^-) - \Delta V_{lt}$ is shown in Figure <ref> for different $\Delta V_{lt}$ and assumed $v_m$ magnitudes.

§ THE EULER'S (TWO-FIXED CENTERS) PROBLEM AND VINTI'S PROBLEM

The same solution method reported above for the constant radial acceleration case and detailed above for the 2D Stark problem, can be used to solve the two-fixed centers, or Euler, problem[7]. While the formulae for the time equation and for the out-of-plane movement turns out to be more complicated, the underlying machinery is still based on the same few relations. In the present contribution we will not go in the details of those formulae, it is here sufficient to know that those formulae exist expressing the explicit solution to the Euler problem in its own pseudo-time (as before the solution with respect to time requires solving numerically a time equation expressed, also, in terms of the Weierstrass functions). We here will only remind that, following Aksenov et al.[20], the solution of the two-fixed centers is linked to Vinti's problem[21], so that as we will briefly outline here, having an explicit solution to Euler's problem entails having an explicit analytical solution for the motion of an artificial satellite around an oblate primary, including the $J_2$ and $J_3$ terms.

Indicating with $\mu_1$ and $\mu_2$ the gravity parameters of the two bodies, we recognize that the system energy is:

\mathcal E = \mathcal T + \mathcal V = \frac{\dot{\mathbf r}^2}2 - \frac{\mu_1}{r_1} - \frac{\mu_2}{r_2}

where, with respect to an inertial system having its origin in the center of mass of the system made by the two attracting bodies alone, we have

\begin{array}{l}
r_1 = \sqrt{x^2+y^2+(z-a_1)^2} \\
r_2 = \sqrt{x^2+y^2+(z-a_2)^2} \\
\end{array}

having indicated with $a_1$ and $a_2$ the distances of the two masses from the system origin. We may then expand the inverse distances appearing in the gravitational potential is series of Legendre polynomials obtaining the expansions:

\begin{array}{l}
\frac 1{r_1} = \frac 1r \sum_{\ell=0}^\infty \left(\frac{a_1}{r}\right)^\ell P_\ell\left(\frac zr\right)\\
\frac 1{r_2} = \frac 1r \sum_{\ell=0}^\infty \left(\frac{a_2}{r}\right)^\ell P_\ell\left(\frac zr\right)\\
\end{array}

and thus proving that the gravitational potential of the two-fixed centers problem has the form:

\begin{equation}
\mathcal V = - \frac{\mu_1+\mu_2}{r}\left(1 + \sum_{\ell=1}^\infty \frac{\gamma_\ell}{r^\ell} P_\ell\left(\frac zr\right)\right)
\end{equation}

\gamma_\ell = \frac{\mu_1 a_1^\ell + \mu_2 a_2^\ell}{\mu_1+\mu_2}

Since our reference system has its origin at the center of mass of the two bodies, we have $\gamma_1 = 0$ and thus we can take the sum from $\ell=2$. Lets compare the expression above to the gravitational potential of an axisymmetric spheroid:

\mathcal V_{sph} = - \frac{\mu}{r}\left(1 - \sum_{\ell=2}^\infty J_\ell \left(\frac{R_e}{r}\right)^\ell P_\ell\left(\frac zr\right)\right)

Assuming $\mu_1$, $\mu_2$, $a_1$ and $a_2$ as complex numbers of the form:

\begin{array}{cc}
\mu_1 = \frac \mu 2 (1+\imath\sigma)& a_1 = c (\sigma + \imath) \\
\mu_2 = \frac \mu 2 (1-\imath\sigma)& a_2 = c (\sigma - \imath)
\end{array}

we are guaranteed that the potential $\mathcal V$ will be real[20]. Writing $J_2 R_e^2 = \gamma_2$ and $J_3 R_e^3 = \gamma_3$ the following conditions are found:

\begin{array}{c}
c^2 (1+\sigma^2) = - J_2 \\
2\sigma c^3  (1+\sigma^2) = -J_3
\end{array}

which allow the first terms of $\mathcal V$ to be equal to the first terms of $\mathcal V_{sph}$. The following term will then be $\gamma_4 = c^4(1+\sigma^2)(1-3\sigma^2)$ and cannot be made equal to $J_4 R_e^4$.

The above developments summarize the original idea from Aksenov et al.[20] who, though, did not develop the explicit solution further stopping at the derivation of the quadratures. The use of the Weierstrassian formalism here introduced allows to find such expressions (as shown in the classic case[7]) and thus to have a fully analytical solution to the problem of an artificial satellite motion around an oblate primary precise up to the $J_3$ term.

§ CONCLUSION

Weierstrass elliptic and related functions express via an elegant formalism solutions to fundamental problems in astrodynamics. These include the constant radial acceleration problem, the Stark problem and the Euler, two-fixed center problem. In all cases the same solution procedure can be applied, resulting in explicit solutions with respect to a Sundmann transformed pseudo-time (or anomaly when possible). The resulting new approach proves to be useful to describe powered fly-bys and the motion of an artificial satellite subject to the gravitational field of an oblate primary including $J_2$ and $J_3$ harmonics (Vinti's problem).

[1]
Biscani, F. and Carloni, S., A first-order secular theory for the
post-Newtonian two-body problem with spin – I. The restricted case, Monthly Notices of the Royal Astronomical Society, Vol. 428, No. 3, Jan.
2013, pp. 2295–2310.

[2]
Gibbons, G. W. and Vyska, M., The application of Weierstrass elliptic
functions to Schwarzschild null geodesics, Classical and Quantum
Gravity, Vol. 29, No. 6, March 2012, pp. 065016.

[3]
Scharf, G., Schwarzschild Geodesics in Terms of Elliptic Functions and
the Related Red Shift, Journal of Modern Physics, Vol. 2, No. 4,
2011, pp. 274–283.

[4]
Hackmann, E., Lämmerzahl, C., Kagramanova, V., and Kunz, J.,
Analytical solution of the geodesic equation in Kerr-(anti-) de
Sitter space-times, Physical Review D, Vol. 81, No. 4, Feb. 2010,
pp. 044020.

[5]
Izzo, D. and Biscani, F., Explicit Solution to the Constant Radial
Acceleration Problem, Journal of Guidance, Control, and Dynamics,
Vol. 38, No. 4, 2014, pp. 733–739.

[6]
Biscani, F. and Izzo, D., The Stark problem in the Weierstrassian
formalism, Monthly Notices of the Royal Astronomical Society,
Vol. 439, No. 1, 2014, pp. 810–822.

[7]
Biscani, F. and Izzo, D., A complete and explicit solution to the
three-dimensional problem of two fixed centres, Monthly Notices of the
Royal Astronomical Society, Vol. 455, No. 4, 2016, pp. 3480–3493.

[8]
Bianchi, L., Lezioni sulla teoria delle funzioni di variabile complessa e
delle funzioni ellittiche, Enrico Spoerri, 1901.

[9]
Byrd, P. F. and Friedman, M. D., Handbook of elliptic integrals for
engineers and scientists, Vol. 67, Springer Science & Business Media,

[10]
Olver, F. W. J., Lozier, D. W., Boisvert, R. F., and Clark, C. W., editors,
NIST Handbook of Mathematical Functions, Cambridge University
Press, New York, NY, 2010, Print companion to [22].

[11]
Whittaker, E. T. and Watson, G. N., Modern analysis, CUP,
Cambridge, 1927.

[12]
Hatten, N. and Russell, R. P., Comparison of three Stark problem
solution techniques for the bounded case, Celestial Mechanics and
Dynamical Astronomy, Vol. 121, No. 1, 2015, pp. 39–60.

[13]
Beth, A., Garnier, P., Toublanc, D., Dandouras, I., and Mazelle, C.,
Theory for planetary exospheres: I. Radiation pressure effect on
dynamical trajectories, arXiv preprint arXiv:1502.06701, 2015.

[14]
Beth, A., Garnier, P., Toublanc, D., Dandouras, I., and Mazelle, C.,
Theory for planetary exospheres: II. Radiation pressure effect on
exospheric density profiles, arXiv preprint arXiv:1503.08122, 2015.

[15]
Biscani, F. and Dario, I., On the implementation of the Weierstrassian
functions, In preparation, 2015.

[16]
Battin, R. H., An introduction to the mathematics and methods of
astrodynamics, AIAA, 1999.

[17]
Lantoine, G. and Russell, R. P., Complete closed-form solutions of the
Stark problem, Celestial Mechanics and Dynamical Astronomy,
Vol. 109, No. 4, 2011, pp. 333–366.

[18]
Tannery, J. and Molk, J., Éléments de la Theorié des
Fonctions Elliptiques, Vol. 4, Gauthier Villars Et Fils Imprimeurs,
Paris, 1893.

[19]
, I. S. and Ryzhik, I. M., Table of Integrals, Series, And Products,
Academic Press, 2007.

[20]
Aksenov, E., Grebenikov, E., and Demin, V., The generalized problem of
motion about two fixed centers and its application to the theory of
artificial earth satellites, Soviet Astronomy, Vol. 7, 1963,
pp. 276.

[21]
Vinti, J. P., New approach in the theory of satellite orbits, Physical Review Letters, Vol. 3, No. 1, 1959, pp. 8.

[22]
NIST Digital Library of Mathematical Functions,
http://dlmf.nist.gov/, Release 1.0.10 of 2015-08-07, Online companion to


[1,2]Juan I. Pisula

1]Institute for Biomedical Informatics, Faculty of Medicine and University
Hospital Cologne, University of Cologne, Germany

2]Center for Molecular Medicine Cologne (CMMC), Faculty of Medicine and
University Hospital Cologne, University of Cologne, Germany

3]Cologne Excellence Cluster on Cellular Stress Responses in Aging- Associated
Diseases (CECAD), University of Cologne, Germany

# Language models are good pathologists: using attention-based sequence
reduction and text-pre-trained transformers for efficient WSI classification

<EMAIL_ADDRESS>Katarzyna Bozek<EMAIL_ADDRESS>[  [  [

###### Abstract

In digital pathology, Whole Slide Image (WSI) analysis is usually formulated
as a Multiple Instance Learning (MIL) problem. Although transformer-based
architectures have been used for WSI classification, these methods require
modifications to adapt them to specific challenges of this type of image data.
Among these challenges is the amount of memory and compute required by deep
transformer models to process long inputs, such as the thousands of image
patches that can compose a WSI at $\times 10$ or $\times 20$ magnification. We
introduce SeqShort, a multi-head attention-based sequence shortening layer to
summarize each WSI in a fixed- and short-sized sequence of instances, that
allows us to reduce the computational costs of self-attention on long
sequences, and to include positional information that is unavailable in other
MIL approaches. Furthermore, we show that WSI classification performance can
be improved when the downstream transformer architecture has been pre-trained
on a large corpus of text data, and only fine-tuning less than 0.1% of its
parameters. We demonstrate the effectiveness of our method in lymph node
metastases classification and cancer subtype classification tasks, without the
need of designing a WSI-specific transformer nor doing in-domain pre-training,
keeping a reduced compute budget and low number of trainable parameters.

##  1 Main

Figure 1: Proposed method. a) After a typical MIL pre-processing step, our
SeqShort layer summarises the amount of patches into a small, oredered
sequence of feature vectors which are then classified by some popular deep
transformer model that was pre-trained on an extensive text corpus. b)
Detailed view of the SeqShort layer, where a set of learnable vectors query
the relevant information in the WSI patches via a Multihead Attention
operation. c) The computational cost of a forward pass of a deep transformer
classifier is considerably reduced when our SeqShort layer is used (measured
with the fvcore library by FAIR).

Tranformers [44] have brought several breakthroughs to the disciplines of
natural language processing (NLP) and computer vision (CV). Their capacity to
link information across sequences of vector embeddings, either visual features
or vectorized words, allowed to capture the structure and meaning necessary
for machine translation [44, 31, 41, 7], question-answering [9, 19, 22, 34],
image classification [11, 23, 49] and segmentation [23, 51], and even multi-
modal tasks such as text-to-image generation [36, 39].

Concurrently in the field of digital pathology, the popularization of Multiple
Instance Learning (MIL) [10, 29] approaches for Whole Slide Image (WSI)
analysis allowed for the fast adoption of transformer models in this domain.
By considering each WSI as a set of feature vectors of smaller tissue patches,
this type of data is a natural input to transformer architectures. However,
although transformer-based, these methods are always adapted to the
idiosyncrasies of MIL and histopathology. Given gigapixel image size, out-of-
the-box Vision Transformers (ViTs) [11] are too memory demanding. Diverse
shapes of WSIs and removal of patches consisting of background, artifacts, pen
marker lines, do not allow for a straightforward implementation of local or
windowed attention [37, 26] . Fixed or learnable positional embeddings are
commonly removed or replaced by alternative new proposals of special
positional encoding [40, 21, 50]. These challenges in WSI analysis do not
allow digital pathology researchers and practitioners to take advantage of
common practices and state-of-art models achieved in CV and NLP. To overcome
these shortcomings, we base our work in two observations.

• The redundancy of information present in full-sequence self-attention
operations can be exploited to reduce the computational cost of large inputs
in deep transformer models. Wang et al. [45] base their Linformer model on the
observation that an attention matrix can be approximated with a matrix of
lower rank. The works of Liu et al. [23] and Dai et al. [6] propose to
construct hierarchical representations instead of maintaining full-length,
token-level resolution. The observations made by Clark et al. [4] about the
importance of the [SEP] token and neighboring tokens have inspired several
methods of local and sparse attention [2, 48, 38, 33]. Comprising thousands of
image patches, a WSI can be transformed to a long sequence of vector
embeddings, and we hypothesise that such findings in the transformer
literature are valid to histopathology data as well, and are necessary to
allow for processing of massive in size WSIs.

• Text pre-trained transformers have been proven successful in non-language
related tasks. Recent works have shown that language models pre-trained on
large unstructured text corpora not only perform strongly in various
downstream NLP tasks, but in several tasks outside of this domain, ranging
from solving math problems [5], to zero-shot action plan generation [15]. We
refer to [15] for an extensive enumeration of such works. In the context of
CV, Ilarco et al. [16] showed that text representations of frozen language
models are predictive of visual representations of their corresponding object.
More recently, Lu et al. [24] demonstrate that pre-trained language models can
perform strongly in image classification, numerical computation, and protein
fold prediction when less than 0.1% of their parameters are fine-tuned. We
adhere to the claims that language modeling pre-training can be leveraged to
perform different, out-of-domain tasks, such as WSI classification.

In this work we bridge the gap between state-of-art transformer models and
digital pathology. To allow for processing of thousands of image patches from
a single WSI, we propose SeqShort, a multi-head attention (MHA) input layer
that reduces long input sequences to a fixed-size short sequence that can be
processed by any well known transformer model. Furthermore, we show that
classification performance is increased when this technique is combined with a
frozen pre-trained language model. This way, we construct a deep, yet
computationally inexpensive model that requires a reduced set of trainable
parameters, and performs well in digital pathology tasks, without the need of
handcrafting yet another transformer variant.

##  2 Results

We showcase our results in compressing the visual information of WSIs with our
sequence reduction technique in order to be classified with standard
transformer models, which were pre-trained on text data in language modeling
tasks. We train several architectures and find that text pre-training improves
classification performance in deep transformer models. Furthermore, we show
that positional information is takent into account by the transformer
classifier when its input is treated as an ordered sequence, instead of an
unordered collection of small image patches as commonly done in other MIL
algorithms.

We then examine how our SeqShort layer works to better understand how visual
information in the WSIs is aggregated. We find that only a small subset of
image patches per WSI is relevant to produce our compressed sequence
representations. In addition, we show that although these representations act
as summaries of the WSIs, a simple extension of the attention rollout
algorithm [1] can trace the output of the transformer classifiers back to each
individual image patch, providing an explanation of the classification
outcome.

###  2.1 WSI classification

We measure the performance of our method on three different classification
tasks: Lymph Node Metastases (LNM) classification (Normal vs Metastases);
Invasive Breast Carcinoma (IBC) subtype classification (Invasive Ductal
Carcinoma vs Invasive Lobular Carcinoma); and Renal Cell Carcinoma (RCC)
subtype classification (Papillary Cell Carcinoma vs Chromophobe Cell Carcinoma
vs Clear Cell Carcinoma). For the LNM classification task we use the dataset
provided by the CAMELYON16 grand challenge (https://camelyon16.grand-
challenge.org/), keeping 10% of the training samples as a validation set, and
evaluating on the grand challenge test set. For the cancer subtyping tasks, we
use WSIs collected from The Cancer Genome Atlas (TCGA)
(https://www.cancer.gov/tcga), and follow the same stratified 10-fold
crossvalidation as [25, 3].

We use 256$\times$256 image patches cropped from the WSIs both at $\times 10$
and $\times 20$ magnification. As a data scarcity ablation, we train the
models using the complete datasets or just 25% of the samples. AUROC is used
as classification performance metric. We compare our method against several
weakly supervised architectures. All networks in the comparison work at a
single magnification at a time, and are agnostic of how their input features
vectors were produced. We use an EfficientNetV2-L [42] pre-trained on ImageNet
[8] as patch-level feature extraction network. In this experiment we use a
frozen RoBERTa-base [22] model as MIL classifier, and only fine-tune its layer
normalization layer. Our LM-based model outperforms the rest of the methods in
CAMELYON16 classification, and shows competitive performance in cancer subtype
classification, without being designed for such digital pathology tasks, and
having been pre-trained with data of a different modality (Table 1).

###  2.2 Language modeling pre-training improves WSI classification

We compare the performance of different frozen text-pre-trained transformers
and a baseline transformer encoder trained from scratch. Given our compute and
memory constrains, the SeqShort layer was required in order to train these
models. All the models under test have 12 layers of 12 attention heads and 768
hidden units, resulting in comparable transformer size. The baseline model,
BERT-base, and RoBERTa-base have identical architecture and only differ in
text-pre-training dataset and language modeling task. This experiment was done
at $\times 20$ magnification, using the IBC task dataset. Except for ALBERT-
base [19], every model outperforms the baseline (Table 2).

###  2.3 The role of positional information

Because of varying WSI shapes and sizes, there are no straightforward
solutions for including positional information in the aggregation of instances
of MIL algorithms, particularly with transformer-based approaches. With our
method, since every WSI is reduced to a fixed length sequence of feature
vectors, it is possible to add to them a sequence of learned positional
embeddings, which is a common practice in typical transformers for CV and NLP
tasks. We demonstrate that such simple operation has a positive effect on
performance with a 0.017 increase in AUROC when training at $\times 10$
magnification and a 0.038 increase in AUROC when training at $\times 20$
magnification (Table 3).

In contrast to the ViT architecture, in which the positional embeddings learn
representations indicative of the 2D location of a patch in the image, in our
method, each WSI is treated as an unordered bag of patches that the SeqShort
layer aggregates into a new sequence. Rather than referring to the original
location of the patches in a WSI, positional information in our model makes
reference to the order of the sequence generated by the SeqShort layer. Our
results show that indeed there is relevant information in this order that
plays a role in classification. This experiment was done on our IBC
classification dataset, using RoBERTa-base as the transformer classifier.

###  2.4 Insights into sequence summarization

Figure 2: WSI summarization inspection. A whole slide image, and attention
heatmaps (A,B,C) produced by three different query vectors in SeqShort are
shown. Although different queries show attentions distributed over a broader
or narrower set of patches, the most important instances remain constant among
the heatmaps. The bottom left plot shows the relative entropy from the
attention distributions of the learned queries to uniform attention (averaged
accross the dataset, and sorted for ease of visualization)

We probe the SeqShort layer to examine how a WSI is summarized. We measure the
relative entropy between the attention distributions produced by the different
learned query vectors in SeqShort and a uniform attention distribution.
Relative entropy values close to zero indicate that such queries pay overall
the same amount of attention to all the input patches, whereas higher relative
entropy values suggest that such queries pay more attention to a reduced
subset of input instances. A plot of these values can be shown in Fig. 2. We
do this measurement with every sample of the IBC dataset test splits at
$\times$20 magnification, and average the results.

Fig. 2 shows an example WSI and its respective attention heatmaps of three
different query vectors with different relative entropy. Even though their
attention distributions are spread over various-sized image areas, high-
attention patches are shared across the heatmaps, indicating their importance
over the rest.

###  2.5 Explanations of classification outcome

Figure 3: Attention rollout heatmaps. Left: original WSIs. Right: their
corresponding attention rollout heatmaps. Although the SeqShort layer reduces
the number of feature vectors that the downstream transformer has to process,
it is still possible to backtrack the contribution of each individual image
patch to the classification decision using the most common transformer
explainability method.

Attention heatmaps from the previous experiment illustrate the workings of the
SeqShort layer of our model: they provide insights into how the patches of a
WSI are taken into account to synthesize the intermediate output of our
method.

To illustrate how the model aggregates the patch representations throughout
its forward pass, we apply attention rollout [1] to generate heatmaps that
provide insights into the overall model decision process. We modify the base
case of the recursive definition of attention rollout to take into
consideration that SeqShort is the first layer of the complete model. Our
modified attention rollout is then defined as:

$\boldsymbol{\tilde{A}_{i}}=\begin{cases}\boldsymbol{A_{i}}\cdot\boldsymbol{\tilde{A}_{i-1}}&\text{if
$i>0$}\\\ \begin{bmatrix}\boldsymbol{0}\\\
\boldsymbol{A_{i}}\end{bmatrix}&\text{if $i=0$},\end{cases}$ (1)

where $\boldsymbol{A_{i}}$ is the attention matrix of layer $i$, and
$\boldsymbol{0}$ is the zero vector in row space, to take into account that
the [CLS] token was not present in the MHA operation of SeqShort. Example
heatmaps are shown in Fig. 3.

##  3 Conclussion

In this work we present a standalone layer for sequence reduction aimed to
overcome common challenges in WSI classification with state-of-art
transformers. Our results show that there is redundant information in the set
of instances of a WSI, and that classification is possible by aggregating this
set into a smaller sequence. We find that the summarized sequences have
intrinsic positional information that can be leveraged by the downstream
classifier. We validate our findings by examining the way our SeqShort layer
aggregates the instances of a WSI.

Our sequence reduction method allows us to use powerful deep transformer
architectures, previously not possible to use in digital pathology, while
requiring only a fraction of their computational cost. Furthermore, we show
that performance is increased when these transformers are pre-trained with
language modeling task on an extensive corpus of text data. Our work supports
the view that the feature representations learned by LMs can be transferred to
other modalities and domains.

##  4 Methods

###  4.1 Sequence Shortening

Existing methods [6, 30, 23] for sequence reduction are not suitable for MIL
WSI problems. Since there is no spatial information about instances in an
unordered bag, concatenating neighboring feature vectors or taking their
strided average is meaningless, as the order of the patches in a bag is
arbitrary. Methods that employ a linear projection for dimensionality
reduction after instance concatenation or sequence reshaping are not
applicable to WSIs either, as they require a fixed and known input shape.

Here we propose using MHA for sequence shortening. Given
$\boldsymbol{X}\in\mathbb{R}^{M\times d}$ the sequence of $M$ $d$-dimensional
feature vectors of non-overlapping WSI tiles, we introduce our SeqShort input
layer that generates a new sequence $\boldsymbol{X_{S}}\in\mathbb{R}^{S\times
h}$ with an MHA layer:

$\begin{split}\boldsymbol{X_{S}}&=\mbox{MHA}(Q=\boldsymbol{Q_{l}},K=\boldsymbol{X},V=\boldsymbol{X})+\boldsymbol{Q_{l}}\\\
&=\mbox{Concat}(head_{1},...,head_{k})\boldsymbol{W^{O}}+\boldsymbol{Q_{l}},\end{split}$
(2)

with

$\begin{split}head_{i}=\mbox{Attention}(Q=\boldsymbol{Q_{l}}\boldsymbol{W_{i}^{Q}},K=\boldsymbol{X}\boldsymbol{W_{i}^{K}},V=\boldsymbol{X}\boldsymbol{W_{i}^{V}}),\\\
\mbox{Attention}(Q,K,V)=\mbox{softmax}(QK^{T}/\sqrt{d_{h}})V,\end{split}$ (3)

where $\boldsymbol{Q_{l}}\in\mathbb{R}^{S\times h}$ is a learnable sequence of
$S$ $h$-dimensional query vectors, the matrices $\boldsymbol{W}$ are learnable
linear projections, $d_{h}$ is a scaling factor commonly set as the layer’s
hidden dimension, and $k$ is the number of attention heads of the layer. Both
$S$ and $h$ are hyperparameters independent of the shape of the original
sequence $\boldsymbol{X}$, and it is $S$ which defines the output sequence
length of the MHA operation in SeqShort.

This MHA operation has a sorting effect: independent of the arrangement of the
patch feature vectors in $\boldsymbol{X}$, the first row of
$\boldsymbol{X_{S}}$ aggregates the instances that the first query vector in
$\boldsymbol{Q_{l}}$ agrees with the most; the second row of
$\boldsymbol{X_{S}}$ aggregates the instances that the second query vector in
$\boldsymbol{Q_{l}}$ agrees with the most, and so on. This enables to
incorporate positional information in our model based on a new interpretation:
instead of thinking of the original arrangement of instances in the WSI 2D
space, we consider the order of the rows of $\boldsymbol{X_{S}}$ as the
available positional information possible to encode.

The resulting time complexity of the MHA operation performed by our input
layer is $O(n)$ because of the fixed-size $\boldsymbol{Q_{l}}$, and since
SeqShort is a single layer, the main compute load lies in the subsequent
deeper transformer model in our pipeline. Although our method does not change
the time complexity of the multi-head self-attention (MHSA) layers of the
transformer itself, by performing sequence reduction, the amount of FLOPs it
requires becomes constant with respect to the original number of WSI patches.
The result is an overall considerable reduction of computational cost (Fig.
1).

###  4.2 Transformer Models

We propose the use of popular NLP transformer architectures for sequence
classification. These models have not been applied before in weakly-supervised
histopathology tasks given their computational cost when dealing with
thousands of instances in a single WSI, whereas in NLP sequences typically
comprise less than 512 tokens. The sequence shortening method that we
introduce above allows us to overcome the computational cost problem.

Inspired by the success of pre-trained language models in different tasks
outside NLP, we propose the use of a frozen, language-modeling pre-trained
transformer as MIL classifier. This is motivated by the hypothesis that the
MHSA layers of a transformer language model learn to capture the
interdependencies among the elements of sequences, independent of the original
data modality or domain. We follow [24] and only fine-tune the layer
normalization layers of the model, reducing the amount of trainable parameters
in our transformer encoder from 85M to 36K (only 0.04% of the total amount).

In our experiments we find that a BERT-base encoder [9] pre-trained with the
masked language modeling task of Robustly optimized BERT pre-training Approach
(RoBERTa) [22] on a corpora of more than 160GB of uncompressed text comprised
by BookCorpus [52], CC-News [27], OpenWebText [13] and Stories [43]. We
discard the vocabulary embeddings lookup table of RoBERTa-base as it is not
needed for weakly supervised image classification.

###  4.3 Complete Pipeline

Our pipeline is illustrated in Figure 1. As a pre-processing step, we extract
non-overlapping tissue tiles of $256\times 256$ pixels from each WSI. Tissue
segmentation is done as in [14] This step is done at $\times$20 and $\times$10
magnification for different experiments. We then generate the instance-level
feature vectors using an EfficientNetV2-L [42] pre-trained on ImageNet [8].

The complete weakly supervised architecture that performs classification on
the bag of instance vectors is composed of the SeqShort layer and a
transformer language model. We set the vector embedding dimension of SeqShort
to $h=768$ (the hidden dimension of the transformers under study), and $k=4$
attention heads. For the lymph node classification task we set the output
length of SeqShort to $S=511$, and for the cancer subtyping tasks, this was
set to $S=256$. A learnable [CLS] token is concatenated to the output of
SeqShort, and added a sequence of learnable positional embeddings. The last
hidden representation of [CLS] is the input of a multilayer perceptron (MLP)
classification head. Altogether, our model comprises a total of 3.3M trainable
parameters.

###  4.4 Implementation and Training

The method was implemented in Python, using PyTorch [32] as deep learning
back-end. The pre-trained weights of EfficientNetV2-L and RoBERTa were
downloaded from Torchvision [28] and HuggingFace [46], respectively. Training
of our models was done with the aid of PyTorch-Lightning [12], on a single
NVIDIA Tesla V100 GPU. The code of this project is available at
https://github.com/bozeklab/lmagp/.

All our models were trained for 200 epochs. For the lymph node classification
task the first 5 epochs were used as learning rate warm-up stage, followed by
one cycle of a cosine schedule, with a maximum learning rate of $1\times
10^{-4}$, and batch size of 16. For the cancer subtyping tasks, the warm-up
stage lasted 10 epochs, followed by two cycles of a cosine schedule, with a
maximum learning rate of $5\times 10^{-5}$, and batch size of 32. Adam [18]
was used as optimization algorithm.

###  4.5 Datasets description

####  4.5.1 Lymph Node Metastases classification

For this task we used the dataset provided by the CAMELYON16 grand challenge
(https://camelyon16.grand-challenge.org/) which comprises 400 Hematoxylin and
Eosin (H&E) stained WSIs of sentinel lymph nodes of breast cancer patients,
scanned by 3DHISTECH and Hamamatsu scanners at $\times 40$ at the Radboud
University Medical Center and the University Medical Center Utrecht,
Netherlands. The grand challenge dataset is divided in a train set of 270 WSIs
(160 normal slides, and 110 slides containing metastases), and a test set of
129 WSIs (80 normal slides, 49 slides containing metastases). In our
experiments, we divided the provided train set in 90%/10% stratified splits
for training and validation, respectively.

####  4.5.2 Invasive Breast Carcinoma subtype classification

We use a subset of 1,038 H&E stained WSIs from the TCGA-BRCA project within
The Cancer Genome Atlas repository (https://www.cancer.gov/tcga). Out of the
1,038 slides, 889 were of patients with Invasive Ductal Carcinoma, and 149
were of patients with Invasive Lobular Carcinoma. We follow the study design
in [25, 3] and evaluate the models using stratified 10-fold crossvalidation on
patient level.

####  4.5.3 Renal Cell Carcinoma subtype classification

We use 918 H&E stained WSIs of Renal Cell Carcinoma cases from the TCGA
repository. Out of these samples, 289 were of Chromophobe Cell Carcinoma
patients, 118 were of Papillary Cell Carcinoma patients, and 498 were of Clear
Cell Carcinoma patients, coming from the TCGA-KICH, TCGA-KIRP and TCGA-KIRC
projects, respectively. We follow the same study design as in the IBC subtype
classification task, and evaluate the models using stratified 10-fold
crossvalidation on patient level.

##  5 Acknowledgments

Both K.B. and J.I.P. were hosted by the Center for Molecular Medicine Cologne
throughout this research. K.B. and J.I.P. were supported by the BMBF program
Junior Group Consortia in Systems Medicine (01ZX1917B) and BMBF program for
Female Junior Researchers in Artificial Intelligence (01IS20054).

## References

  *   * Abnar and Zuidema [2020] Abnar S, Zuidema WH (2020) Quantifying attention flow in transformers. In: Jurafsky D, Chai J, Schluter N, et al (eds) ACL 2020. Association for Computational Linguistics, pp 4190–4197
  * Beltagy et al [2020] Beltagy I, Peters ME, Cohan A (2020) Longformer: The long-document transformer. CoRR abs/2004.05150
  * Chen et al [2022] Chen RJ, Chen C, Li Y, et al (2022) Scaling vision transformers to gigapixel images via hierarchical self-supervised learning. In: CVPR 2022. IEEE, pp 16123–16134
  * Clark et al [2019] Clark K, Khandelwal U, Levy O, et al (2019) What does BERT look at? an analysis of bert’s attention. In: Linzen T, Chrupala G, Belinkov Y, et al (eds) ACL 2019. Association for Computational Linguistics, pp 276–286
  * Cobbe et al [2021] Cobbe K, Kosaraju V, Bavarian M, et al (2021) Training verifiers to solve math word problems. CoRR abs/2110.14168
  * Dai et al [2020] Dai Z, Lai G, Yang Y, et al (2020) Funnel-transformer: Filtering out sequential redundancy for efficient language processing. In: Larochelle H, Ranzato M, Hadsell R, et al (eds) NeurIPS 2020
  * Dehghani et al [2018] Dehghani M, Gouws S, Vinyals O, et al (2018) Universal transformers. 10.48550/ARXIV.1807.03819
  * Deng et al [2009] Deng J, Dong W, Socher R, et al (2009) Imagenet: A large-scale hierarchical image database. In: CVPR 2009, Ieee, pp 248–255
  * Devlin et al [2019] Devlin J, Chang MW, Lee K, et al (2019) BERT: Pre-training of deep bidirectional transformers for language understanding. In: NAACL-HLT. Association for Computational Linguistics, Minneapolis, Minnesota, pp 4171–4186
  * Dietterich et al [1997] Dietterich TG, Lathrop RH, Lozano-Pérez T (1997) Solving the multiple instance problem with axis-parallel rectangles. Artif Intell 89(1-2):31–71
  * Dosovitskiy et al [2021] Dosovitskiy A, Beyer L, Kolesnikov A, et al (2021) An image is worth 16x16 words: Transformers for image recognition at scale. In: ICLR
  * Falcon and The PyTorch Lightning team [2019] Falcon W, The PyTorch Lightning team (2019) PyTorch Lightning. 10.5281/zenodo.3828935, URL https://github.com/Lightning-AI/lightning
  * Gokaslan and Cohen [2019] Gokaslan A, Cohen V (2019) Openwebtext corpus. http://Skylion007.github.io/OpenWebTextCorpus
  * Graham et al [2019] Graham S, Vu QD, Raza SEA, et al (2019) Hover-net: Simultaneous segmentation and classification of nuclei in multi-tissue histology images. Medical Image Analysis 58:101563
  * Huang et al [2022] Huang W, Abbeel P, Pathak D, et al (2022) Language models as zero-shot planners: Extracting actionable knowledge for embodied agents. In: Chaudhuri K, Jegelka S, Song L, et al (eds) ICML 2022, PMLR, vol 162. PMLR, pp 9118–9147
  * Ilharco et al [2021] Ilharco G, Zellers R, Farhadi A, et al (2021) Probing contextual language models for common ground with visual representations. In: Toutanova K, Rumshisky A, Zettlemoyer L, et al (eds) NAACL-HLT 2021. Association for Computational Linguistics, pp 5367–5377
  * Ilse et al [2018] Ilse M, Tomczak JM, Welling M (2018) Attention-based deep multiple instance learning. In: Dy JG, Krause A (eds) ICML 2018, PMLR, vol 80. PMLR, pp 2132–2141
  * Kingma and Ba [2015] Kingma DP, Ba J (2015) Adam: A method for stochastic optimization. In: Bengio Y, LeCun Y (eds) 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, URL http://arxiv.org/abs/1412.6980
  * Lan et al [2020] Lan Z, Chen M, Goodman S, et al (2020) ALBERT: A lite BERT for self-supervised learning of language representations. In: ICLR
  * Li et al [2021a] Li B, Li Y, Eliceiri KW (2021a) Dual-stream multiple instance learning network for whole slide image classification with self-supervised contrastive learning. In: CVPR 2021. Computer Vision Foundation / IEEE, pp 14318–14328
  * Li et al [2021b] Li H, Yang F, Zhao Y, et al (2021b) DT-MIL: deformable transformer for multi-instance learning on histopathological image. In: de Bruijne M, Cattin PC, Cotin S, et al (eds) MICCAI 2021, Part VIII, LNCS, vol 12908. Springer, pp 206–216
  * Liu et al [2019] Liu Y, Ott M, Goyal N, et al (2019) Roberta: A robustly optimized BERT pretraining approach. CoRR abs/1907.11692
  * Liu et al [2021] Liu Z, Lin Y, Cao Y, et al (2021) Swin transformer: Hierarchical vision transformer using shifted windows. In: ICCV. IEEE, pp 9992–10002
  * Lu et al [2022] Lu K, Grover A, Abbeel P, et al (2022) Frozen pretrained transformers as universal computation engines. In: AAAI-IAAI-EAAI 2022. AAAI Press, pp 7628–7636
  * Lu et al [2021] Lu MY, Williamson DF, Chen TY, et al (2021) Data-efficient and weakly supervised computational pathology on whole-slide images. Nature biomedical engineering 5(6):555–570
  * Lv et al [2022] Lv Z, Yan R, Lin Y, et al (2022) Joint region-attention and multi-scale transformer for microsatellite instability detection from whole slide images in gastrointestinal cancer. In: Wang L, Dou Q, Fletcher PT, et al (eds) MICCAI 2022, Part II, Lecture Notes in Computer Science, vol 13432. Springer, pp 293–302
  * Mackenzie et al [2020] Mackenzie JM, Benham R, Petri M, et al (2020) Cc-news-en: A large english news corpus. In: d’Aquin M, Dietze S, Hauff C, et al (eds) CIKM 2020. ACM, pp 3077–3084
  * Marcel and Rodriguez [2010] Marcel S, Rodriguez Y (2010) Torchvision the machine-vision package of torch. In: Bimbo AD, Chang S, Smeulders AWM (eds) Proceedings of the 18th International Conference on Multimedia 2010, Firenze, Italy, October 25-29, 2010\. ACM, pp 1485–1488, 10.1145/1873951.1874254, URL https://doi.org/10.1145/1873951.1874254
  * Maron and Lozano-Pérez [1997] Maron O, Lozano-Pérez T (1997) A framework for multiple-instance learning. In: Jordan MI, Kearns MJ, Solla SA (eds) NeurIPS. The MIT Press, pp 570–576
  * Nawrot et al [2022] Nawrot P, Tworkowski S, Tyrolski M, et al (2022) Hierarchical transformers are more efficient language models. In: Carpuat M, de Marneffe M, Ruíz IVM (eds) NAACL 2022. Association for Computational Linguistics, pp 1559–1571
  * Ott et al [2018] Ott M, Edunov S, Grangier D, et al (2018) Scaling neural machine translation. In: WMT 2018. Association for Computational Linguistics, Brussels, Belgium, pp 1–9
  * Paszke et al [2019] Paszke A, Gross S, Massa F, et al (2019) Pytorch: An imperative style, high-performance deep learning library. In: Wallach HM, Larochelle H, Beygelzimer A, et al (eds) Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pp 8024–8035, URL https://proceedings.neurips.cc/paper/2019/hash/bdbca288fee7f92f2bfa9f7012727740-Abstract.html
  * Qiu et al [2020] Qiu J, Ma H, Levy O, et al (2020) Blockwise self-attention for long document understanding. In: Cohn T, He Y, Liu Y (eds) EMNLP 2020, Findings of ACL, vol EMNLP 2020. Association for Computational Linguistics, pp 2555–2565
  * Radford et al [2018] Radford A, Narasimhan K, Salimans T, et al (2018) Improving language understanding by generative pre-training. Preprint
  * Radford et al [2019] Radford A, Wu J, Child R, et al (2019) Language models are unsupervised multitask learners. OpenAI blog 1(8):9
  * Ramesh et al [2021] Ramesh A, Pavlov M, Goh G, et al (2021) Zero-shot text-to-image generation. In: Meila M, Zhang T (eds) ICML, PMLR, vol 139. PMLR, pp 8821–8831
  * Reisenbüchler et al [2022] Reisenbüchler D, Wagner SJ, Boxberg M, et al (2022) Local attention graph-based transformer for multi-target genetic alteration prediction. In: Wang L, Dou Q, Fletcher PT, et al (eds) MICCAI 2022, Part II, LNCS, vol 13432\. Springer, pp 377–386
  * Roy et al [2021] Roy A, Saffar M, Vaswani A, et al (2021) Efficient content-based sparse attention with routing transformers. Trans Assoc Comput Linguistics 9:53–68
  * Saharia et al [2022] Saharia C, Chan W, Saxena S, et al (2022) Photorealistic text-to-image diffusion models with deep language understanding. CoRR abs/2205.11487
  * Shao et al [2021] Shao Z, Bian H, Chen Y, et al (2021) Transmil: Transformer based correlated multiple instance learning for whole slide image classification. In: Ranzato M, Beygelzimer A, Dauphin YN, et al (eds) NeurIPS 2021, pp 2136–2147
  * So et al [2019] So D, Le Q, Liang C (2019) The evolved transformer. In: ICML 2019, PMLR, pp 5877–5886
  * Tan and Le [2021] Tan M, Le QV (2021) Efficientnetv2: Smaller models and faster training. In: Meila M, Zhang T (eds) ICML 2021, PMLR, vol 139. PMLR, pp 10096–10106
  * Trinh and Le [2018] Trinh TH, Le QV (2018) A simple method for commonsense reasoning. CoRR abs/1806.02847
  * Vaswani et al [2017] Vaswani A, Shazeer N, Parmar N, et al (2017) Attention is all you need. NeurIPS 2017 30
  * Wang et al [2020] Wang S, Li BZ, Khabsa M, et al (2020) Linformer: Self-attention with linear complexity. CoRR abs/2006.04768
  * Wolf et al [2019] Wolf T, Debut L, Sanh V, et al (2019) Huggingface’s transformers: State-of-the-art natural language processing. CoRR abs/1910.03771. URL http://arxiv.org/abs/1910.03771, 1910.03771
  * Yang et al [2019] Yang Z, Dai Z, Yang Y, et al (2019) Xlnet: Generalized autoregressive pretraining for language understanding. NeurIPS 2019 32
  * Zaheer et al [2020] Zaheer M, Guruganesh G, Dubey KA, et al (2020) Big bird: Transformers for longer sequences. In: Larochelle H, Ranzato M, Hadsell R, et al (eds) NeurIPS 2020
  * Zhai et al [2022] Zhai X, Kolesnikov A, Houlsby N, et al (2022) Scaling vision transformers. In: CVPR. IEEE, pp 1204–1213
  * Zhao et al [2022] Zhao Y, Lin Z, Sun K, et al (2022) SETMIL: spatial encoding transformer-based multiple instance learning for pathological image analysis. In: Wang L, Dou Q, Fletcher PT, et al (eds) Medical Image Computing and Computer Assisted Intervention - MICCAI 2022 - 25th International Conference, Singapore, September 18-22, 2022, Proceedings, Part II, LNCS, vol 13432. Springer, pp 66–76, URL https://doi.org/10.1007/978-3-031-16434-7_7
  * Zheng et al [2021] Zheng S, Lu J, Zhao H, et al (2021) Rethinking semantic segmentation from a sequence-to-sequence perspective with transformers. In: CVPR. Computer Vision Foundation / IEEE, pp 6881–6890
  * Zhu et al [2015] Zhu Y, Kiros R, Zemel RS, et al (2015) Aligning books and movies: Towards story-like visual explanations by watching movies and reading books. In: ICCV 2015\. IEEE Computer Society, pp 19–27

Table 1: Performance of different MIL algorithms in the different slide-level classification tasks. Best and the second best classification results are in bold and underlined, respectively. Method | x10 magnification | x20 magnification
---|---|---
25% train set | 100% train set | 25% train set | 100% train set
Lymph Node Metastases classification
ABMIL [17] | 0.501 | 0.664 | 0.516 | 0.616
CLAM [25] | 0.511 | 0.692 | 0.516 | 0.673
DS-MIL [20] | 0.468 | 0.695 | 0.441 | 0.640
TransMIL [40] | 0.529 | 0.629 | 0.470 | 0.723
Ours | 0.627 | 0.772 | 0.642 | 0.865
Invasive Breast Carcinoma subtype classification
ABMIL [17] | 0.542 ± 0.107 | 0.571 ± 0.088 | 0.551 ± 0.103 | 0.554 ± 0.107
CLAM [25] | 0.811 ± 0.055 | 0.850 ± 0.039 | 0.697 ± 0.056 | 0.791 ± 0.082
DS-MIL [20] | 0.779 ± 0.075 | 0.892 ± 0.045 | 0.711 ± 0.084 | 0.819 ± 0.082
TransMIL [40] | 0.864 ± 0.063 | 0.896 ± 0.048 | 0.782 ± 0.094 | 0.856 ± 0.064
Ours | 0.874 ± 0.052 | 0.901 ± 0.049 | 0.765 ± 0.099 | 0.863 ± 0.047
Renal Cell Carcinoma subtype classification
ABMIL [17] | 0.724 ± 0.077 | 0.795 ± 0.040 | 0.697 ± 0.077 | 0.758 ± 0.044
CLAM [25] | 0.965 ± 0.013 | 0.969 ± 0.025 | 0.961 ± 0.013 | 0.974 ± 0.010
DS-MIL [20] | 0.941 ± 0.047 | 0.971 ± 0.001 | 0.926 ± 0.025 | 0.963 ± 0.001
TransMIL [40] | 0.962 ± 0.015 | 0.980 ± 0.001 | 0.971 ± 0.010 | 0.980 ± 0.001
Ours | 0.942 ± 0.019 | 0.974 ± 0.011 | 0.952 ± 0.017 | 0.977 ± 0.013
Table 2: Performance of different Language Models in IBC subtype classification, at $\times 20$ magnification. Baseline model is a transformer encoder of the same size but trained from scratch. Language Model | AUROC
---|---
Baseline | 0.784 ± 0.082
XLNet-base [47] | 0.819 ± 0.090
GPT2-small [35] | 0.827 ± 0.079
BERT-base [9] | 0.849 ± 0.058
ALBERT-base [19] | 0.747 ± 0.118
RoBERTa-base [22] | 0.863 ± 0.047
Table 3: Effect of including positional information on classification performance of IBC subtyping.  Magnification | Pos. emebddings | AUROC
---|---|---
$\times 20$ | No | 0.825 ± 0.052
Yes | 0.863 ± 0.047
$\times 10$ | No | 0.884 ± 0.062
Yes | 0.901 ± 0.049