# Computing Scattering Cross Sections For Spherically Symmetric Potentials

Anil Khachi1

1Department of Physics,
St. Bede’s College, 171002,
Himachal Pradesh, India

Note: This paper is accepted in American Journal of Physics (26/04/2024)

This paper introduces students and instructors to the use of Scilab for
calculating phase shifts from phenomenological potentials in nuclear, atomic,
and molecular scattering, beginning with a Riccati-type nonlinear differential
equation derived from the time-independent Schr$\ddot{\text{o}}$dinger
equation. For spherically symmetric potentials, this equation can be easily
solved using Xcos, a Scilab toolbox. [1][2] Scilab is open source software
that is used for numerical computations and can be freely downloaded for
Windows, Linux and Mac OS.[3] Xcos is a Scilab toolbox dedicated to the
modeling and simulation of hybrid dynamic systems [4]. Xcos provides users
with the ability to model complex dynamical systems using a block diagram
editor (GUI based), offering a step-by-step procedure for learning to solve
complex differential equations interactively.

The normal method for calculating scattering cross sections involves solving
the TISE in the interior ($V\neq 0$) and exterior regions and matching
boundary conditions to find the partial wave amplitudes, $C_{\ell,m}$. [5][6]

The Numerov method can also be used to obtain the phase shift [5], but it
requires knowledge of the wavefunction. In contrast, using the Variable Phase
Approach (VPA), it is not necessary to first calculate radial wave functions
and then to find the phases according to their asymptotic behaviour. The VPA
is a valuable teaching tool because it shows how the scattering phase shift is
affected by the presence of a potential. In Calogero’s words: “It is important
to acquire a physical feeling for the relation between the potential and the
phase shift; this feeling, we believe, is provided by the variable phase
approach”. [7]

The VPA was originally proposed by Morse and Allis [8], and later came to be
known as the phase function method [7]. The VPA has been applied to understand
various nuclear scattering reactions like np and $\alpha-\alpha$ [1][9]. In
this approach the TISE for a particle with wavenumber k and orbital angular
momentum $\ell$ undergoing scattering is transformed into a first order non-
linear Ricatti-type equation.

The TISE for a particle with energy E, reduced mass $\mu$, and orbital angular
momentum $\ell$ undergoing scattering by a potential V(r) is given by

$\frac{\hbar^{2}}{2\mu}\bigg{[}\frac{d^{2}}{dr^{2}}+\bigg{(}k^{2}-\frac{\ell(\ell+1)}{r^{2}}\bigg{)}\bigg{]}u_{\ell}(k,r)=V(r)u_{\ell}(k,r).$
(1)

This linear homogeneous second order differential equation can be transformed
into Riccati type phase function equation [7], containing phase shift
information with initial condition $\delta_{\ell}(k,0)=0$, given by

$\small{\frac{d\delta(k,r)}{dr}=-\frac{U(r)}{k}\big{[}\cos\delta_{\ell}(k,r)\hat{j}_{\ell}(kr)-\sin\delta_{\ell}(k,r)\hat{\eta}_{\ell}(kr)\big{]}^{2}}$
(2)

where $U(r)=V(r)/(\hbar^{2}/2\mu)$, $\hat{j_{\ell}}(kr)$ and
$\hat{\eta_{\ell}}(kr)$ are Riccati-Bessel and Riccati-Neumann functions
respectively.

The dimensionless phase shift $\delta(k)$ is obtained from the phase function
$\delta(k,r)$ through the limiting process:
$\delta(k)=\lim_{r\to\infty}\delta_{\ell}(k,r)$. It should be noted that the
VPA applies only to a restricted class of potentials with a definite behaviour
for $r\rightarrow\infty$. If the potential does not become negligible in the
asymptotic region, then the phase function equation becomes ill-defined. For
example, in interactions involving charged systems such as proton-proton or
$\alpha-\alpha$ interactions, the Coulomb potential does not diminish rapidly
to zero; more detailed analysis has been provided by Laha et al. [10]. For
such cases, a screened Coulomb potential, which behaves like a Coulomb
potential for smaller r and decreases exponentially for large r, such as the
Hultén potential [10], can be utilised for the phase function equation to be
valid [11]. Additionally, it should be noted that obtaining a sufficiently
accurate numerical solution of Eq. 4 requires much more care than finding a
very accurate numerical solution of Eq. 3, as the latter equation is
nonlinear. The VPA also encounters limitations when computing phase shifts for
higher partial waves (typically $\ell>$ 6-8) as the Bessel functions involves
more intricate mathematical expressions.

Figure 1: The upper plots show (a) the attractive Yukawa potential and (b) the
attractive+repulsive Yukawa potential. The lower plots show the resulting
phase shifts for a range of k’ s

Advantages of VPA/PFM over traditional methods of solving TISE are:

  1. 1.

The phase function is sufficient for the complete determination of the wave
function

  2. 2.

At each part of potential we can see how the scattering phase develops in
space

More advantages of VPA have been discussed by Babikov [12]. Figure 6 shows the
phase shift calculation for two spherically symmetric potentials: the
attractive Yukawa potential and a combination of attractive and repulsive
Yukawa interactions. The Xcos code that generated these plots is remarkably
simple and is shown in the supplementary material [13]. One can further obtain
the scattering length a, a parameter that provides an effective size of the
potential and indicates how much the presence of the potential alters the
trajectory at low energies, i.e., near k=0, by simply using the formula
$k\cot\delta_{0}(k,r)\approx-1/a+0.5r_{e}k^{2}$ [14]. More examples along with
their Xcos models are given in the supplementary material [13].

In conclusion we quote Calogero: “VPA is very helpful in most problems of
current physical interest, including those in high-energy physics where a
simple potential model is certainly inadequate to describe the dynamics
completely”.

## Acknowledgements

This paper is dedicated to my father Mr. L. R. Khachi for his love and
affection. I would also like to express my sincere gratitude to the anonymous
reviewers for their invaluable feedback and constructive criticism which have
greatly improved the quality of this paper.

##  Supplementary Material

##  1 Detailed Methodology

The Schr$\ddot{\text{o}}$dinger wave equation for a particle with energy E and
orbital angular momentum $\ell$ undergoing scattering is given by

$\frac{\hbar^{2}}{2\mu}\bigg{[}\frac{d^{2}}{dr^{2}}+\bigg{(}k^{2}-\frac{\ell(\ell+1)}{r^{2}}\bigg{)}\bigg{]}u_{\ell}(k,r)=V(r)u_{\ell}(k,r)$
(3)

where $k=\sqrt{E/(\hbar^{2}/2\mu)}$. The second order differential equation
Eq.3 can been transformed into a pair of first order differential equations:
one is a first order phase function equation i.e $\delta^{\prime}_{\ell}$(k,r)
(Eq.4) and another equation involves both amplitude function $A_{\ell}$(k,r)
and the phase function $\delta_{\ell}$(k,r) called amplitude function
$A^{\prime}_{\ell}$(k,r) [7][12]. The phase function equation given by [7]

$\delta_{\ell}^{\prime}(k,r)=-\frac{2\mu}{\hbar^{2}}\frac{U(r)}{k}\bigg{[}\cos(\delta_{\ell}(k,r))\hat{j}_{\ell}(kr)-\sin(\delta_{\ell}(k,r))\hat{\eta}_{\ell}(kr)\bigg{]}^{2}$
(4)

The phase shift $\delta(k)$ is obtained from the phase function $\delta(k,r)$
through the limiting process:

$\delta(k)=\lim_{r\to\infty}\delta_{\ell}(k,r)$ (5)

In integral form the above equation can be written as

$\delta(k,r)=-\frac{2\mu}{\hbar^{2}}\frac{1}{k}\int_{0}^{r}{U(r)}\bigg{[}\cos(\delta_{\ell}(k,r))\hat{j_{\ell}}(kr)-\sin(\delta_{\ell}(k,r))\hat{\eta_{\ell}}(kr)\bigg{]}^{2}dr$
(6)

Here $\hat{j}_{\ell}$ and $\hat{\eta}_{\ell}$ are Riccati-Bessel and Riccati-
Neumann functions and $U(r)=V(r)/(\hbar^{2}/2\mu)$. Since we have not taken
any physical scattering system in this paper so we have considered
$(\hbar^{2}/2\mu)=1$. Normally $V(r)$ is measured in units of MeV,
$\hbar^{2}/2\mu$ in MeV.fm2 and wavenumber k in fm-1, as a result overall
$\delta(k,r)$ comes out to be dimensionless quantity with unit of radians or
degrees. For $\ell=0$, $\hat{j}_{0}(kr)=\sin(kr)$ and
$\hat{\eta}_{0}(kr)=-\cos(kr)$, hence equation 4 takes following simple form:

$\frac{d\delta_{0}(r)}{dr}=-\frac{V(r)}{k}\sin^{2}(kr+\delta_{0}(r))$ (7)

The value of $\delta_{0}$ at $r=0$ is chosen to be zero, as interaction with
potential is not present, since $u_{0}(r=0)=0$. This Eq. 7 can not be easily
solved using analytical method and hence we take up numerical approach.

For higher partial waves, the Riccati-Bessel and Riccati-Neumann functions
used in VPA/PFM can be easily obtained by using following recurrence formulae:

$\hat{j}_{\ell+1}(kr)=\frac{2\ell+1}{kr}\hat{j_{\ell}}(kr)-{\hat{j}_{\ell-1}}(kr)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(A1)$
$\hat{\eta}_{\ell+1}(kr)=\frac{2\ell+1}{kr}\hat{\eta_{\ell}}(kr)-{\hat{\eta}_{\ell-1}}(kr)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(A2)$

For $\ell$=1 & 2, phase function equation takes following form:

$\delta_{1}^{\prime}(k,r)=-\frac{V(r)}{k}\bigg{[}\frac{\sin(\delta_{1}+(kr))-(kr)\cos(\delta_{1}+(kr))}{(kr)}\bigg{]}^{2}$

$\delta_{2}^{\prime}(k,r)=-\frac{V(r)}{k}\bigg{[}-\sin{\left(\delta_{2}+(kr)\right)}-\frac{3\cos{\left(\delta_{2}+(kr)\right)}}{(kr)}+\frac{3\sin{\left(\delta_{2}+(kr)\right)}}{(kr)^{2}}\bigg{]}^{2}$

Similarly, phase function equations can be obtained for higher partial waves
using equations A1 and A2. Implementation of above VPA equations for all
$\ell$ channels are carried out in Scilab’s Xcos environment [Scilab]

##  2 Results and Discussions

The scattering phase shifts have been obtained for following four spherically
symmetric potentials:

  1. 1.

Kronig-Penney type periodic potential

  2. 2.

Attractive Yukawa

  3. 3.

Repulsive Yukawa

  4. 4.

Combination of attractive and repulsive Yukawa interaction

###  2.0.1 Numerical Implementation of Phase Equation

To solve the phase equation for different $\ell$ channels, one need to specify
the potential U(r) and energies involved through ‘k’. To test the algorithm
and understand the nature of solutions, we have chosen the following three
spherically symmetric potentials:

  * •

Kronig-Penney type periodic potential: The Xcos model for Kronig-Penney type
periodic potential is shown in Fig.2. The potential V(r) is a periodic square
wave as shown in Fig.3, and the corresponding phase shift is shown moving like
stairs above the periodic potential. It can be clearly seen that the phase
shift becomes constant when the potential is constant in accordance with VPA
equation.

Figure 2: Scilab-Xcos scheme for phase shift variation for Kronig-Penney type
periodic potential. Figure 3: Kronig-Penney type periodic potential (black)
and phase shift vs. distance variation (green) for Kronig-Penney type
potential for k=2. Figure 4: Scilab-Xcos scheme for phase shift variation for
attractive ($V(r)=-V_{0}/re^{-r}$) or repulsive potential
($V(r)=+V_{0}/re^{-r}$) w.r.t distance. Figure 5: Phase shift variation with r
for attractive Yukawa potential $V(r)=-10/re^{-r}$ for $\ell=0$.

  * •

Attractive Yukawa: The potential considered is having form
$V(r)=-(10/r)e^{-r}$ and is being solved using VPA for $\ell$ = 0 and 2
partial waves. The Xcos model for attractive/repulsive Yukawa type potential
[16] is given in Fig.4. The $\ell=0$ phase shifts as a function of r are
plotted for values of k = 0.1, 1, 10 and 100 in Fig.5. The SPS reaches a
constant value after the interaction potential goes to zero for values of
$r\geq 2.5fm$. With the increase in k value for $r>1.5$, the SPS decreases,
which is in accordance with the phase equation 4.

Also, with an increase in k value for larger distances, the incident
projectile tends to move away from the attractive core, resulting in a
decrease in SPS with r. Step like behaviour of phase shift in Fig.5 for
smaller k indicates that bound states (here one bound state for r $<$ 1.5 and
no bound state for r $<$ 0.15) are present in accordance with Levinson’s
theorem [15], which states that at k = 0, $\delta=n\pi$, where ‘n’ denotes the
number of bound states. This simply states that as the collision energy tends
to zero (k $\rightarrow$ 0), the phase shift tends to $n\pi$, for bound
states.

For similar attractive potential, with $\ell=2$, no bound state exists because
of the centrifugal barrier, and it can be seen to be pierced by more energetic
particles. The more energetic the scattering particles (for smaller distances
shown in inset), the more they penetrate into the potential region, and
therefore, the larger the scattering phase shift as shown in the inset of Fig.
6. On moving away from the attractive core, there is very weak potential,
almost approaching zero, hence SPS goes down.

Figure 6: Phase shift variation with r for similar attractive Yukawa potential
but for $\ell=2$. Inset shows variation of phase shift for smaller distances.

  * •

Repulsive Yukawa: $V(r)=(20/r)e^{-2r}$ for $\ell=0$: The potential is having
similar form except that it is having shorter range and greater strength. In
case of repulsive potential, the phase shifts are negative which is in
accordance with main VPA equation. In this case, for large, k the particles
face the repulsive core of the potential (r $\approx$ 0.03) and hence phase
gets shifted by large angles. This is shown in inset of Fig. 7. The
accumulation of phase shift begins for smaller value of k and larger value of
r.

Figure 7: $\ell=0$ phase shifts at different k values for Repulsive Yukawa
potential.

  * •

Combination of attractive and repulsive Yukawa:
$V(r)=-(10/r)e^{-r}+(20/r)e^{-2r}$ for $\ell$ = 0, 2, 4 and 6. The Scilab-Xcos
scheme for phase shift variation is shown in Fig.8. Fig.9 depicts some
intricate phase shift plots when both positive and repulsive potential are
present. For $\ell=0$ and k = 5, particles face the repulsive core of the
potential, while in comparison, the k = 1 (less energy) particles, which face
the attractive part of the potential. As a result, the phase shift is negative
for (k = 5, $\ell=0$) and increasing for (k = 1, $\ell=0$). Because the
particles may have encountered the repulsive core at some point in time, a
small portion of the phase shift for (k = 1, $\ell=0$) is negative. We can
observe in the same (Fig.9) that the phase shift changes from negative to
positive for higher values of angular momentum because, as $\ell$ increases,
the particle is forced away from the repulsive core toward the attractive core
of the potential.

Figure 8: Scilab-Xcos scheme for phase shift variation for combination of
attractive and repulsive potential ( $V(r)=-(10/r)e^{-r}+(20/r)e^{-2r}$) i.e.,
w.r.t distance. Figure 9: $\ell=0$ phase shift variation with combination of
attractive and repulsive Yukawa potential $V(r)=-10/re^{-r}+20/re^{-2r}$.

Fig.10 depicts a phase variation within the potential. For a large k value,
i.e., k = 10, the phase shift is negative because of the presence repulsive
wall. For lower k values, i.e., k = 0.02 and 0.05, we have the following
observations:

Figure 10: Variation of $\ell=0$ phase shift with combination of attractive
and repulsive Yukawa potential $V(r)=-10/re^{-r}+20/re^{-2r}$. Variation of
potential is shown in black solid line.

    * –

For r $<$ 0.7 $\rightarrow\delta$=-ve

    * –

For r $>$ 0.7 $\rightarrow\delta$=+ve

In summary, for large values of k, the phase shift gets accumulated for
smaller values of r while for small values of k, the phase shift gets
accumulated for large values of r. All discussed qualitative observations are
in accordance to the main VPA equation. Also for attractive potential
(V(r)$<$0), the phase shift is positive and for repulsive potential
(V(r)$>$0), the phase shift is negative. This property is responsible for the
existance of repulsive core in strong intercations between two protons.

##  3 Conclusion

Xcos is a student-friendly toolbox of Scilab. In this paper, we demonstrate
that a Riccati-type nonlinear differential equation, which is derived from
TISE, can be easily solved using the simple blocks available in Xcos (GUI
environment). In addittion to Kronig-Penney type periodic potential we have
utilized three different Yukawa forms of spherically symmetric potentials: (i)
attractive Yukawa, (ii) repulsive Yukawa, and (iii) attractive and repulsive
Yukawa potentials. Phase shift vs. r has been obtained for the stated
potentials, illustrating how different potentials influence the scattering
outputs. For the first time, an Xcos environment has been created to solve the
Riccati type non-linear differential equation, which students can readily
implement for various phenomenological potentials available in nuclear,
atomic, and molecular fields. The paper can be easily extended to obtain the
variation of phase shift with energy by simply noting the phase shift value
obtained at the asymptotic region (V(r)=0) of the used potential.

## References

  * [1] A. Khachi, L. Kumar, M. R. Ganesh Kumar, and O. S. K. S. Sastri, “Deuteron structure and form factors: Using an inverse potential approach.” Phys. Rev. C 107, 064002 (2023).
  * [2] S. Singh, R. Naghma, J. Kaur and B. Antony, “Calculation of total and ionization cross sections for electron scattering by primary benzene compounds.” J. Chem. Phys., 145, 034309 (2016).
  * [3] Scilab can be downloaded from: http://www.scilab.org.
  * [4] A. K. Verma and R. Verma, Introduction to Xcos: A Scilab Tool for Modeling Dynamical Systems kdp.amazon.com (2020).
  * [5] E. W. Schmid, G. Spitz and W. L$\ddot{o}$sch, “Theoretical Physics on the Personal Computer.” Springer, 189-208 (1990).
  * [6] D. J. Griffiths, Introduction to Quantum Mechanics, Pearson Education (357-362) (1995).
  * [7] F. Calogero. “Variable Phase Approach to Potential Scattering.” Academic, New York (1967).
  * [8] P. M. Morse and W. P. Allis, “The effect of exchange on the scattering of slow electrons from atoms.” Phys. Rev. 44, 269 (1933).
  * [9] A. Khachi, L. Kumar, A. Awasthi and O. S. K. S. Sastri, “ Inverse potentials for all $\ell$ channels of neutron-proton scattering using reference potential approach.” Phys. Scr. 98, 095301 (2023).
  * [10] J. Bhoi and U. Laha, “Hulthen potential models for $\alpha-\alpha$ and $\alpha$-He3 elastic scatterings.” Pramana, 86, 1-6 (2017).
  * [11] Y. Fujiwara and K. Fukukawa. “A Practical Method of Solving Cutoff Coulomb Problems in Momentum Space: Application to the Lippmann-Schwinger Resonating-Group Method and the pd Elastic Scattering.” Prog. Theor. Phys. 128, 301 (2012).
  * [12] V. V. Babikov. “The phase-function method in quantum mechanics.” In: Sov. Phys. Usp. 92, 271 (1967).
  * [13] See Supplemental Material at [url to be inserted by AIPP]
  * [14] G. Darewych and A. E. S. GREEN. “Morse Function and Velocity-Dependent Nucleon-Nucleon Potentials.” Phys. rev. 164, 1324 (1967).
  * [15] N. Levinson. “On the uniqueness of the potential in a Schrodinger equation for a given asymptotic phase”. Danske Videnskab Selskab. Mat. Fys. Medd. 25 (1949).
  * [16] H. Yukawa, S. Sakata, M. Kobayasi, and M. Taketani. “On the Interaction of Elementary Particles. IV.” Progress of Theoretical Physics Supplement 1, 46-71 (1955).


# End-to-End Constrained Optimization Learning: A Survey

James Kotary1 Ferdinando Fioretto1 Pascal Van Hentenryck2&Bryan Wilder3
1Syracuse University
2Georgia Institute of Technology
3Harvard University {jkotary<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper surveys the recent attempts at leveraging machine learning to solve
constrained optimization problems. It focuses on surveying the work on
integrating combinatorial solvers and optimization methods with machine
learning architectures. These approaches hold the promise to develop new
hybrid machine learning and optimization methods to predict fast, approximate,
solutions to combinatorial problems and to enable structural logical
inference. This paper presents a conceptual review of the recent advancements
in this emerging area.

##  1 Introduction

_Constrained optimization_ (CO) has made a profound impact in industrial and
societal applications in numerous fields, including transportation, supply
chains, energy, scheduling, and the allocation of critical resources. The
availability of algorithms to solve CO problems is highly dependent on their
form, and they range from problems that are efficiently and reliably solvable,
to problems that provably have no efficient method for their resolution. Of
particular interest in many fields are _combinatorial optimization_ problems,
which are characterized by discrete state spaces, and whose solutions are
often combinatorial in nature, involving the selection of subsets,
permutations, paths through a network, or other discrete structures to compose
a set of optimal decisions. They are known for their difficulty, and are often
NP-Hard.

Despite their complexity, many CO problems are solved routinely and the AI and
Operational Research communities have devised a wide spectrum of techniques
and algorithms to effectively leverage the problem structure and solve many
hard CO problems instances within a reasonable time and with high accuracy.
While this success has made possible the deployment of CO solutions in many
real-world contexts, the complexity of these problems often prevent them to be
adopted in contexts of repeated (e.g., involving expensive simulations, multi-
year planning studies) or real-time nature, or when they depend in nontrivial
ways on empirical data. However, in many practical cases, one is interested in
solving problem instances sharing similar patterns. Therefore, machine
learning (ML) methods appear to be a natural candidate to aid CO decisions and
have recently gained traction in the nascent area at the intersection between
CO and ML.

Current research areas in the intersection of CO and ML can be categorized
into two main directions: _ML-augmented CO_ and _End-to-End CO learning_. The
former focuses on using ML to aid the decisions performed within an
optimization algorithm used to solve CO problems. The latter involves the
combination of ML and CO techniques to form integrated models which predict
solutions to optimization problems. The survey subdivides this context to into
two main application domains: (1) the fast, approximate prediction of
solutions to CO problems and (2) the integration of data-driven inference with
CO solvers for structure logical inference.

While there exists work surveying ML-augmented CO methods (Bengio et al.,
2020; Lodi and Zarpellon, 2017), the more modern end-to-end CO learning
methods lack of a cohesive and critical analysis. The goal of this survey is
to address this gap and provide a focused overview on the work to-date on end-
to-end CO learning, provide a critical analysis, and pose a set of open
questions and directions.

##  2 Preliminaries: Constrained Optimization

A constrained optimization (CO) problem poses the task of minimizing an
_objective function_ $f:\mathbb{R}^{n}\to\mathbb{R}_{+}$ of one or more
variables $\bm{z}\in\mathbb{R}^{n}$, subject to the condition that a set of
_constraints_ $\mathcal{C}$ are satisfied between the variables:

$\mathcal{O}=\operatorname*{argmin}_{\bm{z}}f(\bm{z})\;\;\text{subject
to}\;\;\bm{z}\in\mathcal{C}.$ (1)

An assignment of values $\bm{z}$ which satisfies $\mathcal{C}$ is called a
_feasible solution_ ; if, additionally $f(\bm{z})\leq f(\bm{w})$ for all
feasible $\bm{w}$, it is called an _optimal solution_.

A well-understood class of optimization problems are _convex_ problems, those
in which the constrained set ${\cal C}$ is a convex set, and the objective $f$
is a convex function. Convex problems have the favorable properties of being
efficiently solvable with strong theoretical guarantees on the existence and
uniqueness of solutions (Boyd et al., 2004).

A particularly common constraint set arising in practical problems takes the
form ${\cal{C}}=\\{\bm{z}\;:\;\bm{A}\bm{z}\leq\bm{b}\\}$, where
$\bm{A}\in\mathbb{R}^{m\times n}$ and $\bm{b}\in\mathbb{R}^{m}$. In this case,
${\cal{C}}$ is a convex set. If the objective $f$ is an affine function, the
problem is referred to as _linear program_ (LP). When a linearly constrained
problem includes a quadratic objective rather than a linear one, the result is
called a _quadratic program_ (QP). If, in addition, some subset of a problem’s
variables are required to take integer values, it is called _mixed integer
program_ (MIP). While LP and QP with convex objectives belong to the class of
convex problems, the introduction of integral constraints
($\bm{x}\in\mathbb{N}^{n}$) results in a much more difficult problem. The
feasible set in MIP consists of distinct points in $\bm{x}\in\mathbb{R}^{n}$,
not only nonconvex but also disjoint, and the resulting problem is, in
general, NP-Hard. Finally, nonlinear programs (NLPs) are optimization problems
where some of the constraints or the objective function are nonlinear. Many
NLPs are nonconvex and can not be efficiently solved (Nocedal and Wright,
2006).

Of particular interest are the _mixed integer linear programs_ (MILPs), linear
programs in which a subset of variables required to take integral values. This
survey is primarily concerned with optimization problems involving linear
constraints, linear or quadratic objective, and either continuous or integral
variables, or a combination thereof.

##### Optimization Solving Methods

A well-developed theory exists for solving convex problems. Methods for
solving LP problems include _simplex_ methods (Dantzig, 1951), _interior
point_ methods (Boyd et al., 2004) and _Augmented Lagrangian_ methods
(Hestenes, 1969; Powell, 1969). Each of these methods has a variant that
applies when the objective function is convex. Convex problems (Boyd et al.,
2004) are in the class of $\bm{P}$, and can be assumed to be reliably and
efficiently solvable in most contexts. For a review on convex optimization and
Lagrangian duality, the interested reader is referred to (Boyd et al., 2004).

MILPs require a fundamentally different approach, as their integrality
constraints put them in the class of NP-Hard problems. _Branch and bound_ is a
framework combining optimization and search, representable by a search tree in
which a LP _relaxation_ of the MILP is formed at each node by dropping
integrality constraints, and efficiently solved using solution methods for
linear programming. Subsequently, a variable $z_{i}$ with (fractional) value
$a_{i}$ in the relaxation is selected for branching to two further nodes. In
the right node, the constraint $z_{i}\geq a_{i}$ is imposed and in the left,
$z_{i}<a_{i}$; bisecting the search space and discarding fractional values in
between the bounds. The solution of each LP relaxtion provides a lower bound
on the MILP’s optimal objective value. When an LP relaxation happens to admit
an integral solution, an upper bound is obtained. The minimal upper bound
found thus far is used, at each node, for pruning.

Finally, _Constraint Programming_ (Rossi et al., 2006) is an additional
effective paradigm adopted to solve industrial-sized MI(L)P and discrete
combinatorial programs.

##  3 Preliminaries: Deep Learning

Supervised deep learning can be viewed as the task of approximating a complex
non-linear mapping from targeted data. Deep Neural Networks (DNNs) are deep
learning architectures composed of a sequence of layers, each typically taking
as inputs the results of the previous layer (LeCun et al., 2015). Feed-forward
neural networks are basic DNNs where the layers are fully connected and the
function connecting the layer is given by $\bm{o}=\pi(\bm{W}\bm{x}+\bm{b}),$
where $\bm{x}\\!\in\\!\mathbb{R}^{n}$ and is the input vector,
$\bm{o}\\!\in\\!\mathbb{R}^{m}$ the output vector,
$\bm{W}\\!\in\\!\mathbb{R}^{m\times n}$ a matrix of weights, and
$\bm{b}\\!\in\\!\mathbb{R}^{m}$ a bias vector. The function $\pi(\cdot)$ is
often non-linear (e.g., a rectified linear unit (ReLU)).

Supervised learning tasks consider datasets
$\bm{\chi}\\!=\\!\\{\bm{x}_{i},\bm{y}_{i}\\}_{i=1}^{N}$ consisting of $N$ data
points with $\bm{x}_{i}\in\mathcal{X}$ being a feature vector and
$\bm{y}_{i}\in\mathcal{Y}$ the associated targets. The goal is to learn a
model ${\cal M}_{\theta}:\mathcal{X}\to\mathcal{Y}$, where $\theta$ is a
vector of real-valued parameters, and whose quality is measured in terms of a
nonnegative, and assumed differentiable, _loss function_
$\mathcal{L}:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}_{+}$. The learning task
minimizes the empirical risk function (ERM):

$\min_{\theta}J({\cal
M}_{\theta},\bm{\chi})=\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}({\cal
M}_{\theta}(\bm{x}_{i}),\bm{y}_{i}).$ (2)

Most of the techniques surveyed in this work use (variants of) DNNs whose
training conforms to the objective above. Other notable classes of deep
learning methods used to solve CO problems are Sequence Models and Graph
Neural Networks (GNNs), which are reviewed below, and Reinforcement Learning
(RL), which differs from supervised learning in not requiring labeled
input/output pairs and concerns with learning a policy that maximizes some
expected reward function. We refer the reader to (Sutton and Barto, 2018) for
an extensive overview of RL.

##### Sequence Models

Recurrent neural networks (RNNs) have proven highly effective for tasks
requiring pattern recognition in sequential data. A basic RNN uses a
_recurrent layer_ to compute a sequence of output arrays
$\bm{y}_{1},\ldots,\bm{y}_{N}$ from a sequence of input arrays
$\bm{x}_{1},\ldots,\bm{x}_{N}$. Each of several hidden units $h_{t}$ depends
on its corresponding input array $\bm{x}_{t}$ along with the previous hidden
unit $h_{t-1}$ in the following manner:

$\bm{h}_{t}=\pi(\bm{W}^{x}\bm{x}_{t}+\bm{W}^{h}\bm{h}_{t-1}+\bm{b})$

where $\pi$ is a nonlinear activation function, $\bm{W}^{x}$,$\bm{W}^{h}$ are
weight matrices, and $\bm{b}$ is an additive bias. $\bm{W}$ and $\bm{b}$ are
learned parameters, applied at each unit of the recurrent layer.

In their basic form, RNNs are known to suffer from vanishing gradient issues:
the propagation of gradients through many hidden units by the chain rule can
lead to loss of information as gradient magnitudes diminish. Modern RNNs use
gated recurrent units (GRU) to form long-short term memory models (LSTM) which
can better preserve relevant information across many hidden layers.

_Sequence-to-Sequence_ models build on recurrent neural networks with the aim
of extending the framework to handle input and output sequences of variable
size. This is accomplished by the introduction of an _endoder-decoder_
architecture in which an encoder LSTM maps an input sequence ${\bf
x}=(\bm{x}_{1},\ldots,\bm{x}_{N})$ to a fixed-length context vector $\bm{c}$,
by way of the hidden encoder states $(e_{1},\ldots,e_{N})$ (typically,
$\bm{c}=e_{N}$). A subsequent decoder LSTM is used to map the context vector
to a sequence of hidden decoder states $(d_{1},\ldots,d_{M})$, which are used
to target a sequence of labeled vectors
$\bm{y}=(\bm{y}_{1},\ldots,\bm{y}_{M})$. That is, some function $g$ is trained
to approximate the likelihood of predicting $\bm{y}_{t}$ given the context
vector $\bm{c}$, the current decoder state and all previous $\bm{y}_{i}$:

$\Pr(\bm{y}_{t}|\\{\bm{y}_{1},\ldots,\bm{y}_{t-1}\\},\bm{c})=g(\bm{y}_{t-1},\bm{d}_{t},\bm{c}),$

where the overall training goal is to model the probability of predicting
sequence $\bf y$. This joint probability is decomposed into a product of
ordered conditionals:

$\Pr({\bf
y})=\Pi_{t=1}^{M}\Pr(\bm{y}_{t}|\\{\bm{y}_{1},\ldots,\bm{y}_{t-1}\\},\bm{c}).$

The primary downside to this approach is the reliance on a single fixed-length
vector to encode the input sequence, resulting in limited expressivity,
particularly when the input sequences observed at test time are of a different
length than those used during training. _Attention_ mechanisms (Bahdanau et
al., 2016) are used to address this shortcoming. They use the same encoder
design, but allow the decoder to adaptively use any subset of the encoder’s
hidden units during decoding. In place of a single context vector, a
_sequence_ of context vectors $\bm{c}_{r}$ (one per decoder unit) is computed
as

$\displaystyle\bm{c}_{r}$
$\displaystyle=\textstyle\sum_{t=1}^{N}\bm{a}_{rt}e_{t}\quad\forall r\leq M$
(3a) $\displaystyle\bm{a}_{rt}$
$\displaystyle=\textstyle\frac{\exp(\bm{u}_{rt})}{\sum_{k=1}^{N}\exp(\bm{u}_{rk})}$
(3b) $\displaystyle\bm{u}_{rt}$
$\displaystyle=\bm{v}^{T}\pi(\bm{W}_{1}e_{t}+\bm{W}_{2}d_{r}),$ (3c)

where $\bm{a}$ above is computed as the softmax of $\bm{u}$, the result of
_alignment_ model (3c) in which $\bm{v}$, $\bm{W}_{1}$, and $\bm{W}_{2}$ are
learnable matrices and $\pi$ is a nonlinear activation function. The values
$\bm{a}_{rt}$ constitute a sequence of _attention_ vectors which measure the
alignment between decoder state $d_{r}$ and encoder state $e_{t}$. Each
attention vector is used as a set of weights to form a linear combination of
encoder states, used as a context vector to inform subsequent predictions.

_Pointer networks_ (Vinyals et al., 2017) are a special use case of attention
mechanisms, used when the intended output of a network is a permutation of its
input. Rather than utilizing the attention vector to assist output sequence
prediction pointer networks treat the attention vector itself as the
prediction. That is,

$\displaystyle\bm{u}_{rt}=\bm{v}^{T}\pi(\bm{W}_{1}e_{t}+\bm{W}_{2}d_{r})\quad\forall
1\leq t\leq N,$
$\displaystyle\Pr(\bm{y}_{t}|\\{\bm{y}_{1},\ldots,\bm{y}_{t-1}\\},\bm{x})=\text{softmax}(\bm{u}_{r}),$

where $\bm{x}$ is again the input sequence and $\bm{y}$ is a permutation of
$\bm{x}$. The attention vector is viewed as a pointer to a particular encoder
state, and can be used in subsequent predictions.

##### Graph Neural Networks

Graph Neural Networks (GNNs) (Scarselli et al., 2008), learn to represent each
node $v$ of a graph with a _state embedding_ $x_{v}$, which depends on
information from its neighboring nodes. The goal of these embeddings is to
provide useful latent information about each node to subsequent components of
a prediction model. Let $x_{v},x_{co[v]},h_{ne[v]},x_{ne[v]}$ represent,
respectively, the features of node $v$, the features of its edges, the state
embeddings, and features of its neighbors. Then consider two functions:

$\displaystyle h_{v}$ $\displaystyle=f(x_{v},x_{co[v]},h_{ne[v]},x_{ne[v]})$
(4a) $\displaystyle o_{v}$ $\displaystyle=g(h_{v},x_{v})$ (4b)

Equivalently to (4a) and (4b),

$\displaystyle H$ $\displaystyle=F(H,X)$ (5a) $\displaystyle O$
$\displaystyle=G(H,X_{N})$ (5b)

where $X$ and $H$ represent the concatentation of all features and all state
embeddings, and $X_{N}$ contains only the node-wise features. When $F$ is a
contractive mapping, the equation has a fixed point solution $H^{*}$ which can
be found by repeated application of (5a). When $F$ and $G$ are interpreted as
feedforward neural networks, their parameters can be trained to produce state
embeddings that are optimized to assist a particular learning task.

Several variants and enhancements to the basic GNN have been proposed; Zhou et
al. (2018) provides a thorough survey of techniques and applications.

Figure 1: Machine Learning and Constrained Optimization.

##  4 Overview of ML and CO

Current research areas in the synthesis of constrained optimization and
machine learning can be categorized into two main directions: _ML-augmented
CO_ , which focuses on using ML to aid the performance of CO problem solvers,
and _End-to-End CO learning (E2E-COL)_ , which focuses on integrating
combinatorial solvers or optimization methods into deep learning
architectures. A diagram illustrating the main research branches within this
area is provided in Figure 1.

The area related with End-to-End CO learning is the focus of this survey and
is concerned with the data-driven prediction of solutions to CO problems. This
paper divides this area into: (1) approaches that develop ML architectures to
predict fast, approximate solutions to predefined CO problems, further
categorized into _learning with constraints_ and _learning on graphs_ , and
(2) approaches that exploit CO solvers as neural network layers for the
purpose of structured logical inference, referred to here as the _Predict-and-
Optimize_ paradigm.

##  5 ML-augmented CO

ML-augmented CO involves the augmentation of already highly-optimized CO
solvers with ML inference models. Techniques in this area draw on both
supervised and RL frameworks to develop more efficient approaches to various
aspects of CO solving for both continuous and discrete combinatorial problems.

In the context of combinatorial optimization, these are broadly categorized
into methods that learn to guide the search decisions in branch and bound
solvers, and methods that guide the application of primal heuristics within
branch and bound. The former include low-cost emulation of expensive branching
rules in mixed integer programming (Khalil et al., 2016; Gasse et al., 2019;
Gupta et al., 2020), prediction of optimal combinations of low-cost variable
scoring rules to derive more powerful ones (Balcan et al., 2018), and learning
to cut (Tang et al., 2020) in cutting plane methods within MILP solvers. The
latter include prediction of the most effective nodes at which to apply primal
heuristics (Khalil et al., 2017), and specification of primal heuristics such
as the optimal choice of variable partitions in large neighborhood search
(Song et al., 2020). The reader is referred to the excellent surveys (Lodi and
Zarpellon, 2017; Bengio et al., 2020) for a thorough account of techniques
developed within ML-augmented combinatorial optimization.

Techniques in this area have also used ML to improve decisions in continuous
CO problems and include learning restart strategies (Gagliolo and Schmidhuber,
2007), learn rules to ignore some optimization variables leveraging the
expected sparsity of the solutions and consequently leading to faster solvers,
(Ndiaye et al., 2017), and learning active constraints (Ng et al., 2018; Prat
and Chatzivasileiadis, 2020) to reduce the size of the problem to be fed into
a CO solver .

##  6 E2E-COL: Predicting CO Solutions

A diverse body of work within the end-to-end CO learning literature has
focused on developing ML architectures to predict fast, approximate solutions
to predefined CO problems _end-to-end_ without the use of CO solvers at the
time of inference, by observing a set of solved instances or execution traces.
These approaches contrasts with those that use ML to augment search-based CO
solvers and configure their subroutines to direct the solver to find solutions
efficiently. This survey categorizes the literature on _predicting CO
solutions_ into (1) methods that incorporate constraints into end-to-end
learning, for the prediction of feasible or near-feasible solutions to both
continuous and discrete constrained optimization problems, and (2) methods
that learn combinatorial solutions on graphs, with the goal of producing
outputs as combinatorial structures from variable-sized inputs. These two
categories, referred to as _learning with constraints_ and _learning CO
solutions_ , respectively, are reviewed next.

###  6.1 Learning with Constraints

The methods below consider datasets
$\bm{\chi}=\\{\bm{x}_{i},\bm{y}_{i}\\}_{i=1}^{N}$ whose inputs $\bm{x}_{i}$
describe some problem instance specification, such as matrix $\bm{A}$ and
vector $\bm{b}$ describing linear constraints in MILPs, and the labels
$\bm{y}_{i}$ describe complete solutions to problem ${\cal O}$ with input
$\bm{x}_{i}$. Notably, each sample may specify a different problem instance
(with different objective function coefficients and constraints).

An early approach to the use of ML for predicting CO problem solutions was
presented by Hopfield and Tank (1985), which used Hopfield Networks (Hopfield
(2018)) with modified energy functions to emulate the objective of a traveling
salesman problem (TSP), and applied Lagrange multipliers to enforce
feasibility to the problem’s constraints. It was shown in Wilson and Pawley
(1988) however, that this approach suffers from several weakness, notably
sensitivity to parameter initialization and hyperparameters. As noted in Bello
et al. (2017) , similar approaches have largely fallen out of favor with the
introduction of practically superior methods.

Frameworks that exploit Lagrangian duality to guide the prediction of
approximate solutions to satisfy the problem’s constraints have found success
in the context of continuous NLPs including energy optimization (Fioretto et
al., 2020b; Velloso and Van Hentenryck, 2020) as well as constrained
prediction on tasks such as transprecision computing and fair classification
(Fioretto et al., 2020a; Tran et al., 2020).

Other end-to-end learning approaches have demonstrated success on the
prediction of solutions to constrained problems by injecting information about
constraints from targeted feasible solutions. Recently, Detassis et al. (2020)
presented an iterative process of using an external solver for discrete or
continuous optimization to adjust targeted solutions to more closely match
model predictions while maintaining feasibility, reducing the degree of
constraint violation in the model predictions in subsequent iterations.

###  6.2 Learning Solutions on Graphs

By contrast to approaches learning solutions to unstructured CO problems, a
variety of methods learn to solve CO cast on graphs. The development of deep
learning architectures, such as sequence models and attention mechanisms, as
well as GNNs, has provided a natural toolset for these tasks.

The survey categorizes these modern approaches broadly based on whether they
rely supervised learning or reinforcement learning. It observes that: (1)
Supervised learning concerns training models that take as input a
representation of a problem instance and output a solution. Target solutions
must be provided for this training task; this is often problematic as it
entails solving, typically, NP-hard problems. (2) Optimization problems have a
native objective function, which in principle can be used in lieu of a loss
function in supervised learning (possibly obviating the need for training
labels) and can serve as a natural reward function in reinforcement learning.

####  6.2.1 Supervised Learning

Vinyals et al. (2017) introduced the _pointer network_ , in which a sequence-
to-sequence model uses an encoder-decoder architecture paired with an
attention mechanism to produce permutations over inputs of variable size.
Noting that certain CO problems provide a compelling use case for their
architecture, the authors test it for solving the Traveling Salesman (TSP),
Delaunay Triangulation, and Convex Hull problem variants. In each, the
solution to a problem instance can be expressed as a single permutation. They
develop a pointer network model to predict approximately optimal solutions by
learning from previously solved instances in a supervised manner, and
demonstrate some ability to generalize over variable-sized problem instances.
For example, in the $2D$ Euclidean TSP, the pointer network’s inputs are the
$2D$ coordinates of each city that must be visited. A predicted permutation
represents a tour over all cities, and each target label is a permutation
representing the pre-computed minimum-length tour over its corresponding input
coordinates. Any city is directly reachable from any other—by a straight line
whose length is equal to the euclidean distance. In general, the approach
introduced in this work applies only to problems whose solutions take the form
of a single permutation and where all permutations are feasible (e.g., the
problem does not feature constraints on which tours are allowed in the TSP).

Despite the subsequent trend toward RL-oriented frameworks, as discussed in
the next section, Nowak et al. (2018) studied a purely supervised learning
method for the general Quadratic Assignment Problem (QAP). The proposed
solution was based on the use of simple graph neural networks trained on
representations of individual problem instances and their targeted solutions.
The TSP is chosen as a test problem, along with two graph matching problems
(all instances of the QAP). Inferences from the model take the form of
approximate permutation matrices, which are converted into feasible tours by a
beam search. The authors note promising accuracy results on small instances,
but the ability for trained models to generalize their performance to larger-
sized instances was not shown.

Joshi et al. (2019) built on this concept by applying a Graph Convolutional
Network model to the $2D$ Euclidean TSP. Aside from the enhanced deep learning
architecture, most aspects of the approach remain similar. However, accurate
results are reported on problem instances much larger than in Nowak et al.
(2018). The authors observe that compared to contemporary reinforcement
learning approaches, their training method is more sample-efficient, but
results in models that do not generalize nearly as well to problem instances
larger than those used in training.

####  6.2.2 Reinforcement Learning

The pointer network architecture of Vinyals et al. (2017) was adopted by Bello
et al. (2017), who proposed to train an approximate TSP solver with
reinforcement learning rather than supervised learning. The transition to RL
was motivated partly by the difficulties associated with obtaining target
solutions that are optimal, and the existence of nonunique optimal solutions
to TSP instances. Rather than generating and targeting precomputed solutions,
the authors present an _actor-critic_ RL framework using expected tour length
$L(\pi|g)$ as the reward signal, where $g$ and $\pi$ represent a graph
(problem instance) and a permutation (a tour over the graph):

$J(\theta|g)=\mathbb{E}_{\pi\sim p_{\theta}(\pi|g)}L(\pi|g).$ (6)

A policy, represented by the underlying pointer network with parameters
$\theta$, is optimized using the well-known REINFORCE rule for computing
policy gradients (Williams, 1992):

$\nabla_{\theta}J(\theta|g)=\mathbb{E}_{\pi\sim
p_{\theta}(\pi|g)}[(L(\pi|g)-b(g)\nabla_{\theta}\log p_{\theta}\;(\pi|g)]$ (7)

The total training objective is then the expectation of $J(\theta|g)$ over a
distribution of graphs. The policy gradient calculation requires a baseline
function $b$ which estimates the expected reward. In this work, an auxiliary
_critic_ network with its own set of parameters is trained to predict the
expected tour length for any graph in supervised fashion, using empirically
observed tours from the most recent policy as training labels. At the time of
inference, two methods are available for producing tours from a trained model:
(1) A set of tours can be sampled from the stochastic policy by varying the
softmax temperature parameter, from which the shortest is selected or (2) an
active search method which modifies the stochastic policy based on solutions
sampled during inference.

Nazari et al. (2018) departed from the RNN encoder scheme used in Bello et al.
(2017) and Vinyals et al. (2017) arguing that when a problem’s defining data
do not adhere to a natural order (e.g. city coordinates in the TSP), it need
not be modeled sequentially. This provides motivation to leave out the encoder
RNN and consider combinatorial problems whose state changes over time,
replacing the role of encoder hidden states with representations of the
problem in each time frame. The authors consider a capacitated vehicle routing
problem (VRP/CVRP) in which demands are associated to points in space and
change over time; a vehicle with limited capacity must satisfy all the demands
by delivering supply from a central depot while minimizing total tour length.
This setting also diverges from those previously discussed in that solutions
no longer take the form of permutations, since an optimal tour may now include
several trips to the supply depot or a demand sink. The static (location
coordinates) and dynamic (demand value) data defining the problem are
concatenated into a single embedding vector $\bm{x}_{t}$ for each time frame
$t$. In each decoder step, the attention layer computes a context vector
$c_{t}$ based only on $\bm{x}_{t}$ and the previous decoder state $d_{t}$. A
final prediction at time $t$ indicates the next visited located based on
$\bm{x}_{t}$ and $c_{t}$, and the embeddings $\bm{x}_{t}$ can then be updated
based on the resulting state-changes. The training consists of an actor-critic
method similar to that of Bello et al. (2017).

Several improvements of these approaches exists. Notably, Kool et al. (2019)
developed a general reinforcement learning framework based on a Graph
Attention Network architecture (Veličković et al., 2017) and trained with the
REINFORCE policy gradients, and present significant improvements in accuracy
on 2D euclidean TSP over Vinyals et al. (2017), Bello et al. (2017), Dai et
al. (2018), and Nowak et al. (2018).

Differently from previous approaches, Ma et al. (2019) focused on CO instances
with hard constraints and consider TSP variants in which certain tours are
considered infeasible. A mutil-level RL framework is proposed in which each
layer of a hierarchy learns a different policy, from which actions can be
sampled. Lower layers learn policies that result in feasible solutions by
using reward functions that bias feasibility, and provide latent variables to
higher layers which learn to optimize a given objective while sampling actions
from lower layers. Layer - wise policy optimization uses the REINFORCE
algorithm, training each layer in turn before fixing its weights and sampling
actions to train the next layer. Each policy is represented by a graph neural
network to produce context embeddings that capture inter-node relationships
from node-wise inputs, and an attention-based decoder to predict distributions
over candidate actions as in prior works. When the direct output $\bm{u}_{i}$
from layer $i$ is predicted at a particular decoding step, it is combined with
the output of the lower layer $i-1$ to ensure mutual compatibility in the
policy distribution

$\bm{p}_{i}=\textrm{softmax}(\bm{u}_{i}+\alpha\bm{u}_{i-1})$

where $\alpha$ is a trainable parameter.

Finally, Dai et al. (2018) adopted a different RL approach based on a greedy
heuristic framework which determines approximate solutions by the sequential
selection of graph nodes. The selection policy is learned using a combination
of Q-learning and fitted Q-iteration, by optimizing the optimization objective
directly. The graph embedding network structure2vec is used to represent the
policy of the learned greedy algorithm. It is used to produce feature
embeddings for each node of the graph, which are updated between actions and
capture relationships between neighboring nodes, given the current state of a
partial solution. It is noted that this framework excels in combinatorial
problems that have a true graphical structure, as opposed to most previous
previously studied applications (primarily the TSP and VRP variants) in the
underlying ’graph’ is typically fully connected. The Minimum Vertex Cover, Max
Cut, and TSP are used as test test problems, and a strong ability to
generalize over various problem sizes is advertised in comparison to prior
approaches.

Figure 2: Predict-and-optimize framework; gradients of a solver output
(solution) must be computed with respect to its input (problem parameters) in
order to maximize empirical model performance.

##  7 E2E-COL: Predict-and-Optimize

A burgeoning topic in the intersection of ML and CO is the fusion of
prediction (ML) and decision (CO) models, in which decision models are
represented by partially defined optimization problems, whose specification is
completed by parameters that are predicted from data. The resulting composite
models employ constrained optimization as a neural network layer and are
trained _end-to-end_ , based on the optimality of their decisions. This
setting is altogether different in motivation to those previously discussed,
in which the goal was to solve predefined CO instances with increased
efficiency. The goal here is the synthesis of predictive and prescriptive
techniques to create ML systems that learn to make decisions based on
empirical data.

The following constrained optimization problem is posed, in which the
objective function $f_{\bm{y}}$ and feasible region ${\cal C}_{\bm{y}}$ depend
on a parameter vector $\bm{y}$:

$\mathcal{O}(\bm{y})=\operatorname*{argmin}_{\bm{z}}f_{\bm{y}}(\bm{z})\;\;\text{subject
to}\;\;\bm{z}\in\mathcal{C}_{\bm{y}}.$ (8)

The goal here is to use supervised learning to predict $\hat{\bm{y}}$ the
unspecified parameters from empirical data. The learning task is performed so
that the optimal solution ${\cal O}(\hat{\bm{y}})$ best matches a targeted
optimal solution ${\cal O}(\bm{y})$, relative to some appropriately chosen
loss function. The empirical data in this setting is defined abstractly as
belonging to a dataset $\chi$, and can represent any empirical observations
correlating with targeted solutions to (8) for some $\bm{y}$. See Figure 2 for
an illustration.

This framework aims to improve on simpler two-stage approaches, in which a
conventional loss function (e.g. MSE) is used to target labeled parameter
vectors $\bm{y}$ that are provided in advance, before solving the associated
CO problem to obtain a decision. Such approaches are suboptimal in the sense
that predictions of $\bm{y}$ do not take into account the accuracy of the
resulting solution ${\cal O}(\bm{y})$ during training.

We note that there are two ways to view the predictions that result from these
integrated models. If $\hat{\bm{y}}$ is viewed as the prediction, then the
calculation of ${\cal O}(\hat{\bm{y}})$ is absorbed into the loss function
$\mathcal{L}(\hat{\bm{y}},\bm{y})$, which targets the provided parameter
vectors. Otherwise, the loss function $\mathcal{L}({\cal
O}(\hat{\bm{y}}),{\cal O}(\bm{y}))$ is considered separately from the decision
model and aims to match computed optimal solutions to targeted ones. One
advantage sought in either case is the opportunity to minimize during training
the ultimate error in the computed optimal objective values
$f_{\hat{\bm{y}}}({\cal O}(\hat{\bm{y}}))$, relative to those of the target
data. This notion of training loss is known as _regret:_

$\textsl{regret}(\hat{\bm{y}},\bm{y})=f_{\hat{\bm{y}}}({\cal
O}(\hat{\bm{y}}))-f_{\bm{y}}({\cal O}(\bm{y})).$

_Otherwise the optimal solution ${\cal O}(\bm{y})$ is targeted and one can use
$\textsl{regret}({\cal O}(\hat{\bm{y}}),{\cal O}(\bm{y}))$, regardless of
whether the corresponding $\bm{y}$ is available._ Depending on the techniques
used, it may be possible to minimize the regret without access to ground-truth
solutions, as in Wilder et al. (2019), since the targeted solutions ${\cal
O}(\bm{y})$ contribute only a constant term to the overall loss. It is worth
mentioning that because the optimization problem in (8) is viewed as a
function, the existence of nonunique optimal solutions is typically not
considered. The implication then is that ${\cal O}(\bm{y})$ is directly
determined by $\bm{y}$.

Training these end-to-end models involves the introduction of external CO
solvers into the training loop of a ML model, often a DNN. Note that
combinatorial problems with discrete state spaces do not offer useful
gradients; viewed as a function, the _argmin_ of a discrete problem is
piecewise constant. _The challenge of forming useful approximations to
$\frac{\partial\mathcal{L}}{\partial y}$ is central in this context and must
be addressed in order to perform backpropagation_. It may be approximated
directly, but a more prevalent strategy is to model $\frac{\partial{\cal
O}(y)}{\partial y}$ and $\frac{\partial\mathcal{L}}{\partial{\cal O}}$
separately, in which case the challenge is to compute the former term by
_differentiation through argmin_. Figure 2 shows the role of this gradient
calculation in context.

###  7.1 Quadratic Programming

One catalyst for the development of this topic was the introduction of
_differentiable optimization layers_ , beginning with Amos and Kolter (2017)
which introduced a GPU-ready QP solver that offers exact gradients for
backpropagation by differentiating the KKT optimality conditions of a
quadratic program at the time of solving, and utilizing information from the
forward pass to solve a linear system for incoming gradients, once the
outgoing gradients are known. Subsequently, Donti et al. (2019) proposed a
_predict-and-optimize_ model in which QPs with stochastic constraints were
integrated in-the-loop to provide accurate solutions to inventory and power
generator scheduling problems specified by empirical data.

###  7.2 Linear Programming

Concurrent with Donti et al. (2019), an alternative methodology for end-to-end
learning with decision models, called _smart predict-and-optimize_ (SPO), was
introduced by Elmachtoub and Grigas (2020), which focused on prediction with
optimization of constrained problems with linear objectives, in which the cost
vector is predicted from data and the feasible region $\mathcal{C}$ is
invariant to the parameter $\bm{y}$:

$\mathcal{O}(\bm{y})=\operatorname*{argmin}_{\bm{z}}\bm{y}^{T}\bm{z}\;\;\text{subject
to}\;\;\bm{z}\in\mathcal{C}.$ (9)

The target data in this work are the true cost vectors $\bm{y}$, and an
inexact subgradient calculation is used for the backpropagation of regret loss
$\mathcal{L}(\hat{\bm{y}},\bm{y})=\hat{y}^{T}({\cal O}(\hat{\bm{y}})-{\cal
O}(\bm{y}))$ on the decision task, by first defining a convex surrogate upper
bound on regret called the _SPO+loss_ , for which it is shown that ${\cal
O}(\bm{y})-{\cal O}(2\hat{\bm{y}}-\bm{y})$ is a useful subgradient. Since this
work is limited to the development of surrogate loss functions on regret from
the optimization task, it does not apply to learning tasks in which the full
solution to an optimization problem is targeted. The paper includes a
discussion justifying the method’s use on problems with discrete constraints
in $\mathcal{C}$, as in combinatorial optimization, but experimental results
are not provided on that topic. It is, however, demonstrated that the approach
succeeds in a case where $\mathcal{C}$ is convex but nonlinear.

Wilder et al. (2019) introduced an alternative framework to predict-and-
optimize linear programming problems, based on exact differentiation of a
smoothed surrogate model. While LPs are special cases of QPs, the gradient
calculation of Amos and Kolter (2017) does not directly apply due to
singularity of the KKT conditions when the objective function is purely
linear. This is addressed by introducing a small quadratic regularization term
to the LP objective $f_{\bm{y}}(\bm{z})=\bm{y}^{T}\bm{z}$ so that the problem
in (8) becomes

$\mathcal{O}(\bm{y})=\operatorname*{argmin}_{\bm{z}}\bm{y}^{T}\bm{z}+\epsilon\|\bm{z}\|^{2}\;\;\text{subject
to}\;\;\bm{A}\bm{z}\leq\bm{b}.$ (10)

The resulting problems approximate the desired LP, but have unique solutions
that vary smoothly as a function of their parameters, allowing for accurate
backpropagation of the result. The integrated model is trained to minimize the
expected optimal objective value across all training samples, which is
equivalent to minimizing the regret loss but without the need for a target
dataset. This work demonstrated success on problems such as the knapsack
(using LP relaxation) and bipartite matching problems where a cost vector is
predicted from empirical data (e.g., historical cost data for knapsack items),
and is shown to outperform two-stage models which lack integration of the LP
problem. We note that although the differentiable QP solving framework of Amos
and Kolter (2017) is capable of handling differentiation with respect to any
objective or constraint coefficient, this work only report results on tasks in
which the _cost_ vector is parameterized within the learning architecture, and
constraints are held constant across each sample. _This limitation is common
to all of the works described below, as well._

Next, Mandi and Guns (2020) introduced an altogether different approach to
obtaining approximate gradients for the argmin of a linear program. An
interior point method is used to compute the solution of a homogeneous self-
dual embedding with early stopping, and the method’s log-barrier term is
recovered and used to solve for gradients in the backward pass. Equivalently,
this can be viewed as the introduction of a log-barrier regularization term,
by analogy to the QP-based approach of Wilder et al. (2019):

$\mathcal{O}(\bm{y})=\operatorname*{argmin}_{\bm{z}}\bm{y}^{T}\bm{z}+\lambda\left(\sum_{i}\ln(z_{i})\right)\;\;\text{subject
to}\;\;\bm{A}\bm{z}\leq\bm{b}.$

Further, the method’s performance on end-to-end learning tasks is evaluated
against the QP approach of Wilder et al. (2019) on LP-based predict-and-
optimize tasks, citing stronger accuracy results on energy scheduling and
knapsack problems with costs predicted from data.

Berthet et al. (2020) detailed an approach based on stochastic perturbation to
differentiate the output of linear programs with respect to their cost
vectors. The output space of the LP problem is smoothed by adding low-
amplitude random noise to the cost vector and computing the expectation of the
resulting solution in each forward pass. This can be done in Monte Carlo
fashion and in parallel across the noise samples. The gradient calculation
views the solver as a black box in this approach, and does not require the
explicit solving of LP for operations that can be mathematically formulated as
LP, but are simple to perform (e.g., sorting and ranking). Results include a
replication of the shortest path experiments presented in (Pogančić et al.,
2019), in which a model integrated with convolutional neural networks is used
to approximate the shortest path through stages in a computer game, solely
from images.

###  7.3 Combinatorial Optimization

Ferber et al. (2020) extended the approach of Wilder et al. (2019) to
integrate MILP within the end-to-end training loop, with the aim of utilizing
more expressive NP-Hard combinatorial problems with parameters predicted from
data. This is done by reducing the MILP with integer constraints to a LP by a
method of cutting planes. In the ideal case, the LP that results from the
addition of cutting planes has the same optimal solution as its mixed-integer
parent. Exact gradients can then be computed for its regularized QP
approximation as in Wilder et al. (2019). Although the LP approximation to
MILP improves with solving time, practical concerns arise when the MILP
problem cannot be solved to completion. Each instance of the NP-Hard problem
must be solved in each forward pass of the training loop, which poses clear
obstructions in terms of runtime. A disadvantage of the approach is that
cutting-plane methods are generally considered to be inferior in efficiency to
staple methods like branch and bound. Improved results were obtained on
portfolio optimization and diverse bipartite matching problems, when compared
to LP-relaxation models following the approach of Wilder et al. (2019).

Mandi et al. (2020) investigates the application of the same SPO approach to
NP-Hard combinatorial problems. Primary among the challenges faced in this
context is, as in (Ferber et al., 2020), the computational cost of solving
hard problems within every iteration of the training. The authors found that
continuous relaxations of common MILP problems (e.g. knapsack) often offer
subgradients of comaparable quality to the full mixed-integer problem with
respect to the SPO loss, so that training end-to-end systems with hard CO
problems can be simplified in such cases by replacing the full problem with an
efficiently solvable relaxation, an approach termed _SPO-relax_. The authors
put continuous relaxations into the broader category of “weaker oracles” for
the CO solver, which also includes approximation algorithms (e.g. greedy
approximation for knapsack). The main results showed that SPO-relax achieves
accuracy competitive with the full SPO approach but with shorter training
times on a handful of discrete problems. The SPO-relax approach was compared
also against the formulation of Wilder et al. (2019) on equivalent
relaxations, but no clear winner was determined.

Pogančić et al. (2019) introduced a new idea for approximating gradients over
discrete optimization problems for end-to-end training, which relies on
viewing a discrete optimization problem as a function of its defining
parameters (in this context coming from previous layers), whose range is
piecewise constant. The only requirement is that the objective be linear. The
gradient calculation combines the outputs of two calls to an optimization
solver; one in the forward pass on initial parameters $\bm{y}$, and one in the
backward pass on perturbed parameters $\bar{\bm{y}}$. The results are used to
construct a piecewise _linear_ function which approximates the original
solver’s output space, but has readily available gradients. Because the
gradient calculation is agnostic to the implementation of the solver, it is
termed ”black-box differentiation”. As such, input parameters to the solver do
not correspond explicitly to coefficients in the underlying optimization
problem. Results on end-to-end learning for the shortest path problem, TSP and
min-cost perfect matching are shown. In each case, the discrete problem’s
specification is implicitly defined in terms of images, which are used to
predict parameters of the appropriate discrete problem through deep networks.
The optimal solutions coming from blackbox solvers are expressed as binary
indicator matrices in each case and matched to targeted optimal solutions
using a Hamming distance loss function.

Finally, Wang et al. (2019) presented a differentiable solver for the MAXSAT
problem, another problem form capable of representing NP-Hard combinatorial
problems. Approximate gradients are formed by first approximating the MAXSAT
problem as a related semidefinite program (SPD), then differentiating its
solution analytically during a specialized coordinate descent method (Wang et
al., 2018) which solves the SDP. The successful integration of MAXSAT into
deep learning is demonstrated with a model trained to solve sudoku puzzles
represented only by handwritten images.

##  8 Challenges and Research Directions

The current state of the art in integrating combinatorial optimization with
end-to-end machine learning shows promise on challenging tasks for which there
was previously no viable approach. Further, it has been demonstrated that a
variety of non-overlapping approaches can be effective. Despite these
encouraging results, a number of challenges remain that must be addressed to
allow an integration that lives up to its full potential. (1) Despite the
variety of approaches, the success of the predict-and-optimize paradigm has
been demonstrated on a relatively limited set of optimization problems and a
majority of reported experimental results focus on linear programming
formulations. Challenges posed by the parametrization of constraints stand in
the way of broader applications, but have not been yet been addressed. (2)
Issues associated with the runtime of combinatorial solvers in-the-loop still
make some potential applications impractical. (3) Additionally, despite being
possible in theory, the role of the CO model in-the-loop has not been
generalized successfully beyond being applied as the final layer of a deep
model. The use of additional layers beyond the solver, or even compositions of
CO solving layers, could potentially lead to new applications if the practical
challenges were to be overcome. (4) In predicting solutions to CO problems,
the current methods cannot reliably guarantee the problem constraints to be
satisfied. This critical shortcoming may be addressed by integrating ML
approaches with methods from the robust optimization literature or by
developing ad-hoc layers to project the predictions onto the feasible space.
(5) While it has been observed in limited contexts Demirović et al. (2019)
that predict-and-optimize frameworks based on optimization layers are
competitive only when the underlying constrained problem is convex, this area
still lacks theoretical results providing guarantees on their viability or
performance. Finally, (6) the need for uniform benchmark experiments and
systematic comparisons between each predict-and-optimization framework is
apparent. Demirović et al. (2019) provided a study comparing the approaches of
Wilder et al. (2019) and Elmachtoub and Grigas (2020), along with problem-
specific approaches, on knapsack problems but did not conclude strongly as to
which method should be preferred. Further, Mandi et al. (2020) reports that
for knapsack problems, SPO performs comparably on the knapsack problem whether
the LP relaxation or the full problem is used, but does not show that this
effect generalize to other NP-Hard problems. This signals a need for
comprehensive studies that test performance on a variety of hard CO problems.

Although the approaches surveyed are still in an early stage of their
development, and have been adopted mainly for academic purposes, we strongly
believe that the integration between combinatorial optimization and machine
learning is a promising direction for the development of new, transformative,
tools in combinatorial optimization and learning.

## Acknowledgments

This research is partially supported by NSF grant 2007164. Its views and
conclusions are those of the authors only and should not be interpreted as
representing the official policies, either expressed or implied, of the
sponsoring organizations, agencies, or the U.S. government.

## References

  * Amos and Kolter [2017] Brandon Amos and J Zico Kolter.  Optnet: Differentiable optimization as a layer in neural networks.  In International Conference on Machine Learning, pages 136–145. PMLR, 2017.
  * Bahdanau et al. [2016] Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate, 2016.
  * Balcan et al. [2018] Maria-Florina Balcan, Travis Dick, Tuomas Sandholm, and Ellen Vitercik.  Learning to branch.  arXiv:1803.10150, 2018.
  * Bello et al. [2017] Irwan Bello, Hieu Pham, Quoc V. Le, Mohammad Norouzi, and Samy Bengio.  Neural combinatorial optimization with reinforcement learning.  arXiv:1611.09940, 2017.
  * Bengio et al. [2020] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost.  Machine learning for combinatorial optimization: a methodological tour d’horizon.  European Journal of Operational Research, 2020.
  * Berthet et al. [2020] Quentin Berthet, Mathieu Blondel, Olivier Teboul, Marco Cuturi, Jean-Philippe Vert, and Francis Bach.  Learning with differentiable perturbed optimizers.  arXiv:2002.08676, 2020.
  * Boyd et al. [2004] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * Dai et al. [2018] Hanjun Dai, Elias B. Khalil, Yuyu Zhang, Bistra Dilkina, and Le Song.  Learning combinatorial optimization algorithms over graphs.  arXiv:1704.01665, 2018.
  * Dantzig [1951] George B Dantzig.  Maximization of a linear function of variables subject to linear inequalities.  Activity analysis of production and allocation, 13:339–347, 1951\.
  * Demirović et al. [2019] Emir Demirović, Peter J Stuckey, James Bailey, Jeffrey Chan, Chris Leckie, Kotagiri Ramamohanarao, and Tias Guns.  An investigation into prediction+ optimisation for the knapsack problem.  In International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research, pages 241–257. Springer, 2019.
  * Detassis et al. [2020] Fabrizio Detassis, Michele Lombardi, and Michela Milano.  Teaching the old dog new tricks: Supervised learning with constraints.  arXiv:2002.10766, 2020.
  * Donti et al. [2019] Priya L. Donti, Brandon Amos, and J. Zico Kolter.  Task-based end-to-end model learning in stochastic optimization.  arXiv:1703.04529, 2019.
  * Elmachtoub and Grigas [2020] Adam N. Elmachtoub and Paul Grigas.  Smart ”predict, then optimize”.  arXiv:1710.08005, 2020.
  * Ferber et al. [2020] Aaron Ferber, Bryan Wilder, Bistra Dilkina, and Milind Tambe.  Mipaal: Mixed integer program as a layer.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2020.
  * Fioretto et al. [2020a] Ferdinando Fioretto, Pascal Van Hentenryck, Terrence W. K. Mak, Cuong Tran, Federico Baldo, and Michele Lombardi.  Lagrangian Duality for Constrained Deep Learning.  In European Conference on Machine Learning and Knowledge Discovery in Databases ECML PKDD, volume 12461 of Lecture Notes in Computer Science, pages 118–135. Springer, 2020.
  * Fioretto et al. [2020b] Ferdinando Fioretto, Terrence W.K. Mak, and Pascal Van Hentenryck.  Predicting ac optimal power flows: Combining deep learning and lagrangian dual methods.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2020.
  * Gagliolo and Schmidhuber [2007] Matteo Gagliolo and Jürgen Schmidhuber.  Learning restart strategies.  In Proceedings of the International Joint Conference on Artificial Intelligence (IJCAI), pages 792–797, 2007.
  * Gasse et al. [2019] Maxime Gasse, Didier Chételat, Nicola Ferroni, Laurent Charlin, and Andrea Lodi.  Exact combinatorial optimization with graph convolutional neural networks.  arXiv:1906.01629, 2019.
  * Gupta et al. [2020] Prateek Gupta, Maxime Gasse, Elias B. Khalil, M. Pawan Kumar, Andrea Lodi, and Yoshua Bengio.  Hybrid models for learning to branch.  arXiv:2006.15212, 2020.
  * Hestenes [1969] Magnus R Hestenes.  Multiplier and gradient methods.  Journal of optimization theory and applications, 4(5):303–320, 1969\.
  * Hopfield and Tank [1985] John Hopfield and D Tank.  Neural computation of decisions in optimization problems.  Biological cybernetics, 52:141–52, 02 1985.
  * Hopfield [2018] J.J. Hopfield.  Neural networks and physical systems with emergent collective computational abilities, pages 7–19.  01 2018.
  * Joshi et al. [2019] Chaitanya K Joshi, Thomas Laurent, and Xavier Bresson.  An efficient graph convolutional network technique for the travelling salesman problem.  arXiv preprint arXiv:1906.01227, 2019.
  * Khalil et al. [2016] Elias Khalil, Pierre Le Bodic, Le Song, George Nemhauser, and Bistra Dilkina.  Learning to branch in mixed integer programming.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 30, 2016.
  * Khalil et al. [2017] Elias B. Khalil, Bistra Dilkina, George L. Nemhauser, Shabbir Ahmed, and Yufen Shao.  Learning to run heuristics in tree search.  In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17, pages 659–666, 2017.
  * Kool et al. [2019] Wouter Kool, Herke van Hoof, and Max Welling.  Attention, learn to solve routing problems!  arXiv:1803.08475, 2019.
  * LeCun et al. [2015] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521:436–444, 2015.
  * Lodi and Zarpellon [2017] Andrea Lodi and Giulia Zarpellon.  On learning and branching: a survey.  Top, 25(2):207–236, Jul 2017.
  * Ma et al. [2019] Qiang Ma, Suwen Ge, Danyang He, Darshan Thaker, and Iddo Drori.  Combinatorial optimization by graph pointer networks and hierarchical reinforcement learning.  arXiv preprint arXiv:1911.04936, 2019.
  * Mandi and Guns [2020] Jayanta Mandi and Tias Guns.  Interior point solving for lp-based prediction+optimisation.  arXiv:2010.13943, 2020.
  * Mandi et al. [2020] Jayanta Mandi, Peter J Stuckey, Tias Guns, et al.  Smart predict-and-optimize for hard combinatorial optimization problems.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), volume 34, pages 1603–1610, 2020.
  * Nazari et al. [2018] Mohammadreza Nazari, Afshin Oroojlooy, Lawrence V Snyder, and Martin Takáč.  Reinforcement learning for solving the vehicle routing problem.  arXiv preprint arXiv:1802.04240, 2018.
  * Ndiaye et al. [2017] Eugene Ndiaye, Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon.  Gap safe screening rules for sparsity enforcing penalties.  The Journal of Machine Learning Research, 18(1):4671–4703, 2017\.
  * Ng et al. [2018] Yeesian Ng, Sidhant Misra, Line A Roald, and Scott Backhaus.  Statistical learning for dc optimal power flow.  In 2018 Power Systems Computation Conference (PSCC), pages 1–7. IEEE, 2018.
  * Nocedal and Wright [2006] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * Nowak et al. [2018] Alex Nowak, Soledad Villar, Afonso S. Bandeira, and Joan Bruna.  Revised note on learning algorithms for quadratic assignment with graph neural networks, 2018.
  * Pogančić et al. [2019] Marin Vlastelica Pogančić, Anselm Paulus, Vit Musil, Georg Martius, and Michal Rolinek.  Differentiation of blackbox combinatorial solvers.  In International Conference on Learning Representations, 2019.
  * Powell [1969] Michael JD Powell.  A method for nonlinear constraints in minimization problems.  Optimization, pages 283–298, 1969.
  * Prat and Chatzivasileiadis [2020] Eléa Prat and Spyros Chatzivasileiadis.  Learning active constraints to efficiently solve bilevel problems.  arXiv preprint arXiv:2010.06344, 2020.
  * Rossi et al. [2006] Francesca Rossi, Peter Van Beek, and Toby Walsh.  Handbook of constraint programming.  Elsevier, 2006.
  * Scarselli et al. [2008] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  IEEE transactions on neural networks, 20(1):61–80, 2008.
  * Song et al. [2020] Jialin Song, ravi lanka, Yisong Yue, and Bistra Dilkina.  A general large neighborhood search framework for solving integer linear programs.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 20012–20023. Curran Associates, Inc., 2020.
  * Sutton and Barto [2018] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * Tang et al. [2020] Yunhao Tang, Shipra Agrawal, and Yuri Faenza.  Reinforcement learning for integer programming: Learning to cut.  In International Conference on Machine Learning, pages 9367–9376. PMLR, 2020.
  * Tran et al. [2020] Cuong Tran, Ferdinando Fioretto, and Pascal Van Hentenryck.  Differentially private and fair deep learning: A lagrangian dual approach.  arXiv:2009.12562, 2020.
  * Veličković et al. [2017] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * Velloso and Van Hentenryck [2020] Alexandre Velloso and Pascal Van Hentenryck.  Combining deep learning and optimization for security-constrained optimal power flow.  arXiv:2007:2007.07002, 2020.
  * Vinyals et al. [2017] Oriol Vinyals, Meire Fortunato, and Navdeep Jaitly.  Pointer networks.  arXiv:1506.03134, 2017.
  * Wang et al. [2018] Po-Wei Wang, Wei-Cheng Chang, and J. Zico Kolter.  The mixing method: low-rank coordinate descent for semidefinite programming with diagonal constraints.  arXiv:1706.00476, 2018.
  * Wang et al. [2019] Po-Wei Wang, Priya Donti, Bryan Wilder, and Zico Kolter.  Satnet: Bridging deep learning and logical reasoning using a differentiable satisfiability solver.  In International Conference on Machine Learning, pages 6545–6554. PMLR, 2019.
  * Wilder et al. [2019] Bryan Wilder, Bistra Dilkina, and Milind Tambe.  Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), volume 33, pages 1658–1665, 2019.
  * Williams [1992] Ronald J Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  Machine learning, 8(3-4):229–256, 1992.
  * Wilson and Pawley [1988] GV Wilson and GS Pawley.  On the stability of the travelling salesman problem algorithm of hopfield and tank.  Biological Cybernetics, 58(1):63–70, 1988.
  * Zhou et al. [2018] Jie Zhou, Ganqu Cui, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph neural networks: A review of methods and applications.  arXiv preprint arXiv:1812.08434, 2018.


* Gibbons et al. (2014) Gibbons S. L. J., Belokurov V., Evans N. W., 2014, MNRAS, 445, 3788
  * Gómez et al. (2015) Gómez F. A., Besla G., Carpintero D. D., Villalobos Á., O’Shea B. W., Bell E. F., 2015, ApJ, 802, 128
  * Gravity Collaboration et al. (2018) Gravity Collaboration et al., 2018, A&A, 615, L15
  * Grillmair & Dionatos (2006) Grillmair C. J., Dionatos O., 2006, ApJ, 643, L17
  * Grillmair et al. (2015) Grillmair C. J., Hetherington L., Carlberg R. G., Willman B., 2015, ApJ, 812, L26
  * Harbeck et al. (2001) Harbeck D., et al., 2001, AJ, 122, 3092
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hattori et al. (2016) Hattori K., Erkal D., Sanders J. L., 2016, MNRAS, 460, 497
  * Hayes et al. (2020) Hayes C. R., et al., 2020, ApJ, 889, 63
  * Helmi (2020) Helmi A., 2020, ARA&A, 58, 205
  * Helmi et al. (1999) Helmi A., White S. D. M., de Zeeuw P. T., Zhao H., 1999, Nature, 402, 53
  * Hendel et al. (2018) Hendel D., et al., 2018, MNRAS, 479, 570
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Holl et al. (2018) Holl B., et al., 2018, A&A, 618, A30
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Husser et al. (2013) Husser T. O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6
  * Hyde et al. (2015) Hyde E. A., et al., 2015, ApJ, 805, 189
  * Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
  * Ibata et al. (2002) Ibata R. A., Lewis G. F., Irwin M. J., Quinn T., 2002, MNRAS, 332, 915
  * Iorio & Belokurov (2021) Iorio G., Belokurov V., 2021, MNRAS, 502, 5686
  * Jethwa et al. (2016) Jethwa P., Erkal D., Belokurov V., 2016, MNRAS, 461, 2212
  * Ji et al. (2021) Ji A. P., et al., 2021, ApJ, 921, 32
  * Jin et al. (2022) Jin S., et al., 2022, arXiv e-prints, p. arXiv:2212.03981
  * Johnston et al. (1999) Johnston K. V., Zhao H., Spergel D. N., Hernquist L., 1999, ApJ, 512, L109
  * Johnston et al. (2002) Johnston K. V., Spergel D. N., Haydn C., 2002, ApJ, 570, 656
  * Jönsson et al. (2020) Jönsson H., et al., 2020, AJ, 160, 120
  * Kallivayalil et al. (2013) Kallivayalil N., van der Marel R. P., Besla G., Anderson J., Alcock C., 2013, ApJ, 764, 161
  * Kirby et al. (2013) Kirby E. N., Cohen J. G., Guhathakurta P., Cheng L., Bullock J. S., Gallazzi A., 2013, ApJ, 779, 102
  * Koleva et al. (2011) Koleva M., Prugniel P., De Rijcke S., Zeilinger W. W., 2011, MNRAS, 417, 1643
  * Koposov (2021) Koposov S., 2021, segasai/minimint: Minimint 0.3.0, doi:10.5281/zenodo.5610692, https://doi.org/10.5281/zenodo.5610692
  * Koposov (2022b) Koposov S., 2022b, segasai/sqlutilpy: sqlutilpy v0.19.0, doi:10.5281/zenodo.6867957, https://doi.org/10.5281/zenodo.6867957
  * Koposov (2022a) Koposov S., 2022a, segasai/stan-splines: v1.0.0, doi:10.5281/zenodo.7193910, https://doi.org/10.5281/zenodo.7193910
  * Koposov & Bartunov (2006) Koposov S., Bartunov O., 2006, in Gabriel C., Arviset C., Ponz D., Enrique S., eds, Astronomical Society of the Pacific Conference Series Vol. 351, Astronomical Data Analysis Software and Systems XV. p. 735
  * Koposov et al. (2010) Koposov S. E., Rix H.-W., Hogg D. W., 2010, ApJ, 712, 260
  * Koposov et al. (2014) Koposov S. E., Irwin M., Belokurov V., Gonzalez-Solares E., Yoldas A. K., Lewis J., Metcalfe N., Shanks T., 2014, MNRAS, 442, L85
  * Koposov et al. (2019) Koposov S. E., et al., 2019, MNRAS, 485, 4726
  * Koppelman et al. (2018) Koppelman H., Helmi A., Veljanoski J., 2018, ApJ, 860, L11
  * Krishnarao et al. (2022) Krishnarao D., et al., 2022, Nature, 609, 915
  * Küpper et al. (2015) Küpper A. H. W., Balbinot E., Bonaca A., Johnston K. V., Hogg D. W., Kroupa P., Santiago B. X., 2015, ApJ, 803, 80
  * Lee et al. (2008a) Lee Y. S., et al., 2008a, AJ, 136, 2022
  * Lee et al. (2008b) Lee Y. S., et al., 2008b, AJ, 136, 2050
  * Lewis et al. (2002) Lewis I. J., et al., 2002, MNRAS, 333, 279
  * Li & S5 Collaboration (2021) Li T., S5 Collaboration 2021, Southern Stellar Stream Spectroscopic Survey: The First Public Data Release, doi:10.5281/zenodo.4695135, https://doi.org/10.5281/zenodo.4695135
  * Li et al. (2019) Li T. S., et al., 2019, MNRAS, 490, 3508
  * Li et al. (2021) Li T. S., et al., 2021, ApJ, 911, 149
  * Li et al. (2022) Li T. S., et al., 2022, ApJ, 928, 30
  * Lilleengen et al. (2023) Lilleengen S., et al., 2023, MNRAS, 518, 774
  * Malhan et al. (2018) Malhan K., Ibata R. A., Martin N. F., 2018, MNRAS, 481, 3442
  * Mateu (2022) Mateu C., 2022, arXiv e-prints, p. arXiv:2204.10326
  * McConnachie (2012) McConnachie A. W., 2012, AJ, 144, 4
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Mendelsohn et al. (2022) Mendelsohn E. J., Newberg H. J., Shelton S., Widrow L. M., Thompson J. M., Grillmair C. J., 2022, ApJ, 926, 106
  * Mercado et al. (2021) Mercado F. J., et al., 2021, MNRAS, 501, 5121
  * Miyamoto & Nagai (1975) Miyamoto M., Nagai R., 1975, PASJ, 27, 533
  * Muraveva et al. (2018) Muraveva T., Delgado H. E., Clementini G., Sarro L. M., Garofalo A., 2018, MNRAS, 481, 1195
  * Myeong et al. (2018) Myeong G. C., Evans N. W., Belokurov V., Sanders J. L., Koposov S. E., 2018, ApJ, 856, L26
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Nidever et al. (2015) Nidever D. L., et al., 2015, AJ, 150, 173
  * Odenkirchen et al. (2003) Odenkirchen M., et al., 2003, AJ, 126, 2385
  * Okamoto et al. (2017) Okamoto S., Arimoto N., Tolstoy E., Jablonka P., Irwin M. J., Komiyama Y., Yamada Y., Onodera M., 2017, MNRAS, 467, 208
  * Pace et al. (2022) Pace A. B., Erkal D., Li T. S., 2022, ApJ, 940, 136
  * Patrick et al. (2022) Patrick J. M., Koposov S. E., Walker M. G., 2022, MNRAS, 514, 1757
  * Peñarrubia et al. (2012) Peñarrubia J., Koposov S. E., Walker M. G., 2012, ApJ, 760, 2
  * Pérez & Granger (2007) Pérez F., Granger B. E., 2007, Computing in Science and Engineering, 9, 21
  * Petersen & Peñarrubia (2021) Petersen M. S., Peñarrubia J., 2021, Nature Astronomy, 5, 251
  * Piatek et al. (2016) Piatek S., Pryor C., Olszewski E. W., 2016, AJ, 152, 166
  * Pietrzyński et al. (2019) Pietrzyński G., et al., 2019, Nature, 567, 200
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Price-Whelan et al. (2016) Price-Whelan A. M., Sesar B., Johnston K. V., Rix H.-W., 2016, ApJ, 824, 104
  * Reid & Brunthaler (2004) Reid M. J., Brunthaler A., 2004, ApJ, 616, 872
  * Reino et al. (2021) Reino S., Rossi E. M., Sanderson R. E., Sellentin E., Helmi A., Koppelman H. H., Sharma S., 2021, MNRAS, 502, 4170
  * Sanders & Binney (2013) Sanders J. L., Binney J., 2013, MNRAS, 433, 1813
  * Sanderson et al. (2015) Sanderson R. E., Helmi A., Hogg D. W., 2015, ApJ, 801, 98
  * Schee et al. (2013) Schee J., Stuchlík Z., Petrásek M., 2013, J. Cosmology Astropart. Phys., 2013, 026
  * Sesar et al. (2013) Sesar B., et al., 2013, ApJ, 776, 26
  * Sesar et al. (2017) Sesar B., et al., 2017, AJ, 153, 204
  * Sharp et al. (2006) Sharp R., et al., 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 62690G (arXiv:astro-ph/0606137), doi:10.1117/12.671022
  * Shipp et al. (2018) Shipp N., et al., 2018, ApJ, 862, 114
  * Shipp et al. (2019) Shipp N., et al., 2019, ApJ, 885, 3
  * Shipp et al. (2021) Shipp N., et al., 2021, ApJ, 923, 149
  * Smolinski et al. (2011) Smolinski J. P., et al., 2011, AJ, 141, 89
  * Sohn et al. (2013) Sohn S. T., Besla G., van der Marel R. P., Boylan-Kolchin M., Majewski S. R., Bullock J. S., 2013, ApJ, 768, 139
  * Stringer et al. (2021) Stringer K. M., et al., 2021, ApJ, 911, 109
  * van der Marel & Kallivayalil (2014) van der Marel R. P., Kallivayalil N., 2014, ApJ, 781, 121
  * van der Marel et al. (2002) van der Marel R. P., Alves D. R., Hardy E., Suntzeff N. B., 2002, AJ, 124, 2639
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 482, 1525
  * Vasiliev & Baumgardt (2021) Vasiliev E., Baumgardt H., 2021, MNRAS, 505, 5978
  * Vasiliev et al. (2021) Vasiliev E., Belokurov V., Erkal D., 2021, MNRAS, 501, 2279
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wan et al. (2020) Wan Z., Guglielmo M., Lewis G. F., Mackey D., Ibata R. A., 2020, MNRAS, 492, 782
  * White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
  * Wu et al. (2008) Wu X., Famaey B., Gentile G., Perets H., Zhao H., 2008, MNRAS, 386, 2199
  * Yoon et al. (2011) Yoon J. H., Johnston K. V., Hogg D. W., 2011, ApJ, 731, 58
  * Zhao et al. (2012) Zhao G., Zhao Y.-H., Chu Y.-Q., Jing Y.-P., Deng L.-C., 2012, Research in Astronomy and Astrophysics, 12, 723

## APPENDIX A

Table 3 contains the list of possible spectroscopic members from ${S}^{5}$,
APOGEE, SDSS and LAMOST selected based on position, proper motion, metallicity
and radial velocity as shown on Figure 1. The stars in that table were not
filtered by colour-magnitude, but we provide a column indicating whether the
star is in the CMD mask shown on Figure 3.

GDR3 source_id | $\alpha$ | $\delta$ | $\phi_{1}$ | $\phi_{2}$ | RV | $V_{GSR}$ | $\sigma_{V}$ | $\mu_{\phi,1}$ | $\mu_{\phi,2}$ | CMD
---|---|---|---|---|---|---|---|---|---|---
| deg | deg | deg | deg | origin | km s-1 | km s-1 | mas yr-1 | mas yr-1 | selected
3888500490279502464 | 155.87168 | 15.86127 | 83.3819 | 1.9210 | S5 | 121.2 | 1.3 | 0.36 | 1.23 | 1
6451721201411909760 | 318.52474 | -62.09448 | -48.9244 | 1.8850 | S5 | -203.3 | 3.8 | 2.65 | 1.23 | 1
6558441247408890240 | 327.06713 | -52.19549 | -59.8363 | 2.1710 | S5 | -191.8 | 0.2 | 2.18 | 1.18 | 1
Table 3: Possible spectroscopic OC member stars from S5, APOGEE, LAMOST and
SDSS datasets used in the analysis. The full table is provided in
supplementary materials. The source_id is Gaia DR3 identifier. The star
coordinates are are given in equatorial and $\phi_{1},\phi_{2}$ coordinate
systems. We also provide the radial velocity in the GSR frame with the
uncertainty as well as the survey that is the source of the radial velocity.
The proper motions are provided in the $\phi_{1},\phi_{2}$ frame. The last
column indicates if the star was selected through with the CMD mask shown on
Figure 3.

## APPENDIX B:

Here, we provide the table of measurements from the spline models used to
extract the stream 6-D track. The full versions of these tables are also
provided as FITS files in supplementary materials, and can be used as inputs
to natural cubic splines (i.e. CubicSpline with bc_type=’natural’ option from
scipy.interpolate in Python programming language). The stream distance
measurements are in Table 4, while distance spread measurements are in Table
5. The radial velocity measurements are provided in Table 6. The stream
velocity dispersion is provided in Table 7. The proper motion stream track
measurements are in Table 8. The intrinsic dispersions in proper motions are
provided in Table 9. The stream track on the sky is provided in Table 10. The
stream width measurements are in Table 11. The stream surface brightness
measurements are in Table 12. The values of nuisance noise parameters used for
the stream modelling are given in Table 13.

$\phi_{1}$ | $m-M$ | $\sigma_{m-M}$
---|---|---
deg | mag | mag
-85.0 | 19.42 | 0.09
-70.4 | 18.33 | 0.02
-55.7 | 17.65 | 0.03
Table 4: Measurements of the stream distance modulus through spline models. We provide the distance modulus and it uncertainty at the set of $\phi_{1}$ points. The full table is included in the supplementary materials. $\phi_{1}$ | $\ln{\mathcal{S}}_{m-M}$ | $\sigma_{\ln{\mathcal{S}},{m-M}}$
---|---|---
deg | mag | mag
-85.0 | -3.10 | 0.60
-44.0 | -2.59 | 0.25
-3.0 | -2.57 | 0.48
Table 5: Measurements of the intrinsic spread of stream distance moduli through spline models. We provide the logarithm of $m-M$ spread and its uncertainty. The full table is included in the supplementary materials. $\phi_{1}$ | $V_{GSR}$ | $\sigma_{V}$
---|---|---
deg | km s-1 | km s-1
-80.0 | -132.4 | 5.7
-70.0 | -171.8 | 2.1
-60.0 | -187.8 | 1.1
Table 6: Measurements of stream radial velocity track through spline models. We provide the radial velocity in the GSR frame and its uncertainty at the set of $\phi_{1}$ points. The full table is available in the supplementary materials. $\phi_{1}$ | $\ln{\mathcal{S}}_{V}$ | $\sigma_{\ln{\mathcal{S}},V}$
---|---|---
deg |  |
-80.0 | 1.4 | 0.5
-40.0 | 1.6 | 0.3
40.0 | 1.2 | 0.2
80.0 | 1.7 | 0.2
120.0 | 1.9 | 0.4
Table 7: Measurements of the radial velocity dispersion through spline models. We provide measurements of the natural log of the velocity dispersion in $\mathrm{\,km\,s^{-1}}$ and its uncertainty at the set of $\phi_{1}$ points. The full table is included in the supplementary materials. $\phi_{1}$ | $\mu_{1}$ | $\sigma_{\mu,1}$ | $\mu_{2}$ | $\sigma_{\mu,2}$
---|---|---|---|---
deg | mas yr-1 | mas yr-1 | mas yr-1 | mas yr-1
-80.0 | 1.08 | 0.07 | 0.54 | 0.06
-70.0 | 1.65 | 0.03 | 0.92 | 0.03
-60.0 | 2.15 | 0.02 | 1.15 | 0.01
Table 8: Measurement of stream proper motion through spline models. We provide proper motions in $\phi_{1}$ and $\phi_{2}$ (not corrected for Solar reflex motion) and their uncertainties at the set of $\phi_{1}$ points. The full table is included in the supplementary materials. $\phi_{1}$ | $\ln{\mathcal{S}}_{\mu,1}$ | $\sigma_{\ln{\mathcal{S}},\mu,1}$ | $\ln{\mathcal{S}}_{\mu,2}$ | $\sigma_{\ln{\mathcal{S}},\mu,2}$
---|---|---|---|---
deg |  |  |  |
-80.0 | -3.4 | 1.0 | -4.1 | 1.1
-50.0 | -3.2 | 0.3 | -3.4 | 0.5
-25.0 | -2.0 | 0.3 | -3.2 | 0.5
Table 9: Measurement of stream proper motion dispersion. We provide measurements of the natural logarithm of the proper motion dispersion in $\phi_{1}$ and $\phi_{2}$ in $\mathrm{\,mas\,yr^{-1}}$ together with the uncertainties. The full table is included in the supplementary materials. $\phi_{1}$ | $\phi_{2}$ | $\sigma_{\phi,2}$
---|---|---
deg | deg | deg
-90.0 | -1.79 | 2.80
-80.0 | 0.70 | 0.38
-70.0 | 1.69 | 0.22
Table 10: Measurement of stream track on the sky through spline models. We provide the $\phi_{2}$ coordinates with the corresponding uncertainty for a set of $\phi_{1}$ locations. The full table is included in the supplementary materials. $\phi_{1}$ | $\ln\Sigma_{\phi,2}$ | $\sigma_{\ln\Sigma,\phi,2}$
---|---|---
deg |  |
-90.0 | -0.13 | 0.28
-60.0 | 0.14 | 0.13
-30.0 | -0.29 | 0.19
Table 11: Measurements of stream width on the sky through spline models. We provide the value of the natural logarithm of stream width in degrees and its uncertainty at at the set of $\phi_{1}$ positions. The full table is included in supplementary materials. $\phi_{1}$ | $\ln I$ | $\sigma_{\ln I}$
---|---|---
deg |  |
-90.0 | -11.83 | 2.99
-80.0 | -4.14 | 0.36
-70.0 | -3.12 | 0.20
Table 12: Measurements of stream surface brightness from spline models (see Eq. 5). We provide the logarithm of surface brightness and its uncertainty at the set of $\phi_{1}$ positions. The units of surface brightness are stars per bin. The full table is included in the supplementary materials. Parameter | Spherical | Prolate | Oblate
---|---|---|---
$\sigma_{\phi_{2}}$ (deg) | $0.74^{+0.16}_{-0.17}$ | $0.41^{+0.19}_{-0.17}$ | $0.39^{+0.17}_{-0.15}$
$\sigma_{\rm DM}$ ($\mathrm{\,mas\,yr^{-1}}$) | $0.09^{+0.04}_{-0.03}$ | $0.10^{+0.03}_{-0.03}$ | $0.07^{+0.03}_{-0.02}$
$\sigma_{\mu_{1}^{*}}$ ($\mathrm{\,mas\,yr^{-1}}$) | $0.05^{+0.03}_{-0.02}$ | $0.04^{+0.02}_{-0.02}$ | $0.05^{+0.02}_{-0.02}$
$\sigma_{\mu_{2}}$ ($\mathrm{\,mas\,yr^{-1}}$) | $0.05^{+0.02}_{-0.01}$ | $0.04^{+0.02}_{-0.02}$ | $0.02^{+0.02}_{-0.01}$
$\sigma_{v_{r}}$ ($\mathrm{\,km\,s^{-1}}$) | $4.93^{+3.71}_{-2.70}$ | $2.51^{+2.18}_{-1.35}$ | $3.18^{+4.84}_{-1.76}$
Table 13: Values of nuisance parameters for the additional scatter in proper
motions, radial velocities and stream track determined as part of the model
presented in Section 5.

## APPENDIX C:

In supplementary materials we also provide fiducial splines used to select
likely members (shown on Figure 1). The files contain sets of points
$\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2},V_{gsr}$ and distance defining
fiducial stream tracks. The actual tracks shown on Figure 1 can be obtained by
linear interpolation of provided points.

## APPENDIX D:

In the supplementary materials, we include Stan models used in this paper.

  * •

rvfit.stan is the radial velocity model used in Section 3.5

  * •

pmfit.stan is the proper motion model used in Section 3.6

  * •

distfit.stan is the distance modulus model for RR Lyrae used in Section 3.3.

  * •

fehmodel.stan is the model for the [Fe/H] gradients used in Section 4.2.

  * •

density_pm_fast_reduce.stan is the proper motion/spatial distribution model
described in Section 3.7. The model is parallisable across multiple cores.

## APPENDIX E:

In Figure 27 we show the evolution of the OC stream over the past 1 Gyr in
energy versus angular momentum.

Figure 27: Evolution of the OC stream in energy and angular momentum. From
left to right, the panels show the energy versus angular momentum for stream
particles at look-back times of $1000$ Myr (i.e. well before the LMC has an
effect), $420$ Myr ($\sim 50$ Myr before the LMC’s closest approach), $320$
Myr ($\sim 50$ Myr after the LMC’s closest approach), and the present day. A
movie showing the evolution of energy and angular momentum is available here.

## APPENDIX F:

The residuals for our best-fit model presented in Figure 14 are shown in
Figure 28. The residuals are mostly small with the exception of the rightmost
RV measurement and the points at the very left edge of the stream that are the
most violently affected by the LMC interaction. Grey points show the points
that were not part of the fit.

Figure 28: Residuals of the best-fit model shown in Figure 14. From top to
bottom, the panels show the residuals of the track on the sky, the distance
modulus, the two proper motions, and the radial velocity. This fit has a
$\chi^{2}$ per degree of freedom of 2.5, suggesting that our best fit model is
not fully describing the data.