Exploring the Power of Graph Neural Networks in Solving Linear
  Optimization Problems
Recently, machine learning, particularly message-passing graph neural
networks (MPNNs), has gained traction in enhancing exact optimization
algorithms. For example, MPNNs speed up solving mixed-integer optimization
problems by imitating computational intensive heuristics like strong branching,
which entails solving multiple linear optimization problems (LPs). Despite the
empirical success, the reasons behind MPNNs' effectiveness in emulating linear
optimization remain largely unclear. Here, we show that MPNNs can simulate
standard interior-point methods for LPs, explaining their practical success.
Furthermore, we highlight how MPNNs can serve as a lightweight proxy for
solving LPs, adapting to a given problem instance distribution. Empirically, we
show that MPNNs solve LP relaxations of standard combinatorial optimization
problems close to optimality, often surpassing conventional solvers and
competing approaches in solving time.
Chendi Qian, Didier Chételat, Christopher Morris
# Exploring the Power of Graph Neural Networks in Solving Linear Optimization Problems

###### Abstract

Recently, machine learning, particularly message-passing graph neural networks (MPNNs), has gained traction in enhancing exact optimization algorithms. For example, MPNNs speed up solving mixed-integer optimization problems by imitating computational intensive heuristics like strong branching, which entails solving multiple linear optimization problems (LPs). Despite the empirical success, the reasons behind MPNNs' effectiveness in emulating linear optimization remain largely unclear. Here, we show that MPNNs can simulate standard interior-point methods for LPs, explaining their practical success. Furthermore, we highlight how MPNNs can serve as a lightweight proxy for solving LPs, adapting to a given problem instance distribution. Empirically, we show that MPNNs solve LP relaxations of standard combinatorial optimization problems close to optimality, often surpassing conventional solvers and competing approaches in solving time.

## 1 Introduction

Recently, there has been a surge of interest in training _message-passing graph neural networks_ (MPNNs) to imitate steps of classical algorithms, such as for shortest-path problems Cappart et al. (2021); Velickovic et al. (2020). As many of those problems can be formulated as linear optimization problems (LPs), it is natural to ask whether MPNNs could be trained to solve general LPs, at least approximately.

Another recent line of research makes this question particularly intriguing. In integer linear optimization, state-of-the-art solvers all rely on the branch-and-bound algorithm, in which one must repeatedly select variables, subdividing the search space. The best-known heuristic for variable selection is known as "strong branching," which entails solving LPs to score the variables. This heuristic is, unfortunately, too computationally expensive to use in practice. However, in recent years, there has been a collection of works (Gasse et al., 2019; Gupta et al., 2020, 2022; Nair et al., 2020; Seyfi et al., 2023) that have proposed to use MPNNs to imitate strong branching, with impressive success. No theoretical explanation has ever been put forward to explain this success. However, perhaps the most straightforward explanation would be that MPNNs implicitly learn to imitate the LP solving underlying strong branching.

Chen et al. (2023) provide a first step towards an explanation. In this work, the authors propose to encode LPs as bipartite graphs in the spirit of (Gasse et al., 2019) and show that MPNNs, in principle, can learn to predict the optimal solution of an LP instance up to arbitrary small \(\varepsilon\) concerning the supremum norm. They also provide some small-scale experiments suggesting that MPNNs can learn to approximate LP solutions surprisingly well. While a step in the right direction, their theoretical result heavily relies on invoking the universal approximation theorem for multi-layer perceptrons Cybenko (1992); Leshno et al. (1993), and therefore does not explain why modestly-sized MPNNs could be particularly effective LP solvers in practice.

In this paper, we present instead a more specific explanation. We show that several variants of _interior-point methods_ (IPMs) Gondzio (2012); Nocedal and Wright (2006), an established class of polynomial-time algorithms for solving LPs, can be interpreted as an MPNN with a specific architecture and choice of parameters that takes as input a graph representation of the LP. Specifically, for two common IPM variants, we show that a sequence of standard MPNN steps can emulate a single iteration of the algorithm on a tripartite (rather than bipartite) graph representation of the LP. This novel theoretical result suggests several conclusions. First, it indicates that MPNNs appear successful at imitating LP solving because they might often be imitating IPMs and that there is a close connection between MPNNs and this specific LP algorithm. Secondly, although the MPNN's architecture in our theoretical result involves many MPNN layers, MPNNs are a natural choice of machine learning model for LP solving. It is likely possible to approximately solve LPs with fewer layers.

To prove the last hypothesis, we train MPNNs with sev
eral layers less than predicted by our theoretical results to imitate the output of a practical IPM algorithm for LP solving, resulting in the _IPM-MPNN_ architecture. Our empirical results show that IPM-MPNNs can lead to reduced solving times compared to a state-of-the-art LP solver with time constraints and competing neural network-based approaches; see Figure 0(b) for an overview of our approach.

_In summary, our findings significantly contribute to the theoretical framework of data-driven exact optimization using MPNNs, and we also showcase the potential of MPNNs in serving as a light-weight proxy for solving LPs._

### Additional related works

In the following, we discuss relevant related work.

**MPNNs** MPNNs (Gilmer et al., 2017; Scarselli et al., 2009) have emerged as a flexible framework for machine learning on graphs and relational data. Notable instances of this architecture include, e.g., Duvenaud et al. (2015); Hamilton et al. (2017); Velickovic et al. (2018), and the spectral approaches proposed in, e.g., Bruna et al. (2014); Defferrard et al. (2016); Kipf and Welling (2017)--all of which descend from early work in Baskin et al. (1997); Kireev (1995); Merkwirth and Lengauer (2005); Micheli (2009); Micheli and Sestito (2005); Scarselli et al. (2009); Sperduti and Starita (1997).

**Machine learning for combinatorial optimization**Bengio et al. (2021) discuss and review machine learning approaches to enhance combinatorial optimization (CO). Concrete examples include the imitation of computationally intensive variable selection rules within the branch-and-cut framework (Khalil et al., 2016; Zarpellon et al., 2020), learning to run (primal) heuristics (Khalil et al., 2017; Chmiela et al., 2021), learning decompositions of large MILPs for scalable heuristic solving (Song et al., 2020), learning to generate cutting planes (Deza and Khalil, 2023) or leveraging machine learning to find (primal) solutions to stochastic integer problems quickly (Bengio et al., 2020)--see Kotary et al. (2021) for a high-level overview of recent advances.

**MPNNs for CO** Many prominent CO problems involve graph or relational structures, either directly given as input or induced by the variable-constraint interactions. Recent progress in using MPNNs to bridge the gap between machine learning and combinatorial optimization is surveyed in Capppart et al. (2021).

Most relevant to the present work, Gasse et al. (2019) proposed to encode the variable-constraint interaction of a mixed-integer linear optimization problem as a bipartite graph and trained MPNNs in a supervised fashion to imitate the costly strong branching heuristic, which entails solving multiple linear optimization problems, within the branch-and-cut framework (Achterberg et al., 2005). Building on that, Gupta et al. (2020) proposed a hybrid branching model using an MPNN at the initial decision point and a light multi-layer perceptron for subsequent steps, showing improvements on pure CPU machines. Subsequently, Nair et al. (2020) expanded the MPNN approach to branching by implementing a GPU-friendly parallel linear programming solver using the alternating direction method of multipliers that allows scaling the strong branching expert to substantially larger instances, also combining this innovation with a novel MPNN approach to diving. For all the above works, it remains largely unclear why MPNNs are good at (approximately) predicting strong branching scores. Moreover, Khalil et al. (2022) used MPNNs to predict the probability of variables being assigned to 0 or 1 in near-optimal solutions of binary-integer linear optimization problems.

Ding et al. (2020) used MPNNs on a tripartite graph consisting of variables, constraints, and a single objective node enriched with hand-crafted node features. The target is to predict the 0-1 values of the so-called _stable variables_, i.e., variables whose assignment does not change over a set of pre-computed feasible solutions. Li et al. (2022) used MPNNs on a bipartite graph together with a pointer network (Bello et al., 2016) to reorder the variables of a given LP instance, resulting in reduced solving time. Fan et al. (2023) leveraged MPNNs to find good initial basis solutions for the Simplex algorithm.

Finally Wu and Lisser (2023) expressed an LP as an ordinary differential equations system whose state solution converges to the LP's optimal solution and trained a feed-forward neural network to approximate this state solution. However, their approach hinges on the need to compute the Jacobian matrix, rendering its training phase computationally costly.

## 2 Background

In the following, we describe the necessary background.

**Notation** Let \(\mathbb{N}\coloneqq\{1,2,3,\dots\}\). For \(n\geq 1\), let \([n]\coloneqq\{1,\dots,n\}\subset\mathbb{N}\). We use \(\{\!\!\{\dots\}\!\}\) to denote multisets, i.e., the generalization of sets allowing for multiple instances for each of its elements. A _graph_\(G\) is a pair \((V(G),E(G))\) with _finite_ sets of _vertices_ or _nodes_\(V(G)\) and _edges_\(E(G)\subseteq\{\{u,v\}\subseteq V(G)\mid u\neq v\}\). For ease of notation, we denote the edge \(\{u,v\}\) in \(E(G)\) by \((u,v)\) or \((v,u)\). Throughout the paper, we use standard notation, e.g., we denote the _neighborhood_ of a node \(v\) by \(N(v)\), and so on; see Appendix A for details. Moreover,
let \(\mathbf{x}\in\mathbb{R}^{1\times d}\), then \(\mathbf{D}(\mathbf{x})\) denotes the diagonal matrix with diagonal \(\mathbf{x}\), \(\mathbf{0}\) and \(\mathbf{1}\) denote the vector of zero and ones, respectively, with an appropriate number of entries. By default, a vector \(\mathbf{x}\in\mathbb{R}^{d}\) is a column vector.

**Linear optimization problems** A _linear optimization problem_ (LP) aims at optimizing a linear function over a feasible set described as the intersection of finitely many half-spaces, i.e., a polyhedron. We restrict our attention to feasible and bounded LPs. Formally, an instance \(I\) of an LP is a tuple \((\mathbf{A},\mathbf{b},\mathbf{c})\), where \(\mathbf{A}\) is a matrix in \(\mathbb{Q}^{m\times n}\), and \(\mathbf{b}\) and \(\mathbf{c}\) are vectors in \(\mathbb{Q}^{m}\) and \(\mathbb{Q}^{n}\), respectively. We aim at finding a vector \(\mathbf{x}^{*}\) in \(\mathbb{Q}^{n}\) that minimizes \(\mathbf{c}^{\mathsf{T}}\mathbf{x}^{*}\) over the _feasible set_

\[F(I)=\{\mathbf{x}\in\mathbb{Q}^{n}\mid \mathbf{A}_{j}\mathbf{x}\leq b_{j}\text{ for }j\in[m]\text{ and } \tag{1}\] \[x_{i}\geq 0\text{ for }i\in[n]\}.\]

In practice, LPs are solved using the Simplex method or polynomial-time IPMs (Nocedal and Wright, 2006).

We now detail the theoretical result that motivates our approach. We first summarize both interior point methods for linear optimization and MPNNs and then prove a theorem that relates the two.

**Interior-point methods for linear optimization** IPMs are algorithms for solving constrained optimization problems. They are particularly efficient for linear optimization, where they were first developed as a (polynomial-time) alternative to the Simplex methods. Variants of the algorithm differ in theoretical guarantees and empirical performance but revolve around the same core approach (Shanno, 2012). In short, the LP to solve is replaced by a perturbed family of problems where a barrier penalty has replaced hard constraints with a parameter \(\mu>0\). IPMs alternate between taking a Newton step to solve this perturbed problem and decreasing this parameter \(\mu>0\), eventually converging to the optimal solution of the original problem.

For concreteness, we present two variants of the approach, an algorithm with theoretical guarantees and a practical algorithm that could be used in practice, following Nocedal and Wright (2006, Chapter 14) and Gondzio (2012), respectively. In the next section, we will show that both algorithms can be connected to MPNNs.

The core idea of IPMs is as follows. First, we consider a perturbed version of the LP (1),

\[\min_{\mathbf{x}\in\mathbb{Q}^{n}}\mathbf{c}^{\mathsf{T}}\mathbf{x}-\mu[\mathbf{1}^{\mathsf{T }}\log(\mathbf{A}\mathbf{x}-\mathbf{b})+\mathbf{1}^{\mathsf{T}}\log(\mathbf{x})] \tag{2}\]

for some \(\mu>0\). By introducing the variables \(s_{i}=\mu/x_{i}\), \(r_{j}=\mathbf{A}_{j}x-b_{j}\) and \(w_{j}=\mu/r_{j}\), the first-order optimality conditions for (2) can be written as a system

\[\mathbf{A}\mathbf{x}^{*}-\mathbf{r}^{*} =\mathbf{b}\] \[\mathbf{A}^{\mathsf{T}}\mathbf{w}^{*}+\mathbf{s}^{*} =\mathbf{c}\] \[x_{i}^{*}s_{i}^{*} =\mu i\in[n],\] \[w_{i}^{*}r_{i}^{*} =\mu j\in[m],\]

with \(\mathbf{x}^{*},\mathbf{w}^{*},\mathbf{s}^{*},\mathbf{r}^{*}\geq\mathbf{0}\). Let \(\sigma\in(0,1)\) be another hyperparameter. The two algorithms start from an initial positive point \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})>\mathbf{0}\), and alternate between computing the Newton step for the perturbed problem (2) at barrier parameter \(\sigma\mu\)

\[\begin{bmatrix}\mathbf{A}&\mathbf{0}&\mathbf{0}&-\mathbf{I}\\ \mathbf{0}&\mathbf{A}^{\mathsf{T}}&\mathbf{I}&\mathbf{0}\\ \mathbf{D}(\mathbf{s})&\mathbf{0}&\mathbf{D}(\mathbf{x})&\mathbf{0}\\ \mathbf{0}&\mathbf{D}(\mathbf{r})&\mathbf{0}&\mathbf{D}(\mathbf{w})\end{bmatrix}\begin{bmatrix}\Delta \mathbf{x}\\ \Delta\mathbf{w}\\ \Delta\mathbf{s}\\ \Delta\mathbf{x}\end{bmatrix}=\begin{bmatrix}\mathbf{b}-\mathbf{A}\mathbf{x}+\mathbf{r}\\ \mathbf{c}-\mathbf{A}^{\mathsf{T}}\mathbf{w}-\mathbf{s}\\ \sigma\mu\mathbf{1}-\mathbf{D}(\mathbf{x})\mathbf{D}(\mathbf{s})\mathbf{1}\\ \sigma\mu\mathbf{1}-\mathbf{D}(\mathbf{w})\mathbf{D}(\mathbf{r})\mathbf{1}\end{bmatrix},\]

and taking a step in that direction with length \(\alpha>0\), such that the resulting point satisfies \((\mathbf{x}^{\prime},\mathbf{w}^{\prime},\mathbf{s}^{\prime},\mathbf{r}^{\prime})>\mathbf{0}\).

The above system can be simplified as follows. First, we can infer that

\[\Delta\mathbf{s} =\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}-\mathbf{s}-\mathbf{D}(\mathbf{x})^{-1}\bm {D}(\mathbf{s})\Delta\mathbf{x}, \tag{3}\] \[\Delta\mathbf{r} =\sigma\mu\mathbf{D}(\mathbf{w})^{-1}\mathbf{1}-\mathbf{r}-\mathbf{D}(\mathbf{w})^{-1}\bm {D}(\mathbf{r})\Delta\mathbf{w}, \tag{4}\]

Figure 1: Overview of our IPM-MPNN framework.
which implies that

\[\mathbf{A}\Delta\mathbf{x}+\mathbf{D}(\mathbf{w})^{-1}\mathbf{D}(\mathbf{r})\Delta\mathbf{w} =\mathbf{b}-\mathbf{A}\mathbf{x}+\sigma\mu\mathbf{D}(\mathbf{w})^{-1}\mathbf{1},\] \[\mathbf{A}^{\mathsf{T}}\Delta\mathbf{w}-\mathbf{D}(\mathbf{x})^{-1}\mathbf{D}(\mathbf{s}) \Delta\mathbf{x} =\mathbf{c}-\mathbf{A}^{\mathsf{T}}\mathbf{w}-\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}.\]

Therefore, \(\Delta\mathbf{x}=\)

\[\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})[\mathbf{A}^{\mathsf{T}}\Delta\mathbf{w} -\mathbf{c}+\mathbf{A}^{\mathsf{T}}\mathbf{w}+\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}] \tag{5}\] \[\mathbf{Q}\Delta\mathbf{w} =\mathbf{b}-\mathbf{A}\mathbf{x}+\sigma\mu\mathbf{D}(\mathbf{w})^{-1}1\] (6) \[+\mathbf{A}\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})[\mathbf{c}-\mathbf{A}^{\mathsf{ T}}\mathbf{w}-\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}]\]

for \(\mathbf{Q}=\mathbf{A}\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})\mathbf{A}^{\mathsf{T}}+\mathbf{D}(\mathbf{ w})^{-1}\mathbf{D}(\mathbf{r})\). Thus, by solving the linear system in Equation (6), we can find \(\Delta\mathbf{w}\), then \(\Delta\mathbf{x}\), \(\Delta\mathbf{r}\), and \(\Delta\mathbf{s}\) through Equations (3)-(5).

The two algorithms we consider only differ in how they compute \(\mu\) and \(\alpha\). The theoretical algorithm recomputes, at every iteration, \(\mu=(\mathbf{x}^{\mathsf{T}}\mathbf{s}+\mathbf{w}^{\mathsf{T}}\mathbf{r})/(n+m)\), and chooses \(\alpha\) to be the largest \(\alpha<1\) such that \(x_{i}^{\prime}s_{i}^{\prime}\geq\gamma(\mathbf{x}^{\mathsf{T}}\mathbf{s}^{\prime}+\bm {w}^{\mathsf{T}}\mathbf{r}^{\prime})/(n+m)\) and \(w_{i}^{\prime}r_{j}^{\prime}\geq\gamma(\mathbf{x}^{\mathsf{T}}\mathbf{s}^{\prime}+\bm {w}^{\mathsf{T}}\mathbf{r}^{\prime})/(n+m)\) for \(i\in[n]\) and \(j\in[m]\) for some hyperparameter \(\gamma\in(0,1]\)(Nocedal and Wright, 2006, Algorithm 14.2). The practical algorithm instead picks \(\mu_{0}=(\mathbf{x}_{0}^{\mathsf{T}}\mathbf{s}_{0}+\mathbf{w}_{0}^{\mathsf{T}}\mathbf{r}_{0}) /(n+m)\) initially, for \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})\) the initial point, and thereafter decreases \(\mu\) as \(\mu^{\prime}=\sigma\mu\) at every iteration, while choosing \(\alpha\) to be \(\alpha=0.99\alpha^{\prime}\) for \(\alpha\) the largest \(\alpha>0\) such that \(x_{i}^{\prime}s_{i}^{\prime}>0\), \(w_{j}^{\prime}r_{j}^{\prime}>0\). The two algorithms are summarized in Algorithms 1 and 2.

Algorithm 1 is guaranteed to converge to an \(\epsilon\)-accurate solution in \(\mathcal{O}((n+m)\log(1/\epsilon))\) iterations (Nocedal and Wright, 2006, Theorem 14.3), that is, to a number of iterations proportional to the problem size. Algorithm 2, in contrast, does not come with any theoretical guarantees but is typical of practical IPM algorithms, which tend to converge in an almost constant number of iterations--usually within 30-40 iterations, irrespective of problem size (Gondzio, 2012; Colombo and Gondzio, 2008).

```
0: An LP instance \((\mathbf{A},\mathbf{b},\mathbf{c})\), a barrier reduction hyperparameter \(\sigma\in(0,1)\), a neighborhood hyperparameter \(\gamma\in(0,1]\) and initial values \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})\) such that \(\mathbf{A}\mathbf{x}_{0}-\mathbf{r}_{0}=\mathbf{b}\), \(\mathbf{A}^{\mathsf{T}}\mathbf{w}_{0}+\mathbf{s}_{0}=\mathbf{c}\), \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})>0\) and \(\min_{i}x_{0}w_{0i}\geq\gamma\mu_{0}\), \(\min_{i}w_{0i}r_{0i}\geq\gamma\mu_{0}\) for \(\mu_{0}=(\mathbf{x}_{0}^{\mathsf{T}}\mathbf{s}_{0}+(\mathbf{w}_{0}^{\mathsf{T}})_{0})/(n+m)\).
1:repeat
2:\(\mu\leftarrow(\mathbf{x}^{T}\mathbf{s}+\mathbf{w}^{\mathsf{T}}\mathbf{r})/(n+m)\)
3: Compute \(\Delta\mathbf{w}\) by solving the linear system \(\mathbf{Q}\Delta\mathbf{w}=\mathbf{b}-\mathbf{A}\mathbf{x}+\sigma\mu\mathbf{D}(\mathbf{w})^{-1}\mathbf{1}\) \(+\mathbf{A}\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})[\mathbf{c}-\mathbf{A}^{\mathsf{T}}\mathbf{w}- \sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}]\) for \(\mathbf{Q}=\mathbf{A}\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})\mathbf{A}^{\mathsf{T}}+\mathbf{D}(\mathbf{w} )^{-1}\mathbf{D}(\mathbf{r})\)
4:\(\Delta\mathbf{x}\leftarrow\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})[\mathbf{A}^{\mathsf{T}} \Delta\mathbf{w}-\mathbf{c}+\mathbf{A}^{\mathsf{T}}\mathbf{w}+\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}]\)
5:\(\Delta\mathbf{s}\leftarrow\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}-\mathbf{s}-\mathbf{D}(\mathbf{x})^{ -1}\mathbf{D}(\mathbf{s})\Delta\mathbf{x}\)
6:\(\Delta\mathbf{r}\leftarrow\sigma\mu\mathbf{D}(\mathbf{w})^{-1}\mathbf{1}-\mathbf{r}-\mathbf{D}(\mathbf{w} )^{-1}\mathbf{D}(\mathbf{r})\Delta\mathbf{w}\)
7: Compute the largest \(\alpha\in(0,1)\) such that \[\min_{i,j}\{(\mathbf{x}+\alpha\Delta\mathbf{x})_{i}(\mathbf{s}+\alpha\Delta\mathbf{s})_{i},(\mathbf{w} +\alpha\Delta\mathbf{w})_{j}(\mathbf{r}+\alpha\Delta\mathbf{r})_{j}\}\] \[\geq\gamma\frac{(\mathbf{x}+\alpha\Delta\mathbf{x})^{\mathsf{T}}(\mathbf{s}+ \alpha\Delta\mathbf{s})+(\mathbf{w}+\alpha\Delta\mathbf{w})^{\mathsf{T}}(\mathbf{r}+\alpha \Delta\mathbf{r})}{n+m}.\]
8: Update \((\mathbf{x},\mathbf{w},\mathbf{s},\mathbf{r})\)\(+=\alpha(\Delta\mathbf{x},\Delta\mathbf{w},\Delta\mathbf{s},\Delta\mathbf{r})\)
9:until convergence of \((\mathbf{x},\mathbf{w},\mathbf{s},\mathbf{r})\)
10:return the point \(\mathbf{x}\), which solves 1.
```

**Algorithm 1** Theoretical IPM for LPs

```
0: An LP instance \((\mathbf{A},\mathbf{b},\mathbf{c})\), a barrier reduction hyperparameter \(\sigma\in(0,1)\), a neighborhood hyperparameter \(\gamma\in(0,1]\) and initial values \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})\) such that \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})>0\) and \(\mu_{0}=(\mathbf{x}_{0}^{\mathsf{T}}\mathbf{s}_{0}+\mathbf{w}_{0}^{\mathsf{T}}\mathbf{r}_{0})/(n+m)\).
1:repeat
2: Compute \(\Delta\mathbf{w}\) by solving the linear system \(\mathbf{Q}\Delta\mathbf{w}=\mathbf{b}-\mathbf{A}\mathbf{x}+\sigma\mu\mathbf{D}(\mathbf{w})^{-1}\mathbf{1}\) \(+\mathbf{A}\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})[\mathbf{c}-\mathbf{A}^{\mathsf{T}}\mathbf{w}- \sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}]\) for \(\mathbf{Q}=\mathbf{A}\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})\mathbf{A}^{\mathsf{T}}+\mathbf{D}(\mathbf{w})^{ -1}\mathbf{D}(\mathbf{r})\).
3:\(\Delta\mathbf{x}\leftarrow\mathbf{D}(\mathbf{s})^{-1}\mathbf{D}(\mathbf{x})[\mathbf{A}^{\mathsf{T}} \Delta\mathbf{w}-\mathbf{c}+\mathbf{A}^{\mathsf{T}}\mathbf{w}+\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}]\)
4:\(\Delta\mathbf{s}\leftarrow\sigma\mu\mathbf{D}(\mathbf{x})^{-1}\mathbf{1}-\mathbf{s}-\mathbf{D}(\mathbf{x})^{ -1}\mathbf{D}(\mathbf{s})\Delta\mathbf{x}\)
5:\
features, and the output of \(\mathsf{MSG}^{(t)}\) to a single vectorial representation.

To adapt the parameters of the above functions, they are optimized end-to-end, usually through a variant of stochastic gradient descent, e.g., Kingma and Ba (2015), together with the parameters of a neural network used for classification or regression. In the following, we define a _message-passing step_ as the application of a message and update function.

## 3 Simulating IPMs via MPNNs

We now show that there exist MPNNs, with specific architecture and choices of parameters, such that Algorithms 1 and 2 can be interpreted as inference over these MPNNs for a specific tripartite graph encoding the LP as input.

**Representing LPs as graphs** Let \(I=(\mathbf{A},\mathbf{b},\mathbf{c})\) be an instance of LP. Similar to the setting of Ding et al. (2020), we model the instances with an (undirected) _weighted tripartite graph_\(G(I):=(V(I),C(I),\{o\},E(I)_{\mathrm{vc}},E(I)_{\mathrm{vo}},E(I)_{\mathrm{ co}})\). Here, the node set \(V(I):=\{v_{i}\mid i\in[n]\}\) represents the variables of \(I\), the node set \(C(I):=\{c_{i}\mid i\in[m]\}\) represents the constraints of \(I\), and the node \(o\) represents the objective. Further, the first edge set \(E(I)_{\mathrm{vc}}\) models the variable-constraint interaction, i.e., \(E(I)_{\mathrm{vc}}:=\{(v_{i},c_{j})\mid A_{ij}\neq 0\}\), where each such edge \((v_{i},c_{j})\) is annotated with the weight \(A_{ij}\). Further, the objective node \(o\) is connected to all other nodes in the graphs, i.e., \(E(I)_{\mathrm{vo}}:=\{(o,v_{i})\mid v_{i}\in V(I)\}\) and \(E(I)_{\mathrm{co}}:=\{(o,c_{i})\mid c_{i}\in C(I)\}\). Each edge \((o,v_{i})\in E(I)_{\mathrm{vo}}\) is annotated with the weight \(c_{i}\). Similarly, each edge \((o,c_{i})\in E(I)_{\mathrm{co}}\) is annotated with the weight \(b_{i}\). The resulting graph is illustrated in Figure 0(a).

**Theoretical results** We now state the main results of this paper. To describe them, first notice that Algorithms 1 and 2 operate by taking an initial point \((\mathbf{x}_{0},\mathbf{w}_{0},\mathbf{s}_{0},\mathbf{r}_{0})\) and a duality measure \(\mu_{0}\), and updating them after every iteration, yielding a sequence of points \((\mathbf{x}_{t},\mathbf{w}_{t},\mathbf{s}_{t},\mathbf{r}_{t})\) and duality measure \(\mu_{t}\) for iterations \(t>0\). The following result shows that Algorithm 1 can be reproduced by a specific MPNN, in the sense that a fixed-depth MPNN can reproduce each of its iterations.

**Theorem 1**.: There exists an MPNN \(f_{\mathsf{MPNN},\mathsf{IPM1}}\) composed of \(\mathcal{O}(m)\) message-passing steps that reproduces an iteration of Algorithm 1, in the sense that for any LP instance \(I=(\mathbf{A},\mathbf{b},\mathbf{c})\) and any iteration step \(f_{\mathsf{MPNN},\mathsf{IPM1}}\) maps the graph \(G(I)\) carrying \([\mathbf{x}_{t},\mathbf{s}_{t}]\) on the variable nodes and \([\mathbf{w}_{t},\mathbf{r}_{t}]\) on the constraint nodes on the variable nodes, \([\mathbf{w}_{t+1},\mathbf{s}_{t+1}]\) on the constraint nodes.

This implies, by composing several instances of \(f_{\mathsf{MPNN},\mathsf{IPM1}}\), that Algorithm 1 can be simulated by an MPNN with a number of layers proportional to the number of iterations taken by the algorithm.

**Proposition 2**.: There exists an MPNN \(f_{\mathsf{MPNN},\mathsf{IPM2}}\) composed of \(\mathcal{O}(m)\) message-passing steps that reproduces each iteration of Algorithm 2, in the sense that for any LP instance \(I=(\mathbf{A},\mathbf{b},\mathbf{c})\) and any iteration step \(t\geq 0\), \(f_{\mathsf{MPNN},\mathsf{IPM2}}\) maps the graph \(G(I)\) carrying \([\mathbf{x}_{t},\mathbf{s}_{t}]\) on the variable nodes, \([\mathbf{w}_{t},\mathbf{r}_{t}]\) on the constraint nodes and \([\mu_{t}]\) on the objective node to the same graph \(G(I)\) carrying \([\mathbf{x}_{t+1},\mathbf{s}_{t+1}]\) on the variable nodes, \([\mathbf{w}_{t+1},\mathbf{r}_{t+1}]\) on the constraint nodes and \([\mu_{t+1}]\) on the objective node.

Similarly, this implies, by composing several instances of \(f_{\mathsf{MPNN},\mathsf{IPM2}}\), that Algorithm 2 can be simulated by an MPNN with a number of layers proportional to the number of iterations taken by the algorithm.

**Implication of the theoretical results** Theorem 1 and Proposition 2 show that MPNNs are, in principle, capable of simulating modern IPMs. That is, they are capable of solving LPs to optimality. Hence, our findings shed light on the recent success of MPNN-based neural architectures by Gasse et al. (2019) and similar approaches, which use MPNNs to mimic strong branching within the branch-and-bound framework for solving mixed-integer linear optimization problems. Moreover, in the following section, we derive MPNN architectures that act as lightweight proxies for solving LPs while being able to adapt to the given problem instance distributions. We verify their effectiveness empirically on real-world LP instances stemming from relaxing mixed-integer linear formulations of well-known combinatorial optimization problems.

## 4 IPM-MPNNs: MPNNs FOR LPs

Inspired by the theoretical alignment of IPMs and MPNNs derived above, we outline our IPM-MPNN framework, allowing for solving LP instances while adapting to a given problem instance distribution.

Given an LP instance \(I\), we now outline an asynchronous MPNN operating on the tripartite graph \(G(I)\); see Section 3. Let \(\mathbf{h}_{c}^{(t)}\in\mathbb{R}^{d}\), \(d>0\), be the node features of a constraint node \(c\in C(I)\) at iteration \(t>0\), and let \(\mathbf{h}_{v}^{(t)}\in\mathbb{R}^{d}\) and \(\mathbf{h}_{v}^{(t)}\in\mathbb{R}^{d}\) be the node features of a variable node \(v\in V(I)\) and the objective node \(o\) at iteration \(t\), respectively. Moreover, let \(\mathbf{e}_{co},\mathbf{e}_{vc},\mathbf{e}_{vo}\) denote the edge weights. Initially, at \(t=0\), we set the node features by applying a linear mapping to the raw node features \(\mathbf{x}_{v}\), \(\mathbf{x}_{c}\) or \(\mathbf{x}_{o}\), extracted from
the instance \(I\); see Section 5 for details.

All three node types are updated in three separate update passes. In the first pass, we update the embeddings of constraint nodes from the embeddings of the variable nodes and of the objective node. That is, let \(c\in C(I)\) be a constraint node and let \(t>0\), then

\[\mathbf{h}_{c}^{(t)}\coloneqq \mathsf{UPD}_{c}^{(t)}\Big{[}\mathbf{h}_{c}^{(t-1)},\mathsf{MSG}_{ \textnormal{o}\to\textnormal{c}}^{(t)}\Big{(}\mathbf{h}_{o}^{(t-1)},\mathbf{e}_{oc} \Big{)},\] \[\mathsf{MSG}_{\textnormal{v}\to\textnormal{c}}^{(t)}\Big{(} \big{\{}\!\{\mathbf{h}_{v}^{(t-1)},\mathbf{e}_{vc}\}\mid v\in N(c)\cap V(I)\!\}\Big{)} \Big{]}.\]

Here, the parameterized message function \(\mathsf{MSG}_{\textnormal{v}\to\textnormal{c}}^{(t)}\) maps a multiset of vectors, i.e., variable node features and corresponding edge features \(\mathbf{e}_{vc}\), to a vector in \(\mathbb{R}^{d}\). Similarly, the parameterized function \(\mathsf{MSG}_{\textnormal{o}\to\textnormal{c}}\) maps the current node features of the objective node and edge features \(\mathbf{e}_{oc}\) to a vector in \(\mathbb{R}^{d}\). Finally, the parameterized function \(\mathsf{UPD}_{c}^{(t)}\) maps the constraint node's previous features, the outputs of \(\mathsf{MSG}_{\textnormal{o}\to\textnormal{c}}^{(t)}\) and \(\mathsf{MSG}_{\textnormal{v}\to\textnormal{c}}^{(t)}\) to vector in \(\mathbb{R}^{d}\).

Next, similarly to the above, we update the objective node's features depending on variable and constraint node features,

\[\mathbf{h}_{\textnormal{o}}^{(t)}\coloneqq \mathsf{UPD}_{\textnormal{o}}^{(t)}\Big{[}\mathbf{h}_{o}^{(t-1)}, \mathsf{MSG}_{\textnormal{c}\to\textnormal{o}}^{(t)}\Big{(}\big{\{}\!\{\mathbf{h} _{c}^{(t)},\mathbf{e}_{co}\mid c\in C(I)\!\}\Big{)}\!\Big{\}},\] \[\mathsf{MSG}_{\textnormal{v}\to\textnormal{o}}^{(t)}\Big{(}\big{\{} \!\{\mathbf{h}_{v}^{(t-1)},\mathbf{e}_{vo}\mid v\in V(I)\!\}\Big{)}\Big{]}.\]

Finally, analogously to the update of the constraint nodes' features, we update the representation of a variable node \(v\in V(I)\) from the constraints nodes and objective node,

\[\mathbf{h}_{v}^{(t)}\coloneqq \mathsf{UPD}_{\textnormal{v}}^{(t)}\Big{[}\mathbf{h}_{v}^{(t-1)}, \mathsf{MSG}_{\textnormal{o}\to\textnormal{v}}^{(t)}\Big{(}\mathbf{h}_{o}^{(t)}, \mathbf{e}_{ov}\Big{)},\] \[\mathsf{MSG}_{\textnormal{c}\to\textnormal{v}}^{(t)}\Big{(} \big{\{}\!\{\mathbf{h}_{c}^{(t)},\mathbf{e}_{cv}\mid c\in N(v)\cap V(C)\!\}\big{\}} \Big{)}\Big{]}.\]

The entire process is executed asynchronously, in that the nodes updated later incorporate the most recent features updates from preceding nodes. We map each variable node feature \(\mathbf{h}_{v}^{(t)}\) to \(\mathsf{MLP}(\mathbf{h}_{v}^{(t)})\in\mathbb{R}\), where \(\mathsf{MLP}\) is a multi-layer perceptron, and concatenate the resulting real numbers over all variable nodes column-wise, resulting in the final prediction \(\mathbf{z}^{(t)}\in\mathbb{R}^{n}\).

In the experiments in Section 5, we probe various message-passing layers to express the various message and update functions. Concretely, we leverage the GCN (Kipf and Welling, 2017), GIN (Xu et al., 2019), and GEN (Li et al., 2020) MPNN layers, respectively; see Appendix C for details.

We train the above IPM-MPNN architecture, i.e., adapt its parameters, in a supervised fashion. Below, we outline the three components constituting our training loss function. To that, let \((\mathbf{A},\mathbf{b},\mathbf{c})\) be an LP instance.

**Variable supervision** As discussed above, our IPM-MPNN aims to simulate the solving steps provided by a standard IPM. Thus, aligned with Theorem 1 and Proposition 2, a perfectly parameterized MPNN is expected to follow the steps without deviation for all \(T\) iterations. We maintain the intermediate outputs \(\mathbf{z}^{(t)}\in\mathbb{R}^{n}\) of each MPNN layer and calculate the mean squared error (MSE) loss between every pair of the expert solution \(\mathbf{y}^{(t)}\in\mathbb{R}^{n}\) and MPNN prediction \(\mathbf{z}^{(t)}\). Moreover, we introduce a step decay factor \(\alpha\in[0,1]\) so that early steps play a less important role, resulting in the following loss function,

\[\mathcal{L}_{\textnormal{var}}\coloneqq\frac{1}{N}\sum_{i=1}^{N} \sum_{t=1}^{T}\alpha^{T-t}\|\mathbf{y}_{i}^{(t)}-\mathbf{z}_{i}^{(t)}\|_{2}^{2}, \tag{7}\]

where \(N\) denotes the number of training samples.

**Objective supervision** We use regularization on the prediction regarding the ground-truth objective values at every step. Empirically, this regularization term helps with convergence and helps finding more feasible solutions that minimize the objective in case the LP instance has multiple solutions. Note that we do not predict the objective directly with our MPNNs but calculate it via \(\mathbf{c}^{\mathsf{T}}\mathbf{z}^{(t)}\) instead. Suppose the ground-truth values are given by \(\mathbf{c}^{\mathsf{T}}\mathbf{y}^{(t)}\), we have

\[\mathcal{L}_{\textnormal{obj}}\coloneqq\frac{1}{N}\sum_{i=1}^{N} \sum_{t=1}^{T}\alpha^{T-t}\Big{[}\mathbf{c}^{\mathsf{T}}\big{(}\mathbf{y}_{i}^{(t)}- \mathbf{z}_{i}^{(t)}\big{)}\Big{]}^{2}. \tag{8}\]

**Constraint supervision** Finally, we aim for the IPM-MPNN to predict an optimal solution regarding the objective value while satisfying all the constraints. To that, we introduce a regularization penalizing constraint violations, i.e.,

\[\mathcal{L}_{\textnormal{cons}}\coloneqq\frac{1}{N}\sum_{i=1}^{N} \sum_{t=1}^{T}\alpha^{T-t}\|\mathsf{ReLU}(\mathbf{A}_{i}\mathbf{z}_{i}^{(t)}-\mathbf{b}_{i} )\|_{2}^{2}. \tag{9}\]

Finally, we combine the above three loss terms into the loss function

\[\mathcal{L}\coloneqq w_{\textnormal{var}}\mathcal{L}_{\textnormal{var}}+w_{ \textnormal{obj}}\mathcal{L}_{\textnormal{obj}}+w_{\textnormal{cons}}\mathcal{L }_{\textnormal{cons}}, \tag{10}\]

where we treat \(w_{\textnormal{var}},w_{\textnormal{obj}}\), and \(w_{\textnormal{cons}}>0\) as hyperparamters.

At test time, given an LP instance \(I\), we construct the tripartite graph \(G(I)\) as outlined above and use the trained MPNN to predict the variables' values.

## 5 Experimental Study

Here, we empirically evaluate the ability of IPM-MPNNs to predict the optimal solutions of LPs. In particular, we aim to answer the following questions.
**Q1** Can MPNNs properly imitate the performance of IPM solvers in practice?

**Q2** How is MPNNs' performance compared with competing neural-network-based solutions?

**Q3** What advantage does our MPNN architecture hold compared with traditional IPM solvers?

**Q4** Do MPNNs generalize to instances larger than seen during training?

Our experimental results are reproducible with the code available at [https://github.com/chendiqian/IPM_MPNN](https://github.com/chendiqian/IPM_MPNN).

**Datasets** We obtain LP instances from mixed-integer optimization instances by dropping the integrality constraints over variables. Following Gasse et al. (2019), we use four classes of problems, namely _set covering_, _maximal independent set_, _combinatorial auction_, and _capacitated facility location_. For each problem type, we generate small- and large-size instances. We describe the datasets and our parameters for dataset generation in Appendix D. From each LP instance \(I=(\mathbf{A},\mathbf{b},\mathbf{c})\), we generate a tripartite graph \(G(I)\), following Section 3, and construct initial node and edge features as follows. We indicate the \(i\)th row of the constraint matrix \(\mathbf{A}\) as \(\mathbf{A}_{i}\) and the \(j\)th column as \(\mathbf{A}_{\cdot,j}\). For a given variable node \(v_{j}\in V(I)\), its initial node features are set to the mean and standard deviation of the column vector \(\mathbf{A}_{\cdot,j}\), resulting in two features. Analogously, for a constraint node \(c_{i}\in C(I)\), we derive initial features from the statistical properties of the row vector \(\mathbf{A}_{i}\). For the objective node \(o\), the features are characterized in a corresponding manner using the vector \(\mathbf{c}\). In the supervised learning regime, it is mandatory to have ground truth labels to guide the model predictions. For our IPM-MPNNs, the outputs of each layer \(t\) are the prediction of the variables' value \(\mathbf{z}^{(t)}\). Consequently, we utilize the intermediate variable values \(\mathbf{y}^{(t)}\), as provided by the solver, to serve as our ground truth. However, to prevent the MPNN from becoming excessively deep, we sample the ground truth steps, i.e., we adopt an equidistant sampling strategy for the solver's steps, ensuring that the number of sampled steps aligns with the depth of the corresponding MPNN. We split the graph datasets into train, validation, and test splits with ratios \(0.8,0.1\), and \(0.1\). We conducted the experiments by evaluating each dataset multiple times, performing three independent runs with distinct random seeds. The reported results represent the average numbers over the test set across these runs and the corresponding standard deviations. We executed all experiments on a single NVIDIA A100 80GB GPU card--see Appendix E for training-related hyperparameters.

**IPM-MPNNs' ability to solve LPs (Q1)** We collect and organize the results of our MPNN method in Table 1. We report the numbers on the four types of relaxed MILP instances, each with small and large sizes. As explained in Section 4, we leverage three types of MPNN layers, GCN (Kipf and Welling, 2017), GIN (Xu et al., 2019), and GEN (Li et al., 2020), as the backbone of our MPNN architectures. We split the table into two main parts, namely the mean absolute relative objective gap

\[\frac{1}{N}\sum_{i=1}^{N}\left|\frac{\mathbf{c}^{\mathsf{T}}\big{(}\mathbf{y}_{i}^{(T )}-\mathbf{z}_{i}^{(T)}\big{)}}{\mathbf{c}^{\mathsf{T}}\mathbf{y}_{i}^{(T)}}\right| \times 100\]

of the last step \(T\) over the test set, and the mean absolute constraint violation of the last step

\[\frac{1}{N}\sum_{i=1}^{N}\frac{1}{m_{i}}\|\mathsf{reLU}(\mathbf{A}_{i}\mathbf{z}_{i}^ {(T)}-\mathbf{b}_{i})\|_{1},\]

where the normalization term \(m_{i}\) is the number of constraints of the \(i\)th instance. As seen, our IPM-MPNN architectures consistently align with the IPM solver

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline \multicolumn{2}{c}{} & \multirow{2}{*}{Th.} & \multirow{2}{*}{MPNN} & \multicolumn{2}{c}{**Small instances**} & \multicolumn{4}{c}{**Large instances**} \\ \cline{4-10} \multicolumn{1}{c}{} & & & \multicolumn{2}{c}{**Index**} & \multicolumn{2}{c}{**Case**} & \multicolumn{2}{c}{**Face**} & \multicolumn{2}{c}{**Index**} & \multicolumn{2}{c}{**Case**} & \multicolumn{2}{c}{**Face**} \\ \hline \multirow{4}{*}{**PNN**} & \multirow{4}{*}{✓} & GEN & **0.319\(\pm\)0.020** & 0.119\(\pm\)0.000 & **0.612\(\pm\)0.049** & **0.549\(\pm\)0.121** & 0.629\(\pm\)0.048 & 0.158\(\pm\)0.033 & **0.306\(\pm\)0.047** & **0.747\(\pm\)0.038** \\  & & GCN & 0.418\(\pm\)0.008 & **0.103\(\pm\)0.006** & 0.682\(\pm\)0.029 & 0.578\(\pm\)0.015 & **0.420\(\pm\)0.047** & **0.049\(\pm\)0.005** & 0.407\(\pm\)0.038 & 0.914\(\pm\)0.141 \\  & & GIN & 0.478\(\pm\)0.038 & 0.146\(\pm\)0.011 & 0.632\(\pm\)0.036 & 0.810\(\pm\)0.227 & 0.711\(\pm\)0.15 & 0.126\(\pm\)0.027 & 0.378\(\pm\)0.052 & 0.911\(\pm\)0.127 \\ \cline{2-10}  & & GEN & 8.310\(\pm\)1.209 & 0.735\(\pm\)0.022 & 1.417\(\pm\)0.029 & 2.976\(\pm\)0.013 & 1.517\(\pm\)0.541 & 0.320\(\pm\)0.008 & 0.51\(\pm\)0.122 & 2.531\(\pm\)0.025 \\  & & \multirow{2}{*}{✓} & GEN & 5.523\(\pm\)0.133 & 0.639\(\pm\)0.000 & 1.350\(\pm\)0.001 & 3.031\(\pm\)0.003 & 6.092\(\pm\)0.044 & 0.298\(\pm\)0.000 & 0.766\(\pm\)0.003 & 2.535\(\pm\)0.034 \\  & & GIN & 5.592\(\pm\)0.1179 & 0.634\(\pm\)0.012 & 1.202\(\pm\)0.146 & 2.996\(\pm\)0.031 & 5.835\(\pm\)1.177 & 2.909\(\pm\)0.006 & 0.810\(\pm\)0.140 & 2.660\(\pm\)0.002 \\ \hline \multirow{4}{*}{**PNN**} & \multirow{4}{*}{✓} & GEN & **0.002\(\pm\)0.0002** & 0.000\(\pm\)0.0000 & 0.003\(\pm\)0.0002 & 0.002\(\pm\)0.001 & 0.009\(\pm\)0.001 & 0.001\(\pm\)0.003 & **0.0004\(\pm\)0.0002** & 0.002\(\pm\)0.001 \\  & & \multirow{2}{*}{✓} & GCN & 0.002\(\pm\)0.001 & **0.0003\(\pm\)0.0001** & 0.002\(\pm\)0.0007 & **0.002\(\pm\)0.000** & 0.001\(\pm\)0.000 & **0.0005\(\pm\)0.0004** & 0.001\(\pm\)0.0005 & **0.001\(\pm\)0.0004** \\  & & \multirow{2}{*}{✓} & GCN & 0.004\(\pm\)0.001 & 0.000\(\pm\)0.0001 & 0.001\(\pm\)0.001 & 0.002\(\pm\)0.000 & **0.0008\(\pm\)0.002** & 0.001\(\pm\)0.001 & 0.002\(\pm\)0.0008 & 0.002\(\pm\)0.007 \\ \cline{1-1} \cline{2-10}  & & GEN & 0.181\(\pm\)0.032 & 0.006\(\pm\)0.0003 & 0.006\(\pm\)0.001 & 0.011\(\pm\)0.004 & 0.309\(\pm\)0.022 & 0.004\(\pm\)0.0002 & 0.006\(\pm\)0.001 & 0.003\(\pm\)0.001 \\ \cline{1-1} \cline{2-10}  & & \multirow{2}{*}{✓} & GCN & 0.207\(\pm\)0.006 & 0.004\(\pm\)0.001 & 0.002\(\pm\)0.01 & 0.006\(\pm\)0.001 & 0.267\(\pm\)0.090 & 0.033\(\pm\)0.0004 & 0.004\(\pm\)0.001 & 0.002\(\pm\)0.0003 \\ \cline{1-1}  & & \multirow{2}{*}{✓} & GEN & 0.211\(\pm\)0.0007 & 0.003\(\pm\)0.0002 & 0.003\(\pm\)0.001 & 0.008\(\pm\)0.002 & 0.236\(\pm\)0.014 & 0.033\(\pm\)0.0004
at the last converged step, with marginal constraint violation. Moreover, the relative objective gaps of our method are all under 1%. The GCN-based IPM-MPNNs perform best on large maximal independent set relaxation instances at \(0.094\pm 0.005\%\). Overall, GEN and GCN perform better than GIN layers among the three graph convolutions. The absolute constraint violations are mostly at the \(1\times 10^{-3}\) level, with the best \((0.0003\pm 0.0001)\) achieved by GCN on small maximal independent set relaxation instances.

**Baselines (Q2)** To answer question **Q2**, we compare IPM-MPNNs to two baselines. First, we compare our IPM-MPNNs to Chen et al. (2023), where the authors proposed encoding LP instances as bipartite graphs. In Table 1, the rows with a cross mark (\(\mathbf{\mathcal{X}}\)) on the left are the bipartite baselines. As seen in this table, our IPM-MPNNs outperform the bipartite architectures in all instances with all types of MPNNs. Most relative objective gaps lie above 1%, except on the maximal independent set and large combinatorial auction relaxation instances. The largest gap is observed on the small set covering relaxation instances by GEN, with the baseline reporting as much as 26.1\(\times\) higher constraint violation number. Hence, our results indicate that IPM-MPNNs's tripartite representation is crucial. We also compare our IPM-MPNNs to a neural-ODE-based approach Wu and Lisser (2023). Since their approach is quite expensive during training by means of both time and GPU memory, we generate 1000 mini-sized instances. Due to the architecture-agnostic property of their approach, we embed our MPNNs in their training pipeline. When presenting the runtime and GPU usage, we use the same batch size as the baseline method on our MPNN approach for a fair comparison, even though our method can scale to a much larger batch size in practice. We report results in Table 2. Taking the GEN layer as an example, our method shows consistently better results than the neural-ODE baseline method by reaching at most 14.2\(\times\) lower relative objective gap, 8.6\(\times\) faster training, and 300.8\(\times\) less memory on the mini-sized capacitated facility location relaxation instances.

**Inference time profiling (Q3)** We also compare IPM-MPNN's performance to exact IPM solvers. Thereto, we compare the time required to solve an instance between traditional solvers, namely SciPy's IPM solver and a Python-based custom-build one, and our IPM-MPNN. We run the solvers and our MPNNs on the test set of each dataset and report the mean and standard deviation in seconds. According to Table 3, our MPNN clearly outperforms the SciPy IPM solver on all large instances. Further, IPM-MPNNs beat our Python-based IPM solver described in Algorithm 2 on all instances. It is worth noting that both IPM solvers exhibit sensitivity to problem sizes. For example, the SciPy solver sees a performance degradation of approximately 65.0\(\times\) when transitioning from small to large set covering problems. In contrast, our MPNN method demonstrates a more consistent behavior across varying problem sizes, showing only 1.2\(\times\) slowdown for the analogous instances. Hence, we can positively answer question **Q3**.

**Size generalization (Q4)** We investigate the possibility of generalizing our pre-trained MPNNs to larger instances than encountered during training. To that, we generate new sets of novel instances, each with the same number of instances as the test set. In Table 11 in the appendix, we list our training instance sizes and a collection of test instance sizes. Taking GEN on set covering relaxation problems as an example, as the in

\begin{table}
\begin{tabular}{l l l l l l l} \hline \hline \multicolumn{2}{c}{**Instances**} & SciPy Slower & Our Slower & GEN & GCN & GCN \\ \hline \multirow{4}{*}{\begin{tabular}{l} Small network \\ Large network \\ \end{tabular} } & \multirow{4}{*}{\begin{tabular}{l} ODE \\ \end{tabular} } & GEN & 14.91\(\pm\)0.003 & 6.22\(\pm\)0.007 & 12.35\(\pm\)0.004 & 20.000\(\pm\)0.009 \\  & & GEN & 14.95\(\pm\)0.005 & 6.71\(\pm\)0.005 & 12.94\(\pm\)0.005 & 20.000\(\pm\)0.002 \\  & & GEN & 15.00\(\pm\)0.005 & 6.74\(\pm\)0.004 & 13.74\(\pm\)0.004 & 21.00\(\pm\)0.009 \\  & & GEN & 25.25\(\pm\)0.005 & 1.47\(\pm\)0.005 & 27.90\(\pm\)0.004 & 1.478\(\pm\)0.005 \\  & & GEN & 27.94\(\pm\)0.005 & 1.47\(\pm\)0.005 & 27.60\(\pm\)0.004 & 1.278\(\pm\)0.004 \\ \hline \multirow{4}{*}{\begin{tabular}{l} Small network \\ \end{tabular} } & \multirow{4}{*}{\begin{tabular}{l} ODE \\ \end{tabular} } & GEN & 007.20\(\pm\)0.006 & 0.08\(\pm\)0.005 & 0.09\(\pm\)0.004 & 0.09\(\pm\)0.001 \\  & & GEN & 00.04\(\pm\)0.005 & 0.08\(\pm\)0.005 & 0.09\(\pm\)0.000 & 0.09\(\pm\)0.006 \\  & & GEN & 00.04\(\pm\)0.005 & 0.04\(\pm\)0.005 & 0.09\(\pm\)0.000 & 0.09\(\pm\)0.006 \\  & & GEN & 00.04\(\pm\)0.005 & 0.04\(\pm\)0.005 & 0.09\(\pm\)0.000 & 0.09\(\pm\)0.006 \\ \hline \multirow{4}{*}{\begin{tabular}{l} Small network \\ \end{tabular} } & \multirow{4}{*}{\begin{tabular}{l} Ours \\ \end{tabular} } & GEN & **00.023\(\pm\)0.000** & 0.05\(\pm\)0.000 & 0.01\(\pm\)0.005 & 0.01\(\pm\)0.003 \\  & & GEN & 00.04\(\pm\)0.000 & 0.003\(\pm\)0.000 & 0.09\(\pm\)0.006 & 0.09\(\pm\)0.006 \\  & & GEN & 00.02\(\pm\)0.005 & **0.005\(\pm\)0.000** & **0.014\(\pm\)0.002** & **0.006\(\pm\)0.000** \\ \hline \multirow{4}{*}{\begin{tabular}{l} Small network \\ \end{tabular} } & \multirow{4}{*}{\begin{tabular}{l} ODE \\ \end{tabular} } & GEN & 47.20\(\pm\)0.005 & 51.25\(\pm\)0.003 & 63.00\(\pm\)0.006 & 96.258 \\  & & GEN & 57.19\(\pm\)0.005 & 80.13\(\pm\)0.006 & 73.29\(\pm\)0.006 & 32.297 \\  & & GEN & 55.91\(\pm\)0.005 & 64.028 & 32.90\(\pm\)0.001 & 62.488 \\ \cline{1-1} \cline{2-6}  & & GEN & 10.11\(\pm\)0.077 & 9.617 & 9.366 & 11.21\(\pm\)0.006 \\  & & GEN & 18.94\(\pm\)0.04 & 6.688 & **7.268** & **8.344** \\ \cline{1-1}  & & GEN & **0.042** & **8.005** & 5.817 & 10.771 \\ \hline \multirow{4}{*}{\begin{tabular}{l} Small network \\ \end{tabular} } & \multirow{4}{*}{
\begin{tabular}{l} ODE \\ \end{tabular} } & GEN & 14.65\(\pm\)0.05 & 25.911 & 23.354 & 14.520 \\  & & GEN & 14.469\(\pm\)0.005 & 34.001 & 23.862 & 10.680 \\  & & GEN & 15.289\(\pm\)0.005 & 30.101 & 13.462 & 22.713 \\ \cline{1-1} \cline{2-6}  & & GEN & **0.091** & 0.088 & 0.01 & 0.148 \\ \cline{1-1} \cline{2-6}  & & GEN & 0.201 & 0.131 & **0.090** & **0.142** \\ \cline{1-1} \cline{2-6}  & & GEN & **0.071** & **0.073** & 0.148 & 0.157 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Comparing IPM-MPNNs’ inference time to SciPy’s IPM implementation and our Python-based IPM solver. We report mean and standard deviation in seconds over three runs. We print the best results per target in bold.
ference size grows, the relative objective gap increases a bit, while the constraint violations are overall worse than on the training set. Notably, for the training size [600,800], which is at least \(1.37\times\) larger than the training instances, the objective gap is merely \(0.28\%\) worse than [500,700] size. For question **Q4**, we can conclude that our pre-trained MPNNs can generalize to unseen, larger instances to some extent.

## 6 Conclusion

In summary, our study establishes a strong connection between (MPNNs and IPMs for LPs. We have shown that MPNNs can effectively simulate IPM iterations, revealing their potential in emulating strong branching within the branch-and-bound framework, as demonstrated by Gasse et al. (2019). In addition, building on this connection, we proposed IPM-MPNNs for learning to solve large-scale LP instances approximately, surpassing neural baselines and exact IPM solvers in terms of solving time across various problem domains. Looking forward, promising avenues for further research involve expanding our theoretical framework to encompass a broader range of convex optimization problems.

#### Acknowledgments

CQ and CM are partially funded by a DFG Emmy Noether grant (468502433) and RWTH Junior Principal Investigator Fellowship under Germany's Excellence Strategy.
