\section{Introduction}

Applying machine learning (ML) techniques to accelerate optimization, also known as \emph{Learning to Optimize (L2O)},
is attracting increasing attention recently.
It has been reported in the literature that L2O shows great potentials on both continuous optimization~\cites{monga2021algorithm,chen2021learning,amos2022tutorial} and combinatorial optimization problems~\cites{bengio2021machine,mazyavkina2021reinforcement}
Many of the L2O works train a parameterized model that takes the optimization problem as input and outputs information useful to the classic algorithms, such as a good initial solution and branching decisions~\cite{nair2020solving}, and some even directly generate an approximate optimal solution~\cite{gregor2010learning}.

In these works, one is building an ML model to approximate the mapping from an explicit optimization instance either to its key properties or directly to its solution. The accuracy of this approximation is called the representation power or expressive power of the model. When the approximation is accurate, the model can solve the problem or provide useful information to guide an optimization algorithm.
This paper tries to address a fundamental but open theoretical problem for linear programming (LP):
\begin{equation}
\label{eq:p0}
\tag{P0}
\begin{aligned}
& \textit{Which neural network can represent LP and predict its key properties and solution?}
\end{aligned}
\end{equation}
To clarify, by solution we mean the optimal solution. We hope the answer to this question paves the way toward answering this question for other optimization types.

\paragraph{\textbf{Linear Programming (LP)}} LP is an important type of optimization problem with a long history and a wide range of applications, such as scheduling~\cite{hanssmann1960linear}, signal processing~\cite{candes2005decoding}, machine learning~\cite{dedieu2022solving}, etc. A general LP problem is defined as:
\begin{equation}\label{MILP}
	\begin{split}
		\min_{x\in\mathbb{R}^n} ~~ c^\top x,\quad
		\textup{s.t.} ~~ Ax\circ b,\ l\leq x\leq u
	\end{split}
\end{equation}
where $A\in\mathbb{R}^{m\times n}$, $c\in\mathbb{R}^n$, $b\in\mathbb{R}^m$, $l\in(\mathbb{R}\cup\{-\infty\})^n$, $u\in(\mathbb{R}\cup\{+\infty\})^n$, and $\circ\in\{\leq, = ,\geq\}^m$.
Any LP problems must follow one of the following three cases~\cite{bertsimas1997introduction}:
\begin{itemize}[leftmargin=*]
    \item Infeasible. The feasible set $\mathcal{X}_{F} := \{x\in \mathbb{R}^n: Ax\circ b,\ l\leq x\leq u\}$ is empty. In another word, there is no point in $\mathbb{R}^n$ that satisfies the constraints in LP \eqref{MILP}.
    \item Unbounded. The feasible set is non-empty, but the objective value can be arbitrarily good, i.e., lower unbounded. For any $R>0$, there exists an $x\in\mathcal{X}_{F}$ such that $c^\top x < -R$.
    \item Feasible and bounded. There exists $x^\ast \in \mathcal{X}_{F}$ such that $c^\top x^\ast \leq c^\top x$ for all $x \in \mathcal{X}_{F}$. Such $x^\ast$ is named as an optimal solution, and $c^\top x^\ast$ is the optimal objective value.
\end{itemize}

Thus, considering \eqref{eq:p0} on the LP problem \eqref{MILP}, an ideal ML model is expected to be able to predict the three key characteristics of LP: \emph{feasibility, boundedness, and one of its optimal solutions} (if exist), by taking the LP features $(A,b,c,l,u,\circ)$ as input.

Such input has strong mathematical structure.
If we swap the position of the $i,j$-th variable in \eqref{MILP}, elements in vectors $b,c,l,u,\circ$ and columns of matrix $A$ will be reordered. The reordered features $(\hat{A},b,\hat{c},\hat{l},\hat{u},\hat{\circ})$ actually represent an exact equivalent LP problem with the original one $(A,b,c,l,u,\circ)$. Such property is named as \emph{permutation invariance}.
If we do not explicitly restrict ML models with a permutation invariant structure, the models may overfit the variable/constraint orders of instances in the training set.
Motivated by this point, we adopt \emph{Graph Neural Networks (GNNs)} that are permutation invariant naturally.

\paragraph{\textbf{GNN in L2O}}GNN is a type of neural network defined on graphs and widely applied in many areas, for example, recommender systems, traffic, chemistry, etc \cites{wu2020comprehensive,zhou2020graph}. Solving optimization problems with GNN is a recent arising topic in the literature \cites{peng2021graph,ijcai2021-595}.
Many graph-related optimization problems, like minimum vertex cover, traveling salesman problem, vehicle routing problem, can be represented and solved approximately with GNN due to their problem structure \cites{khalil2017learning,kool2018attention,joshi2019efficient,drori2020learning}. Besides that, one may solve a general LP or mixed-integer linear programming (MILP) with the help of GNN.
\cite{gasse2019exact} proposed to represent an MILP problem with a bipartite graph, then applied a message-passing GNN\footnote{We will abbreviate message-passing GNN as GNN in this paper for simplicity.} on this graph to guide a classic MILP solver.
Since that, many approaches have been proposed to guide MILP or LP solvers with GNN \cites{nair2020solving,gupta2020hybrid,ding2020accelerating,shen2021learning,gupta2022lookback,khalil2022mip,liu2022learning,paulus2022learning,qu2022improved,li2022learning}.
Although encouraging empirical results have been observed, theoretical foundations are still lack for this approach.
Specifying \eqref{eq:p0}, we ask:
\begin{equation}
\label{eq:p1}
\tag{P1}
   \begin{aligned}
&\textit{Are there GNNs that can predict the feasibility, boundedness,}\\
&\textit{and an optimal solution (if exists), for all LPs?}
\end{aligned}
\end{equation}

\paragraph{\textbf{Related works and contributions}} To answer \eqref{eq:p1}, one needs the theory of \emph{separation power} and \emph{representation power}.
Separation power of a neural network (NN) means its ability to distinguish two different inputs.
In our settings, a NN with strong separation power means that it can outputs different results when it is applied on any two different LPs.
Representation power of NN means its ability to approximate functions of interest. A famous example of representation power is the Universal Approximation Theorem \cites{hornik1989multilayer,cybenko1989approximation}: multilayer perceptron (MLP), a specific type of NN, can universally approximate continuous functions on a compact domain.
The theory of representation power is established upon the theory of separation power. Only functions with strong enough separation power may possess strong representation power.
The power of GNN has been studied in the literature (see \cites{sato2020survey,jegelka2022theory,Li2022} for comprehensive surveys), and some theoretical efforts have been made to represent some graph-related optimization problems with GNN~\cites{sato2019approximation,Loukas2020What}. However, there are still gaps to answer question \eqref{eq:p1} since the relationships between characteristics of LP and properties of graphs are not well established.
Our contributions are listed below:
\begin{itemize}[leftmargin=*]
    \item (Separation Power). In the literature, it has been shown that the separation power of GNN is equal to the WL test\cites{xu2019powerful,azizian2020expressive,geerts2022expressiveness}. However, there exist many pairs of LPs that cannot be distinguished by the WL test. We show that those puzzling LP pairs share the same feasibility, boundedness, and even an optimal solution if exists. Thus, GNN has strong enough separation power.
    \item (Representation Power). To the best of our knowledge, we established the first complete proof
    that GNN can universally represent a broad class of LPs.
    More precisely, we prove that, there exist GNNs that can be arbitrarily close to the following three mappings:
    \begin{itemize}
        \item LP $\to$ feasibility,
        \item LP $\to$ optimal objective value ($-\infty$ if unbounded and $\infty$ if infeasible),
        \item LP $\to$ an optimal solution (if exists),
    \end{itemize}
    although they are not continuous functions and cannot be covered by the literature\cites{keriven2019universal,chen2019equivalence,maron2019provably,maron2019universality,keriven2021universality}.
    \item (Experimental Validation). We design and conduct experiments that demonstrate the power of GNN on representing LP.
\end{itemize}

The rest of this paper is organized as follows. In Section \ref{sec:preliminary}, we provide preliminaries, including related notions, definitions and concepts. In Section \ref{sec:mainthm}, we show our main theoretical results. The sketches of proofs are provided in Section\ref{sec:sketch}.  We validate our results with numerical experiments in Section \ref{sec:numerical} and we conclude this paper with Section \ref{sec:conclusions}.

\section{Preliminaries}
\label{sec:preliminary}
In this section, we present concepts and definitions that will be used throughout this paper. We first describe how to represent an LP with a weighted bipartite graph, then we define GNN on those LP-induced graphs, and finally we further clarify question (\ref{eq:p1}) with strict mathematical definitions.

\subsection{LP represented as weighted bipartite graph}
\label{sec:lp-graph}
Before representing LPs with graphs, we first define the graph that we will adopt in this paper: \emph{weighted bipartite graph}.
A weighted bipartite graph $G = (V\cup W,E)$ consists of a vertex set $V\cup W$ that are divided into two groups $V$ and $W$ with $V\cap W = \emptyset$, and a collection $E$ of weighted edges, where
each edge connects exactly one vertex in $V$ and one vertex in $W$. Note that there is no edge connecting vertices in the same vertex group. $E$ can also be viewed as a function $E:V\times W\rightarrow \mathbb{R}$. We use $\mathcal{G}_{m,n}$ to denote the collection of all weighted bipartite graphs $G = (V\cup W,E)$ with $|V| = m$ and $|W| = n$. We always write $V = \{v_1,v_2,\dots,v_m\}$, $W = \{w_1,w_2,\dots,w_n\}$, and $E_{i,j} = E(v_i,w_j)$, for $i\in\{1,2,\dots,m\},\ j\in\{1,2,\dots,n\}$.

One can equip each vertex with a feature vector. Throughout this paper, we denote $h_i^V\in \mathcal{H}^V$ as the feature vector of
vertex $v_i\in V$ and denote $h_j^W\in\mathcal{H}^W$ as the feature vector of
vertex $w_j\in W$, where $\mathcal{H}^V,\mathcal{H}^W$ are feature spaces.
Then we define $\mathcal{H}^V_m := (\mathcal{H}^V)^m, \mathcal{H}^W_n := (\mathcal{H}^W)^n$ and concatenate all the vertex features together as $H = (h_1^V,h_2^V,\dots,h_m^V,h_1^W,h_2^W,\dots,h_n^W) \in \mathcal{H}^V_m\times\mathcal{H}^W_n$. Finally, a weighted bipartite graph with vertex features is defined as a tuple $(G,H)\in \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times\mathcal{H}^W_n$.

With the concepts described above, one can represent an LP \eqref{MILP} as a graph~\cite{gasse2019exact}. Each vertex in $W$ represents a variable in LP and each vertex in $V$ represents a constraint. The graph topology is defined with the matrix $A$ in the linear constraint.  More specifically, let us set
\begin{itemize}[leftmargin=*]
    \item Vertex $v_i$ represents the $i$-th constraint in $Ax\circ b$, and vertex $w_j$ represents the $j$-th variable $x_j$.
    \item Information of constraints is involved in the feature of $v_i$: $h_i^{V} = (b_i,\circ_i)$.
    \item The space of constraint features is defined as $\mathcal{H}^V := \mathbb{R}\times \{\leq,=,\geq\}$.
    \item Information of variables is involved in the feature of $w_i$: $h_j^{W} = (c_j,l_j,u_j)$.
    \item The space of variable features is defined as $\mathcal{H}^W := \mathbb{R}\times (\mathbb{R}\cup\{-\infty\})\times (\mathbb{R}\cup\{+\infty\})$.
    \item The edge connecting $v_i$ and $w_j$ has weight $E_{i,j} = A_{i,j}$.
\end{itemize}
Then an LP is represented as a graph $(G,H)\in \mathcal{G}_{m,n}\times \mathcal{H}^V_m\times\mathcal{H}^W_n$. In the rest of this paper, we coin such graphs as \emph{LP-induced Graphs} or \emph{LP-Graphs} for simplicity. We present an LP instance and its corresponding LP-graph in Figure~\ref{fig:ex_LPgraph}.

\begin{figure}[htb!]
	\centering
	\begin{minipage}{0.35\textwidth}
\begin{equation*}
	\begin{split}
		\min_{x\in\mathbb{R}^3} ~~ & x_1 + x_2 + 2 x_3,\\
		\text{s.t.} ~~ & x_1 + 2 x_2 + x_3 \geq 1,\\
		& 2 x_1 + x_2 + x_3 \geq 1, \\
		&  x_1\geq -2,\\
		&  x_2\geq -2,\\
		&  -3\leq x_3\leq 3.
	\end{split}
\end{equation*}
\end{minipage}
\begin{minipage}{0.6\textwidth}
	\begin{tikzpicture}[
		constraintb/.style={circle, draw = blue},
		variablex/.style={circle, draw = red},
		]

		\draw (-1,0.7) node[constraintb] (b1) {$v_1$};
		\draw (-1,-0.7) node[constraintb] (b2) {$v_2$};
		\draw (1,1.4) node[variablex] (x1) {$w_1$};
		\draw (1,0) node[variablex] (x2) {$w_2$};
		\draw (1,-1.4) node[variablex] (x3) {$w_3$};
		\draw (-2.5,0.7) node {$h_1^V = (1,\geq)$};
		\draw (-2.5,-0.7) node {$h_2^V = (1,\geq)$};
		\draw (3,1.4) node {$h_1^W = (1,-2,+\infty)$};
		\draw (3,0) node {$h_2^W = (1,-2,+\infty)$};
		\draw (3,-1.4) node {$h_3^W = (2,-3,3)$};

		\draw[-] (x1.west) -- (b1.east) node[pos = 1/2] {$\textcolor{red}{1}$};
		\draw[-] (x2.west) -- (b1.east) node[pos = 1/3] {$\textcolor{red}{2}$};
		\draw[-] (x3.west) -- (b1.east) node[pos = 1/3] {$\textcolor{red}{1}$};
		\draw[-] (x1.west) -- (b2.east) node[pos = 1/3] {$\textcolor{red}{2}$};
		\draw[-] (x2.west) -- (b2.east) node[pos = 1/3] {$\textcolor{red}{1}$};
		\draw[-] (x3.west) -- (b2.east) node[pos = 1/2] {$\textcolor{red}{1}$};

	\end{tikzpicture}
\end{minipage}
\caption{An example of LP-graph}
\label{fig:ex_LPgraph}
\end{figure}

\subsection{Graph neural networks for LP}
Now we define GNNs on LP-graphs. The GNNs in this paper always take an LP-Graph as input and the output has two cases:
\begin{itemize}[leftmargin=*]
    \item The output is a single real number. In this case, GNN is a function $\mathcal{G}_{m,n}\times \mathcal{H}^V_m\times\mathcal{H}^W_n \to \mathbb{R}$ and usually used to predict the properties of the whole graph.
    \item Each vertex in $W$ has an output. Consequently, GNN is a function $\mathcal{G}_{m,n}\times \mathcal{H}^V_m\times\mathcal{H}^W_n \to \mathbb{R}^n$. Since $W$ represents variables in LP, GNN is used to predict properties of each variable in this case.
\end{itemize}
Now we define the GNN structure precisely. First we encode the input features into the embedding space with learnable functions $f_{\mathrm{in}}^V:\mathcal{H}^V\to \mathbb{R}^{d_0}$ and $f_{\mathrm{in}}^W:\mathcal{H}^W\to \mathbb{R}^{d_0}$
\begin{equation}
    \label{eq:gnn-in}
    h_i^{0,V} = f_{\mathrm{in}}^V(h_i^V),~~h_j^{0,W} = f_{\mathrm{in}}^W(h_j^W),~~i=1,2,\dots,m,~j=1,2,\dots,n,
\end{equation}
where $h_i^{0,V},h_j^{0,W} \in \mathbb{R}^{d_{0}}$ are initial embedded vertex features, $d_{0}$ is their dimension.
Then we choose learnable functions $f_l^V,f_l^W:\mathbb{R}^{d_{l-1}} \rightarrow \mathbb{R}^{d_{l}}$ and $g_l^V,g_l^W: \mathbb{R}^{d_{l-1}} \times \mathbb{R}^{d_{l}} \rightarrow \mathbb{R}^{d_{l}}$ and update the hidden states with\footnote{Note that the update rules in \eqref{eq:gnn-update-v} and \eqref{eq:gnn-update-w} follow a message-passing way, where each vertex only collects information from its neighbors. Since $E_{i,j}=0$ if there is no connection between vertices $v_i$ and $w_j$, the sum operator in \eqref{eq:gnn-update-v} can be rewritten as $\sum_{j \in \mathcal{N}(v_i)}$, where $\mathcal{N}(v_i)$ denotes the set of neighbors of vertex $v_i$.}:
\begin{align}
    h_i^{l,V} = g_l^V\bigg( h_i^{l-1,V},\sum_{j=1}^n E_{i,j} f_l^W(h_j^{l-1, W}) \bigg),\quad i=1,2,\dots,m,\label{eq:gnn-update-v}\\
    h_j^{l,W} = g_l^W\bigg( h_j^{l-1,W},\sum_{i=1}^m E_{i,j} f_l^V(h_i^{l-1, V}) \bigg),\quad j=1,2,\dots,n,\label{eq:gnn-update-w}
\end{align}
where $h_i^{l,V},h_j^{l,W} \in \mathbb{R}^{d_{l}}$ are vertex features are layer $l (1 \leq l \leq L)$ and their dimensions are $d_{1},\cdots,d_{L}$ respectively.
The output layer of the single-output GNN is defined with learnable function $f_\text{out}: \mathbb{R}^{d_{L}} \times \mathbb{R}^{d_{L}}\to\mathbb{R}$
\begin{equation}
\label{eq:gnn-out-single}
	y_{\mathrm{out}} = f_\text{out}\bigg(\sum_{i=1}^m h_i^{L,V},\sum_{j=1}^n h_j^{L,W}\bigg).
\end{equation}
The output of the vertex-output GNN is defined with learnable function $f^W_{\text{out}}: \mathbb{R}^{d_{L}} \times\mathbb{R}^{d_{L}} \times \mathbb{R}^{d_{L}} \rightarrow\mathbb{R}$
\begin{equation}
\label{eq:gnn-out-vertex}
y_{\mathrm{out}}(w_j) = f^W_\text{out}\bigg(\sum_{i=1}^m h_i^{L,V},\sum_{j=1}^n h_j^{L,W},h_j^{L,W}\bigg),\quad j = 1,2,\cdots,n.
\end{equation}
We denote collections of single-output GNNs and vertex-output GNNs with $\mathcal{F}_{\text{GNN}}$ and $\mathcal{F}_{\text{GNN}}^W$ respectively.
\begin{equation}
    \label{eq:gnn-sets}
    \begin{aligned}
    \mathcal{F}_{\text{GNN}} = & \{F:\mathcal{G}_{m,n}\times \mathcal{H}^V_m\times\mathcal{H}^W_n \to \mathbb{R}\ |\ F \text{ yields } \eqref{eq:gnn-in},\eqref{eq:gnn-update-v},\eqref{eq:gnn-update-w},\eqref{eq:gnn-out-single} \}, \\
    \mathcal{F}_{\text{GNN}}^W = & \{F:\mathcal{G}_{m,n}\times \mathcal{H}^V_m\times\mathcal{H}^W_n \to \mathbb{R}^n\ |\ F \text{ yields } \eqref{eq:gnn-in},\eqref{eq:gnn-update-v},\eqref{eq:gnn-update-w},\eqref{eq:gnn-out-vertex} \}.
    \end{aligned}
\end{equation}
In practice, all the learnable functions in GNN $f_{\mathrm{in}}^V,f_{\mathrm{in}}^W,f_\text{out}, f^W_{\text{out}}, \{f_l^V,f_l^W,g_l^V,g_l^W\}_{l=0}^{L}$ are usually parameterized with multi-linear perceptrons (MLP). In our theoretical analysis, we assume for simplicity that those functions may take all continuous functions on given domains, following the settings in \cite{azizian2020expressive}*{Section C.1}. Thanks to the universal approximation properties of MLP \cites{hornik1989multilayer,cybenko1989approximation}, one can extend our theoretical results by taking those learnable functions as large enough MLPs.

\subsection{Revisiting question (\ref{eq:p1})}
\label{sec:p1-math-define}

With the definitions of LP-Graph and GNN defined above, we revisit question (\ref{eq:p1}) and give a more precise mathematical description here. First we define three mappings that respectively describe the feasibility, optimal objective value and an optimal solution of an LP.
\paragraph{\textbf{Feasibility mapping}} The feasibility mapping is a classification function
\begin{equation}\label{feas_map}
    \Phi_{\text{feas}}: \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n \rightarrow \{0,1\},
\end{equation}
where $\Phi_{\text{feas}}(G,H) = 1$ if the LP associated with $(G,H)$ is feasible and $\Phi_{\text{feas}}(G,H) = 0$ otherwise.

\paragraph{\textbf{Optimal objective value mapping}} Denote
\begin{equation}\label{obj_map}
    \Phi_{\text{obj}}: \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n \rightarrow \mathbb{R}\cup\{\infty,-\infty\},
\end{equation}
as the optimal objective value mapping, i.e., for any $(G,H)\in \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$, $\Phi_{\text{obj}}(G,H)$ is the optimal objective value of the LP problem associated with $(G,H)$.

\begin{remark}
The optimal objective value of a LP problem can be a real number or $\infty/-\infty$. The ``$\infty$'' case corresponds to infeasible problems, while the ``$-\infty$'' case consists of LP problems whose objective function is unbounded from below in the feasible region. The preimage of all finite real numbers under $\Phi_{\text{obj}}$, $\Phi_{\text{obj}}^{-1}(\mathbb{R})$, actually describes all LPs with finite optimal objective value.
\end{remark}

\begin{remark}\label{rmk:soluLeast2norm}
In the case that a LP problem has finite optimal objective value, it is possible that the problem admits multiple optimal solutions. However, the optimal solution with the smallest $\ell_2$-norm must be unique. In fact, if $x\neq x'$ are two different solutions with $\|x\| = \|x'\|$. Then $\frac{1}{2}(x+x')$ is also an optimal solution due to the convexity of LPs, and it holds that $\|\frac{1}{2}(x+x')\|^2 < \frac{1}{2}\|x\|^2 + \frac{1}{2}\|x'\|^2 = \|x\|^2 = \|x'\|^2$, i.e., $\|\frac{1}{2}(x+x')\| < \|x\| = \|x'\|$, where the inequality is strict since $x\neq x'$. Therefore, $x$ and $x'$ cannot be optimal solutions with the smallest $\ell_2$-norm.
\end{remark}

\paragraph{\textbf{Optimal solution mapping}} For any $(G,H)\in \Phi_{\text{obj}}^{-1}(\mathbb{R})$, we have remarked before that the LP problem associated with $(G,H)$ has a unique optimal solution with the smallest $\ell_2$-norm. Let
\begin{equation}\label{solu_map}
    \Phi_{\text{solu}}: \Phi_{\text{obj}}^{-1}(\mathbb{R})\rightarrow \mathbb{R}^n,
\end{equation}
be the mapping that maps $(G,H)\in \Phi_{\text{obj}}^{-1}(\mathbb{R})$ to the optimal solution with the smallest $\ell_2$-norm.

\paragraph{\textbf{Invariance and Equivariance}}
We denote $S_m,S_n$ as the group consisting of all permutations on vertex groups $V,W$ respectively. In another word, $S_m$ involves all permutations on the constraints of LP and $S_n$ involves all permutations on the variables. In this paper, we say a function $F:\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n \to \mathbb{R}$ is invariant if it satisfies
\begin{equation*}
    F(G,H) = F\left(  (\sigma_V,\sigma_W)\ast (G,H)  \right),~~\forall \sigma_V \in S_m, \sigma_W \in S_n,
\end{equation*}
and a function $F_W:\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n \to \mathbb{R}^n$ is equivariant if it satisfies
\begin{equation*}
    \sigma_W( F_W(G,H)) = F_W \left(  (\sigma_V, \sigma_W)\ast (G,H))  \right),~~\forall \sigma_V \in S_m, \sigma_W \in S_n,
\end{equation*}
where $(\sigma_V, \sigma_W)\ast (G,H)$ is the permuted graph obtained from reordering indices in $(G,H)$ using $(\sigma_V,\sigma_W)$, which is the group action $S_m\times S_n$ on $\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$. One can check that $\Phi_{\text{feas}}$ and $\Phi_{\text{obj}}$ and any $F\in\mathcal{F}_{\text{GNN}}$ are invariant, $\Phi_{\text{solu}}$ and any $F_W \in \mathcal{F}_{\text{GNN}}^W$ are equivariant.

Question \eqref{eq:p1} actually asks: are there GNNs that can accurately approximate $\Phi_{\text{feas}}$, $\Phi_{\text{obj}}$ and $\Phi_{\text{solu}}$? More precisely, does there exist $F \in \mathcal{F}_{\mathrm{GNN}}$ that is arbitrarily close to $\Phi_{\text{feas}}$ and $\Phi_{\text{obj}}$? And does there exist $F_W \in \mathcal{F}_{\mathrm{GNN}}^W$ that is arbitrarily close to $\Phi_{\text{solu}}$?

\section{Main results}
\label{sec:mainthm}

This section presents our main theorems that answer question \eqref{eq:p1}. As we state in the introduction, representation power is built upon separation power in our paper. We first present with the following theorem that GNN has strong enough separation power to represent LP.
\begin{theorem}\label{thm:sameGNN2sameFeas}
    For any two LP problems respectively corresponding to $(G,H), (\hat{G},\hat{H})\in \mathcal{G}_{m,n}\times\mathcal{H}^V_m \times\mathcal{H}^W_n$, if $F(G,H) = F(\hat{G},\hat{H}),\ \forall~F\in\mathcal{F}_{\text{GNN}}$, then they share some common characteristics:
    \begin{itemize}[leftmargin=8mm]
        \item[(i)] Both LP problems are feasible or both are infeasible, i.e., $\Phi_{\text{feas}}(G,H) = \Phi_{\text{feas}}(\hat{G},\hat{H})$.
        \item[(ii)] The two LP problems have the same optimal objective value, i.e., $\Phi_{\text{obj}}(G,H) = \Phi_{\text{obj}}(\hat{G},\hat{H})$.
        \item[(iii)] If both problems are feasible with finite optimal objective value, then they have the same optimal solution with the smallest $\ell_2$-norm, up to a permutation.
    \end{itemize}
    Furthermore, if $F_W(G,H) = F_W(\hat{G},\hat{H}),\ \forall~F_W\in\mathcal{F}_{\text{GNN}}^W$, then (iii) holds without taking permutations, i.e., $\Phi_{\text{solu}}(G,H) = \Phi_{\text{solu}}(\hat{G},\hat{H})$.
\end{theorem}
This theorem demonstrates that the function spaces $\mathcal{F}_{\mathrm{GNN}}$ and $\mathcal{F}_{\mathrm{GNN}}^W$ are rich enough to distinguish the characteristics of LP. Given two LP instances $(G,H), (\hat{G},\hat{H})$, as long as their feasibility or boundedness are different, there must exist $F \in \mathcal{F}_{\mathrm{GNN}}$ that can distinguish them: $F(G,H)\neq F(\hat{G},\hat{H})$. Moreover, as long as their optimal solutions with the smallest $\ell_2$-norm are different, there must exist $F_W \in \mathcal{F}_{\mathrm{GNN}}^W$ that can distinguish them: $F_W(G,H)\neq F_W(\hat{G},\hat{H})$.

With Theorem \ref{thm:sameGNN2sameFeas} served as a foundation, we can prove that GNN universally approximates the three mappings $\Phi_{\text{feas}}$, $\Phi_{\text{obj}}$ and $\Phi_{\text{solu}}$. Before presenting those results, we first define some concepts of the space $\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$.

\paragraph{\textbf{Topology and measure}} Let us define the topology and the measure on $\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$, where $\mathcal{H}^V = \mathbb{R} \times\{\leq, =, \geq\}$ and $\mathcal{H}^W = \mathbb{R} \times (\mathbb{R}\cup\{-\infty\})\times (\mathbb{R}\cup\{+\infty\})$. Throughout this paper, we consider $\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times\mathcal{H}^W_n$ as a topology space with product topology and measurable space with product measure. It's enough to define the topology and measure of each part separately. Since  each graph in this paper (without vertex features) can be represented with matrix $A\in\mathbb{R}^{m \times n}$, the graph space $\mathcal{G}_{m,n}$ is isomorphic to the Euclidean space $\mathbb{R}^{m \times n}$: $\mathcal{G}_{m,n} \cong \mathbb{R}^{m\times n}$ and we equip $\mathcal{G}_{m,n}$ with the standard Euclidean topology and the standard Lebesgue measure. The real spaces $\mathbb{R}$ in $\mathcal{H}^W$ and $\mathcal{H}^V$ are also equipped with the standard Euclidean topology and Lebesgue measure.
All the discrete spaces $\{\leq, =, \geq\}$, $\{-\infty\}$, and $\{+\infty\}$ have the discrete topology, and all unions are disjoint unions.
We equip those spaces with a discrete measure $\mu(S) = |S|$, where $|S|$ is the number of elements in a discrete space $S$.
This finishes the whole definition and we denote $\text{Meas}(\cdot)$ as the measure on $\mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$.

\begin{theorem}\label{thm:GNNapproxPhi_feas}
    Let $X\subset \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$ be measurable with finite measure. For any $\epsilon>0$, there exists some $F\in\mathcal{F}_{\text{GNN}}$, such that
	\begin{equation*}
	\label{eq:GNNapproxPhi_feas}
		\text{Meas}\left(\left\{(G,H)\in X:\mathbb{I}_{F(G,H)>1/2}\neq \Phi_{\text{feas}}(G,H)\right\}\right)<\epsilon,
	\end{equation*}
    where $\mathbb{I}_{\cdot}$ is the indicator function, i.e., $\mathbb{I}_{F(G,H)>1/2} = 1$ if $F(G,H)>1/2$ and $\mathbb{I}_{F(G,H)>1/2} = 0$ otherwise.
\end{theorem}
This theorem shows that GNN is a good classifier for LP instances in $X$ as long as $X$ has finite measure. If we use $F(G,H)>1/2$ as the criteria to predict the feasibility, the classification error rate is controlled by $\epsilon / \text{Meas}(X)$, where $\epsilon$ can be arbitrarily small. Furthermore, we show that GNN can perfectly fit any dataset with finite samples, which is presented in the following corollary.

\begin{corollary}\label{cor:GNNapproxPhi_feas}
    For any $\mathcal{D}\subset \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$ with finite instances, there exists $F\in\mathcal{F}_{\text{GNN}}$ that
    \begin{equation*}
        \mathbb{I}_{F(G,H)>1/2} = \Phi_{\text{feas}}(G,H),\quad\forall~(G,H)\in \mathcal{D}.
    \end{equation*}
\end{corollary}

Besides feasibility, GNN can also approximate the optimal objective value and the optimal solution with the smallest $\ell_2$-norm up to any given precision. The following theorems provide formal statements.

\begin{theorem}\label{thm:GNNapproxPhi_obj}
    Let $X\subset \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$ be measurable with finite measure. For any $\epsilon>0$, there exists $F_1\in\mathcal{F}_{\text{GNN}}$ such that
        \begin{equation}
        \label{eq:GNNapproxPhi_obj_1}
            \text{Meas}\left(\left\{(G,H)\in X:\mathbb{I}_{F_1(G,H)>1/2}\neq \mathbb{I}_{\Phi_{\text{obj}}(G,H) \in \mathbb{R}}\right\}\right)<\epsilon.
        \end{equation}
 For any $\epsilon,\delta>0$, there exists $F_2\in\mathcal{F}_{\text{GNN}}$ such that
	    \begin{equation}
	    \label{eq:GNNapproxPhi_obj_2}
		    \text{Meas}\left(\left\{(G,H)\in X \cap \Phi_{\text{obj}}^{-1}(\mathbb{R}):| F_2(G,H) -  \Phi_{\text{obj}}(G,H)|>\delta\right\}\right)<\epsilon.
    	\end{equation}
\end{theorem}
Recall the definition of $\Phi_{\text{obj}}$ in \eqref{obj_map} that it can take $\{\pm \infty\}$ as its value. Thus, $\Phi_{\text{obj}}(G,H) \in \mathbb{R}$ means the LP corresponding to $(G,H)$ is feasible and bounded with a finite optimal objective value, and inequality \eqref{eq:GNNapproxPhi_obj_1} illustrates that GNN can identify those feasible and bounded LPs among the whole set $X$, up to a given precision $\epsilon$. Inequality \eqref{eq:GNNapproxPhi_obj_2} shows that GNN can also approximate the optimal value. The measure of the set of LP instances of which the optimal value cannot be approximated with $\delta$-precision is controlled by $\epsilon$.
The following corollary gives the results on dataset with finite instances.

\begin{corollary}\label{cor:GNNapproxPhi_obj}
    Let $\mathcal{D}\subset \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$ be a finite dataset. There exists $F_1\in\mathcal{F}_{\text{GNN}}$ such that
        \begin{equation*}
            \mathbb{I}_{F_1(G,H)>1/2} = \mathbb{I}_{\Phi_{\text{obj}}(G,H) \in \mathbb{R}},\quad\forall~(G,H)\in \mathcal{D}.
        \end{equation*}
   For any $\delta>0$, there exists $F_2\in\mathcal{F}_{\text{GNN}}$, such that
	    \begin{equation*}
		    | F_2(G,H) -  \Phi_{\text{obj}}(G,H)| < \delta,\quad\forall~(G,H)\in \mathcal{D}\cap \Phi_{\text{obj}}^{-1}(\mathbb{R}).
    	\end{equation*}
\end{corollary}

Finally, we show that GNN is able to represent the solution mapping $\Phi_{\text{solu}}$.
\begin{theorem}\label{thm:GNNapproxPhi_solu}
    Let $X\subset \Phi_{\text{obj}}^{-1}(\mathbb{R})\subset \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$ be measurable with finite measure and be closed under actions in $S_m\times S_n$. For any $\epsilon,\delta>0$, there exists some $F_W\in\mathcal{F}_{\text{GNN}}^W$, such that
	\begin{equation*}
		\text{Meas}\left(\left\{(G,H)\in X:\| F(G,H) -  \Phi_{\text{solu}}(G,H)\|>\delta\right\}\right)<\epsilon.
	\end{equation*}
\end{theorem}

\begin{corollary}\label{cor:GNNapproxPhi_solu}
    Let $\mathcal{D}\subset \Phi_{\text{obj}}^{-1}(\mathbb{R})\subset \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times \mathcal{H}^W_n$ be a finite dataset. For any $\delta>0$, there exists some $F_W\in\mathcal{F}_{\text{GNN}}^W$, such that
	\begin{equation*}
		\| F(G,H) -  \Phi_{\text{solu}}(G,H)\|<\delta,\quad\forall~(G,H)\in\mathcal{D}.
	\end{equation*}
\end{corollary}

\section{Sketch of proof}
\label{sec:sketch}
In this section, we will present a sketch of our proof lines and provide examples to show the intuitions. The full proof lines are presented in the appendix.

\paragraph{\textbf{Separation power}}
Separation power of GNN is closely related to the graph isomorphism problem. To identify whether two given graphs are isomorphic,  one may apply the Weisfeiler-Lehman (WL) test \cite{weisfeiler1968reduction}, a widely studied method that can distinguish a broad class of graphs.
It has been shown in the literature that the separation power of GNN is actually equal to the standard WL test. Thus, we present some concepts of the WL test before analyzing GNN.

\begin{algorithm}[htb!]
	\caption{The WL test for LP-Graphs (denoted by $\mathrm{WL}_{\mathrm{LP}}$)}\label{alg:WL}
	\begin{algorithmic}[1]
		\Require A graph instance $(G,H) \in \mathcal{G}_{m,n}\times \mathcal{H}_{m}^V\times\mathcal{H}_{n}^W$ and iteration limit $L>0$.
		\State Initialize with $C_i^{0,V} = \text{HASH}_{0,V}(h_i^V)$,  $C_j^{0,W} = \text{HASH}_{0,W}(h_j^W)$.
		\For{$l=1,2,\cdots,L$}
		\State $C_i^{l,V} = \text{HASH}_{l, V} \left(C_i^{l-1,V}, \sum_{j=1}^n E_{i,j}\text{HASH}_{l,W}'\left(C_j^{l-1, W}\right)\right)$.
		\State $C_j^{l,W} = \text{HASH}_{l, W} \left(C_j^{l-1,W}, \sum_{i=1}^n E_{i,j} \text{HASH}_{l,V}'\left(C_i^{l-1, V} \right)\right)$.
		\EndFor
		\State \textbf{return} The multisets
		containing all colors $\{\{ C_i^{L,V} \}\}_{i=0}^m, \{\{ C_j^{L,W} \}\}_{j=0}^n$.
	\end{algorithmic}
\end{algorithm}

To apply the WL test on LP-graphs, we describe a modified WL test in Algorithm~\ref{alg:WL}, that is slightly different with the standard WL test.
In Algorithm~\ref{alg:WL}, the hash functions are injective mappings, and multisets are collections of elements that allow multiple appearance of the same element. $\{\{\}\}$ will be used throughout this paper to denote the multisets.
We denote Algorithm~\ref{alg:WL} by $\mathrm{WL}_{\mathrm{LP}}(\cdot)$, and we say that two LP-graphs $(G,H),(\hat{G},\hat{H})$ \emph{can be distinguished by Algorithm~\ref{alg:WL}} if and only if there exist a positive integer $L$ and injective hash functions $\{\mathrm{HASH}_{l,V},\mathrm{HASH}_{l,W}\}_{l=0}^L$ such that
\begin{equation*}
    \mathrm{WL}_{\mathrm{LP}}\big((G,H),L\big) \neq \mathrm{WL}_{\mathrm{LP}}\big((\hat{G},\hat{H}),L\big).
\end{equation*}

Unfortunately, similar with the standard WL test (see e.g. \cite{douglas2011weisfeiler} and references therein), Algorithm~\ref{alg:WL} does not perfectly solve the graph isomorphism problem on all LP-graphs. There exist infinitely many pairs of  non-isomorphic LP-graphs that cannot be distinguished by Algorithm~\ref{alg:WL}. Figure~\ref{fig:not_iso_WL} provide such an example.

\begin{figure}[h]

 \begin{minipage}{0.16\textwidth}
		\begin{tikzpicture}[
			constraintb/.style={circle, draw = blue, scale = 0.8},
			variablex/.style={circle, draw = red, scale = 0.8},
			scale = 0.7,
			]

			\draw (-2,1.5) node[constraintb] (x1) {$v_1$};
			\draw (-2,0.5) node[constraintb] (x2) {$v_2$};
			\draw (-2,-0.5) node[constraintb] (x3) {$v_3$};
			\draw (-2,-1.5) node[constraintb] (x4) {$v_4$};
			\draw (0,1.5) node[variablex] (b1) {$w_1$};
			\draw (0,0.5) node[variablex] (b2) {$w_2$};
			\draw (0,-0.5) node[variablex] (b3) {$w_3$};
			\draw (0,-1.5) node[variablex] (b4) {$w_4$};

			\draw[-] (x1.east) -- (b1.west);
			\draw[-] (x2.east) -- (b2.west);
			\draw[-] (x3.east) -- (b3.west);
			\draw[-] (x4.east) -- (b4.west);
			\draw[-] (x1.east) -- (b2.west);
			\draw[-] (x2.east) -- (b3.west);
			\draw[-] (x3.east) -- (b4.west);
			\draw[-] (x4.east) -- (b1.west);
		\end{tikzpicture}
    \end{minipage}
\begin{minipage}{0.26\textwidth}
\begin{small}
\begin{equation*}
\begin{split}
        \min ~& x_1 + x_2 + x_3 + x_4, \\
        \text{s.t.} ~ & x_1 + x_2 = 1, \\
        & x_2 + x_3 = 1, \\
        & x_3 + x_4 = 1, \\
        & x_4 + x_1 = 1,\\
        & x_j\geq 1,\ 1\leq j\leq 4.
    \end{split}
\end{equation*}
\end{small}
\end{minipage}
\begin{minipage}{0.26\textwidth}
\begin{small}
\begin{equation*}
    \begin{split}
        \min ~ & x_1 + x_2 + x_3 + x_4, \\
        \text{s.t.} ~ & x_1 + x_2 \leq 1, \\
        & x_2 + x_3 \leq 1, \\
        & x_3 + x_4 \leq 1, \\
        & x_4 + x_1 \leq 1,\\
        & x_j\leq 1,\ 1\leq j\leq 4.
    \end{split}
\end{equation*}
\end{small}
\end{minipage}
\begin{minipage}{0.26\textwidth}
\begin{small}
\begin{equation*}
    \begin{split}
        \min ~& x_1 + x_2 + x_3 + x_4, \\
        \text{s.t.} ~ & x_1 + x_2 = 1, \\
        & x_2 + x_3 = 1, \\
        & x_3 + x_4 = 1, \\
        & x_4 + x_1 = 1,\\
        & x_j\leq 1,\ 1\leq j\leq 4.
    \end{split}
\end{equation*}
\end{small}
\end{minipage}

\begin{minipage}{0.16\textwidth}
		\begin{tikzpicture}[
			constraintb/.style={circle, draw = blue, scale = 0.8},
			variablex/.style={circle, draw = red, scale = 0.8},
			scale = 0.7,
			]

			\draw (-2,1.5) node[constraintb] (x1) {$v_1$};
			\draw (-2,0.5) node[constraintb] (x2) {$v_2$};
			\draw (-2,-0.5) node[constraintb] (x3) {$v_3$};
			\draw (-2,-1.5) node[constraintb] (x4) {$v_4$};
			\draw (0,1.5) node[variablex] (b1) {$w_1$};
			\draw (0,0.5) node[variablex] (b2) {$w_2$};
			\draw (0,-0.5) node[variablex] (b3) {$w_3$};
			\draw (0,-1.5) node[variablex] (b4) {$w_4$};

			\draw[-] (x1.east) -- (b1.west);
			\draw[-] (x2.east) -- (b2.west);
			\draw[-] (x3.east) -- (b3.west);
			\draw[-] (x4.east) -- (b4.west);
			\draw[-] (x1.east) -- (b2.west);
			\draw[-] (x2.east) -- (b1.west);
			\draw[-] (x3.east) -- (b4.west);
			\draw[-] (x4.east) -- (b3.west);
		\end{tikzpicture}
		\label{fig:not_iso_WLb}
    \end{minipage}
\begin{minipage}{0.26\textwidth}
\begin{small}
\begin{equation*}
    \begin{split}
        \min ~& x_1 + x_2 + x_3 + x_4, \\
        \text{s.t.} ~ & x_1 + x_2 = 1, \\
        & x_2 + x_1 = 1, \\
        & x_3 + x_4 = 1, \\
        & x_4 + x_3 = 1,\\
        & x_j\geq 1,\ 1\leq j\leq 4.
    \end{split}
\end{equation*}
\end{small}

\end{minipage}
\begin{minipage}{0.26\textwidth}
\begin{small}
\begin{equation*}
    \begin{split}
        \min ~ & x_1 + x_2 + x_3 + x_4, \\
        \text{s.t.} ~ & x_1 + x_2 \leq 1, \\
        & x_2 + x_1 \leq 1, \\
        & x_3 + x_4 \leq 1, \\
        & x_4 + x_3 \leq 1,\\
        & x_j\leq 1,\ 1\leq j\leq 4.
    \end{split}
\end{equation*}
\end{small}
\end{minipage}
\begin{minipage}{0.26\textwidth}
\begin{small}
\begin{equation*}
    \begin{split}
        \min ~& x_1 + x_2 + x_3 + x_4, \\
        \text{s.t.} ~ & x_1 + x_2 = 1, \\
        & x_2 + x_1 = 1, \\
        & x_3 + x_4 = 1, \\
        & x_4 + x_3 = 1,\\
        & x_j\leq 1,\ 1\leq j\leq 4.
    \end{split}
\end{equation*}
\end{small}

\end{minipage}
\caption{LP-graphs that cannot be distinguished by the WL test since neighbor information of $\{v_{i}\}$ and $\{w_j\}$ in the two graphs is equal. Based on the graph pair, we construct three pairs of LPs that are both infeasible, both unbounded, both feasible bounded with the same optimal solution, respectively.}
\label{fig:not_iso_WL}
\end{figure}

However, for two LP-graphs that cannot be distinguished by Algorithm~\ref{alg:WL}, they must share some common characteristics even if they are not isomorphic. For example, let us consider the $6$ LP problems in Figure \ref{fig:not_iso_WL}.
In each of the three columns, the two non-isomorphic LP problems cannot be distinguished by the WL test. It can be checked that the two problems in the same column share some common characteristics. More specifically, both problems in the first column are infeasible; both problems in the second column are feasible but unbounded; both problems in the third column are feasible and bounded with $(1/2,1/2,1/2,1/2)$ as the optimal solution with the smallest $\ell_2$-norm. Actually, this phenomenon does not only happen on the above instances, but also serves as an universal principle for all LP instances. We summarize the results in the following theorem:

\begin{theorem}\label{thm:LP_sameWL2sameproperties-main-txt}
If $(G, H),(\hat{G},\hat{H})\in \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times\mathcal{H}^W_n$ are not distinguishable by Algorithm \ref{alg:WL}, then
\[\Phi_{\text{feas}}(G,H) = \Phi_{\text{feas}}(\hat{G},\hat{H}),\quad \Phi_{\text{obj}}(G,H) = \Phi_{\text{obj}}(\hat{G},\hat{H}).\]
Furthermore, if $(G, H),(\hat{G},\hat{H})\in \Phi_{\text{obj}}^{-1}(\mathbb{R})$, it holds that $\Phi_{\text{solu}}(G,H) = \sigma_W( \Phi_{\text{solu}}(\hat{G},\hat{H}) )$ for some $\sigma_W\in S_n$.
\end{theorem}

In other words, the above theorem guarantees the sufficient power of the WL test for separating LP problems with different characteristics, including feasibility, optimal objective value, and optimal solution with smallest $\ell_2$-norm (up to permutation).
This combined with the following theorem, which states the equivalence of the separation powers of the WL test and GNNs, yields that GNNs also have sufficient separation power for LP-graphs in the above sense.

\begin{theorem}\label{thm:separate_GNN-main-txt}
    For any $(G, H),(\hat{G},\hat{H})\in \mathcal{G}_{m,n}\times\mathcal{H}^V_m\times\mathcal{H}^W_n$, the followings are equivalent:
    \begin{itemize}[leftmargin=8mm]
        \item[(i)] $(G,H)$ and $(\hat{G},\hat{H})$ are not distinguishable by Algorithm \ref{alg:WL}.
        \item[(ii)] $F(G,H) = F(\hat{G},\hat{H}),\ \forall~F\in\mathcal{F}_{\text{GNN}}$.
        \item[(iii)] For any $F_W\in\mathcal{F}_{\text{GNN}}^W$, there exists $\sigma_W\in S_n$ such that $F_W(G,H) = \sigma_W( F_W(\hat{G},\hat{H}))$.
    \end{itemize}
\end{theorem}
Theorem \ref{thm:separate_GNN-main-txt} extends the results in \cites{xu2019powerful,azizian2020expressive,geerts2022expressiveness} to the case with the modified WL test (Algorithm \ref{alg:WL}) and LP-graphs.

\paragraph{\textbf{Representation power}}
Based on the separation power of GNNs, one is able to investigate the representation/approximation power of GNNs. To prove Theorems \ref{thm:GNNapproxPhi_feas}, \ref{thm:GNNapproxPhi_obj}, and \ref{thm:GNNapproxPhi_solu}, we first determine the closure of $\mathcal{F}_{\text{GNN}}$ or $\mathcal{F}_{\text{GNN}}^W$ in the space of invariant/equivariant continuous functions with respect to the $\sup$-norm, which is also named as the \emph{universal approximation}. The result of $\mathcal{F}_{\text{GNN}}$ is stated as follows, where $\mathcal{C}(X,\mathbb{R})$ is the collection/algebra of all real-valued continuous function on $X$. The result of $\mathcal{F}_{\text{GNN}}^W$ can be found in the appendix.

\begin{theorem}\label{thm:uniapprox_scal}
    Let $X\subset \mathcal{G}_{m,n}\times \mathcal{H}^V_m\times \mathcal{H}^W_n$ be a compact set. For any $\Phi\in\mathcal{C}(X,\mathbb{R})$ that satisfies $\Phi(G,H) = \Phi(\hat{G},\hat{H})$ for all $(G,H),(\hat{G},\hat{H})$ that are not distinguishable by Algorithm \ref{alg:WL},
    and for any $\epsilon>0$, there exists $F\in \mathcal{F}_{\text{GNN}}$ such that
    \begin{equation*}
        \sup_{(G,H)\in X} |\Phi(G,H) - F(G,H)| < \epsilon.
    \end{equation*}
\end{theorem}

Theorem~\ref{thm:uniapprox_scal} can be viewed as an LP-graph version of results in \cites{azizian2020expressive, geerts2022expressiveness}. Roughly speaking, graph neural networks can approximate any invariant continuous function with arbitrarily small error. Although our target mappings $\Phi_{\text{feas}},\Phi_{\text{obj}},\Phi_{\text{solu}}$ are not continuous, we prove (in appendix) that they are measurable. Applying Lusin's theorem \cite{evans2018measure}*{Theorem 1.14}, we show that GNN can be arbitrarily close to the target mappings except for a small domain.

\section{Numerical Experiments}
\label{sec:numerical}

We present the numerical results that validate our theoretical results in this section. We generate LP instances with $m = 10$ and $n = 50$ that are possibly infeasible or feasible and bounded.
To check whether GNN can predict feasibility, we generate three data sets with $100,500,2500$ independent LP instances respectively, and call the solver wrapped in \texttt{scipy.optimize.linprog} to get the feasibility, optimal value and an optimal solution for each generated LP.
To generate enough feasible and bounded LPs to check whether GNN can approximate the optimal objective value and optimal solution, we follow the same approach as before to generate LP randomly and discard those infeasible LPs until the number of LPs reach our requirement.
We train GNNs to fit the three LP characteristics by minimizing the distance between GNN-output and those solver-generated answers.
The building and the training of the GNNs are implemented using \texttt{TensorFlow} and the codes are modified from\footnote{Our codes are modified from \url{https://github.com/ds4dm/learn2branch} and we will release codes.} \cite{gasse2019exact}. We set $L=2$ for all GNNs and those learnable functions $f_{\mathrm{in}}^V,f_{\mathrm{in}}^W,f_\text{out}, f^W_{\text{out}}, \{f_l^V,f_l^W,g_l^V,g_l^W\}_{l=0}^{L}$ are all parameterized with MLPs. Details can be found in the appendix. Our results are reported in Figure \ref{fig:experiments-lp}.

\begin{figure}[h]
    \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{lp_fea.pdf}
    \caption{Feasibility}
    \label{fig:feas}
    \end{subfigure}
    \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{lp_obj.pdf}
    \caption{Optimal objective value}
    \label{fig:obj}
    \end{subfigure}
    \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{lp_sol.pdf}
    \caption{Optimal solution}
    \label{fig:solu}
    \end{subfigure}
    \caption{GNN can approximate $\Phi_{\text{feas}}$, $\Phi_{\text{obj}}$, and $\Phi_{\text{solu}}$}
    \label{fig:experiments-lp}
\end{figure}

All the errors reported in Figure \ref{fig:experiments-lp} are training errors since generalization is out of the scope of this paper. In Figure \ref{fig:feas}, the ``rate of errors'' means the proportion of instances with $\mathbb{I}_{F(G,H)>1/2}\neq \Phi_{\text{feas}}(G,H)$. This metric exactly equals to zeros as long as the number of parameters in GNN is large enough, which directly validates Corollary \ref{cor:GNNapproxPhi_feas}. With the three curves in Figure \ref{fig:feas} combined together, we conclude that such principle does not violate as the number of samples increases. This consists with Theorem \ref{thm:GNNapproxPhi_feas}.
Mean squared errors in Figures \ref{fig:obj} and \ref{fig:solu} are respectively defined as $\mathbb{E}_{(G,H)}|F(G,H)-\Phi_{\text{feas}}(G,H)|^2$ and $\mathbb{E}_{(G,H)}\|F(G,H)-\Phi_{\text{feas}}(G,H)\|^2$. Therefore, Figures \ref{fig:obj} and \ref{fig:solu} validates Theorems \ref{thm:GNNapproxPhi_obj} and \ref{thm:GNNapproxPhi_solu} respectively.
Note that all the instances used in Figure \ref{fig:obj} are feasible and bounded. Thus, Figure \ref{fig:obj} actually only validates \eqref{eq:GNNapproxPhi_obj_2} in Theorem \ref{thm:GNNapproxPhi_obj}. However, due to the fact that feasibility of an LP is equal to the boundedness of its dual problem, one may dualize each LP and use the conclusion of Figure \ref{fig:feas} to validate \eqref{eq:GNNapproxPhi_obj_1} in Theorem \ref{thm:GNNapproxPhi_obj}.

\section{Conclusions}
\label{sec:conclusions}

In this work, we show that graph neural networks, as well as the WL test, have sufficient separation power to distinguish linear programming problems with different characteristics. In addition, GNNs can approximate LP feasibility, optimal objective value, and optimal solution with arbitrarily small errors on compact domains or finite datasets. These results guarantee that GNN is a proper class of machine learning models to represent linear programs, and hence contribute to the theoretical foundation in the learning-to-optimize community. Future directions include investing the representation power of graph neural networks for mixed-integer linear programming (MILP), which has been observed with promising experimental results in the literature.

\bibliographystyle{amsxport}