1 Introduction↩︎

Bundle pricing is a widely adopted strategy across industries such as e-commerce, digital subscriptions, and retail. It refers to the practice where a firm provides combinations (i.e., “bundles”) of products or services at discounted prices, supplementing the traditional component pricing (CP) strategy where products are only sold separately. For instance, brands like Tula Skincare actively employ a mixed bundling strategy by offering individual items alongside an impressive range of curated sets and starter kit bundles.¹ As illustrated in Figure 1, these bundles explicitly highlight both the percentage and absolute dollar savings to clearly communicate value and incentivize larger purchases. In the e-commerce sector, KIND Snacks reported a 24% increase in average order value (AOV) after introducing “build-your-own” bundles that empower customers to actively select their own personalized product combinations, while Peet’s Coffee drove a 27% increase in new subscribers within a single year through targeted subscription bundles. To implement this strategy effectively, a firm needs to solve a complex optimization problem: determining which subsets of products to offer and setting their corresponding prices to maximize total profit, under the constraint that heterogeneous customers self-select the options that maximize their own surplus.

The difficulty of the bundle pricing problem stems from two intertwined challenges: the combinatorial explosion of the search space and the complex customer dynamics. In particular, with \(n\) products, the number of possible bundles grows exponentially (\(2^n\)), and customers strategically choose the bundle that maximizes their individual surplus. This self-selection behavior creates highly coupled constraints across all offered bundles, as a price change in one bundle can drastically alter the demand for others. Classical formulations, such as the mixed bundling (MB) model by [1], rigorously capture these customer surplus constraints but become computationally intractable as \(n\) increases (e.g., greater than 15). To address these computational bottlenecks, researchers have proposed various approximation strategies aimed at improving scalability. A notable example is bundle size pricing [2], which simplifies the problem by enforcing identical prices for all bundles of the same size. However, because such methods rely on strong simplifying assumptions—such as completely ignoring product heterogeneity—both their runtime efficiency and final profit performance remain fundamentally limited in large-scale settings.

In recent years, neural networks have demonstrated remarkable success in “Learning to Optimize” (L2O) for combinatorial problems, such as mixed-integer linear programming (MILP) solving via neural branching [3] or neural diving [4]. Motivated by these advances, one might consider leveraging these neural-network-based MILP solvers to accelerate the solution process. However, they share a fundamental dependency on the explicit graph representation of the MILP’s variable-constraint matrix. In the bundle pricing context, the number of potential bundles \(2^n\) leads to a catastrophic memory bottleneck. For any large product size \(n\), the resulting bipartite graph becomes too massive to be stored or processed by a GNN, rendering these neural heuristics computationally infeasible at the very stage of model initialization.

To overcome these challenges, we investigate how to design a scalable GNN framework that learns to identify high-quality bundle structures directly from data, without constructing the explicit exponential graph. Specifically, we aim to address the following research questions:

1. Representation: How can we leverage graph neural networks to extract features from the bundle pricing problem without explicitly constructing the exponentially large variable-constraint graph?

2. Search Space Reduction: Given the neural network predictions, how can we design inference policies that effectively prune the intractable \(2^n\) search space into a manageable subset?

3. Scalability: How does the performance of the neural-network-based approach scale with problem size?

To address the first research question, we propose a GNN-based framework that learns to identify the latent structure of the bundle pricing problem and hence provides high-quality solutions. Leveraging the structural properties of this problem, we construct a bipartite graph representing these customer-product interactions, which scales linearly with the number of segments and products. We then train a GNN on this compact representation to extract latent structural features characterizing customer-product affinities. These features are decoded to predict inclusion-probability matrix, representing the probability that each customer segment will purchase each product in the optimal solution. Theoretically, we show that, under mild distinguishability conditions, the edge-output GNN class can represent the optimal product-assignment mapping and make the binary cross-entropy loss arbitrarily small.

To address the second research question, we design pruning-based inference policies that translate the GNN output into a restricted bundle family. Fixed cutoff pruning (FCP) forms one candidate bundle per segment by retaining products whose predicted probabilities exceed a cutoff. Progressive cutoff pruning (PCP) is more conservative: for each segment, it ranks retained products by predicted probability and constructs a sequence of nested prefix bundles. In both cases, the seller then solves the [1]’s mixed-bundling formulation over the retained candidate family for the joint optimization of bundle offerings and pricing. To further refine the solution quality, we incorporate a GNN-guided local search method that strategically adds or drops products from the assigned bundles.

To address the third research question, we evaluate the framework across problem sizes. Specifically, the GNN’s weight-sharing capability enables us to train the model on tractable small-scale instances (where optimal labels are efficiently computable) and directly apply it to intractable large-scale ones. Noticing that the GNN’s scalability is limited, we further introduce an iterative self-improvement procedure: a GNN trained on small optimal labels is first used to generate high-quality near-optimal labels on larger instances, and these near-optimal labels are then used to train a larger-scale model. This procedure expands the training regime without requiring exact large-scale mixed-bundling labels. Numerical experiments are provided to illustrate its effectiveness.

Our contributions are threefold. First, we introduce a compact learning representation for mixed bundling. Instead of predicting over the exponential bundle space, we learn the segment-product projection of an optimal mixed-bundling solution on a customer-product graph. Second, we develop a GNN-guided pruning-then-optimization framework, together with a GNN-guided local-search procedure to further refine the candidate family, and an iterative self-improvement framework to further enhance the GNN performance on large-scale instance. Third, we provide theoretical and computational evidence for the proposed framework. Theoretically, we establish an expressiveness result for edge-output GNNs under mild conditions. Computationally, we show that the proposed policies achieve near-optimal profits on small instances and outperform scalable bundle-pricing benchmarks on larger instances.

2 Literature Review↩︎

This work is broadly related to three streams of research: bundle pricing, neural-network-based optimization, and bundle recommendation.

Bundle pricing. The study of bundle pricing has a long history in economics and operations research communities. Early work in the two-product setting established that bundling can be more profitable than separate sales, including pure bundling (PB), component pricing (CP), and mixed bundling (MB) [5]–[7]. Later research extends these insights to large-scale settings. A stream of work analyzes simplified mechanisms such as PB and bundle-size pricing (BSP) [2], [8]–[10]. In particular, PB can approximate the revenue of MB when the number of products grows large under zero costs and independent additive valuations [8], while later work studies large-scale bundling with positive production costs [9]. Empirical and theoretical studies show that BSP can outperform CP and PB in certain conditions and can achieve asymptotic optimality when product costs are homogeneous, though it may deteriorate under cost heterogeneity [2], [10]–[12]. More recently, novel mechanisms have been proposed to address different limitations of classical bundling schemes. Pure bundling with disposal for cost (PBDC) mitigates the adverse effect of high production costs [13]; loot box mechanisms study randomized adaptive bundling schemes [14]; and component pricing with bundle-size discounts (CPBSD) subsumes CP, PB, BSP, PBDC, and two-part tariffs as special cases and inherits constant-factor guarantees under additive independent valuations [15]. Complementary to mechanism design, a computational line of work models bundle pricing directly as a mixed-integer program. The seminal MB formulation explicitly models customer self-selection and price subadditivity constraints [1], though its scalability is limited by the exponential growth of candidate bundles. Our work partially follows this computational tradition: we retain the mixed bundling framework while using graph-based learning to construct a restricted candidate bundle family before invoking the MILP solver, thereby combining the structural rigor of MB with practical scalability.

Neural-network-based optimization. Recently, the use of machine learning to accelerate the solution of optimization problems has attracted increasing attention. In MILP, a seminal work by [3] encodes an instance as a constraint–variable bipartite graph and trains a GCN to imitate strong branching decisions. Subsequent studies learn primal heuristics, branching, diving, cut selection, separator configuration, and local-search decisions within B&B or branch-and-cut [4], [16]–[21]. A related predict-then-optimize line uses neural predictions to restrict or guide the subsequent search: [22] predict marginal variable probabilities and search around the predicted solution, [23] alternate prediction and correction to fix reliable variables, and [24] use probabilistic multi-variable branching to partition MIP feasible regions. Recent generative and differentiable approaches further expand this direction: [25] formulate ILP solving as an unsupervised differentiable learning problem, and [26] use flow matching to model the joint continuous–integer solution distribution of MILPs. In the context of linear programming, [27] propose GCN-based initial basis selection to warm-start simplex and column generation, while [28] train GCNs to imitate expert pivot policies, shortening pivot paths. On the theory side, some researchers study the capability of GNNs (see [29] for a comprehensive review). [30] show that GNNs can represent key properties of linear programs, including feasibility, boundedness, and optimal solutions. For MILPs, [31] characterize the expressive limitations of standard message-passing GNNs and identify conditions under which these limitations can be overcome. Our theoretical result for GNNs is different in an important way: rather than representing and predicting the solution of the full mixed-bundling MILP over the exponential bundle-level decision space, we learn a polynomial-size product-assignment projection of the optimal solution, with one edge-level prediction for each segment-product pair.

Machine learning has also been directly integrated into the core structure of operations research problems. Attention-based and reinforcement-learning models have been developed for routing problems such as TSP and VRP [32], [33]. [34] integrate first-order methods with a self-supervised neural network framework to efficiently solve large-scale, constrained assortment optimization problems across various complex choice models. [35] leverage the generalizability of GCNs to solve large-scale NP-hard constrained assortment optimization problems under the mixed multinomial logit choice model. In a similar vein, we apply GNNs to the bundle pricing problem. Instead of constructing the exponential bundle-level formulation, our GNN learns segment–product interaction patterns on a compact customer-product graph and predicts product-assignment probabilities, which are then converted into a restricted candidate bundle family. This preserves [1]’s pricing and self-selection framework while using GNN predictions to accelerate large-scale instances.

Bundle recommendation. Our work is also related to bundle recommendation; we refer readers to [36], [37], [38], [39] and references therein for comprehensive reviews and representative methods. Recommendation research typically predicts user preference scores or personalized rankings for bundles from historical user–item, user–bundle, and bundle–item interactions, and is usually evaluated by ranking metrics such as Recall and NDCG. This differs from bundle pricing in two essential ways. First, recommendation methods usually treat price as exogenous or absent and therefore do not optimize a market-wide price menu. Second, they do not impose self-selection constraints across customer segments. Introducing pricing makes the problem substantially harder: in a full mixed-bundling formulation, the seller may need a price variable for each of the exponential number of possible bundles, and each price can change all customers’ surplus-maximizing choices.

3 Problem Definition↩︎

We study a bundle pricing problem where a seller offers \(n\) products to \(m\) customer segments. Let \([n]=\{1,\dots,n\}\) be the product set, \([m]=\{1, \dots, m\}\) the customer set, and \(\mathfrak F:=2^{[n]}\) the full family of product bundles, including the empty bundle. We use \(\mathcal{K}:=\{0,1,\dots,L\}\) with \(L=2^n-1\), as the corresponding bundle index set. We fix a bijection \(\mathcal{A}:\mathcal{K}\to\mathfrak F,\) where \(\mathcal{A}(b)\subseteq[n]\) denotes the product set represented by bundle index \(b\). A candidate bundle family is denoted by \(\mathfrak B\subseteq\mathfrak F,\) and its associated index set is \(\mathcal{I}_{\mathfrak B} := \{b\in\mathcal{K}: \mathcal{A}(b)\in\mathfrak B\}.\) The full mixed-bundling problem corresponds to \(\mathfrak B=\mathfrak F\), or equivalently \(\mathcal{I}_{\mathfrak B}=\mathcal{K}\).

Each customer segment \(k \in [m] = \{1, \dots, m\}\) represents a cohort of homogeneous consumers, characterized by a proportion \(\alpha_k\) and a distinct valuation vector \(\mathbf{v}_k=\{R_{kb}\mid b\in\mathcal{K}\},\) where \(R_{kb}\) denotes segment \(k\)’s maximum willingness-to-pay for the product set \(\mathcal{A}(b)\). The index \(0\) corresponds to the empty bundle, representing the outside option, and we normalize \(R_{k0}=0\) for \(k\in[m]\). For a candidate bundle family \(\mathfrak B\), the seller chooses the prices \(\{p_b\mid b\in\mathcal{I}_{\mathfrak B}\}\) to maximize the total profit, which can be defined as \(\sum_{k=1}^{m}\sum_{b\in\mathcal{I}_{\mathfrak B}} (p_b-c_{kb})\cdot \alpha_k \theta_{kb}.\) Here, \(c_{kb}\) is the serving cost of selling bundle index \(b\) to segment \(k\), and \(\theta_{kb}\in\{0,1\}\) is a binary decision variable denoting whether segment \(k\) selects bundle index \(b\).

The demand for bundles follows standard assumptions: each segment \(k\) has a valuation \(R_{kb}\) for bundle index \(b\). After observing all displayed bundles in \(\mathfrak B\), each buyer selects exactly one bundle index \(b\in\mathcal{I}_{\mathfrak B}\) that maximizes its individual surplus, defined by \(\Delta_{kb}=R_{kb}-p_b,\) provided that the surplus is non-negative. We follow the assumption of free disposal of unwanted products within a bundle [1]. Under this setting, customers can discard unwanted items at no cost. Consequently, if a superset bundle \(b_2\in\mathcal{I}_{\mathfrak B}\) were priced lower than its subset bundle \(b_1\in\mathcal{I}_{\mathfrak B}\), i.e., \(\mathcal{A}(b_1)\subseteq \mathcal{A}(b_2)\) and \(p_{b_2}<p_{b_1},\) then rational customers desiring the product set \(\mathcal{A}(b_1)\) could simply purchase bundle \(b_2\), discard the items in \(\mathcal{A}(b_2)\setminus \mathcal{A}(b_1)\), and effectively obtain \(\mathcal{A}(b_1)\) at a lower price. Thus, to preclude such arbitrage opportunities, bundle prices should be non-decreasing with respect to set inclusion.

Under these assumptions, for a candidate bundle family \(\mathfrak B\), the seller’s profit optimization problem can be formulated as follows: \[\begin{align} \max \quad & \sum_{k=1}^m \sum_{b \in \mathcal{I}_{\mathfrak B}} (p_b - c_{kb}) \cdot \alpha_k \, \theta_{kb} \\ \text{s.t.} \quad & \theta_{kb} \leq \mathbb{1}\!\left( b \in \arg\max_{b' \in \mathcal{I}_{\mathfrak B}} \{R_{kb'} - p_{b'}\} \;\wedge\; R_{kb} - p_b \ge 0 \right), \quad \forall k,\;b \in \mathcal{I}_{\mathfrak B} \\[6pt] & \sum_{b \in \mathcal{I}_{\mathfrak B}} \theta_{kb} = 1, \quad \forall k \\[6pt] & p_b \ge 0, \quad \forall b \in \mathcal{I}_{\mathfrak B}. \end{align}\] Here, \(\mathbb{1}(\cdot)\) is the indicator function that enforces buyers’ surplus maximization choice.

In the following discussion, we set \(\mathcal{A}(0)=\emptyset\). For any bundle index \(b\in\mathcal{K}\), we define the valuation of segment \(k\) as \(R_{kb} = f\!\left(\sum_{j\in \mathcal{A}(b)}u_{kj}\right),\) where \(u_{kj}\) is the utility of product \(j\) for segment \(k\), and \(f(\cdot)\) is an increasing concave function with \(f(0)=0\). For any bundle index \(b\in\mathcal{K}\setminus \{0\}\), we assume that the serving cost is \(c_{kb}= \sum_{j\in \mathcal{A}(b)} c_j^{\mathit u}+c_k^s,\) where \(c_k^s\) is the cost associated with serving customer segment \(k\) (e.g., the shipping cost) and \(c_j^{\mathit u}\) is the unit cost of product \(j\). For the outside option, \(c_{k0}=0\) for all \(k \in [m]\). The concavity of \(f(\cdot)\) captures the economic consensus of diminishing marginal utility, and the cost structure captures the principle of economies of scale: since the fixed cost is spread across all items in the bundle, the average cost per item decreases as more products are added.

To solve this problem, [1] formulate a mixed-integer linear program (MILP) with an exponential number of variables. For the sake of completeness, we describe this formulation in the following.

3.1 Mixed Bundling Formulation [1]↩︎

Consider a candidate bundle family \(\mathfrak B\subseteq\mathfrak F\) with associated index set \(\mathcal{I}_{\mathfrak B} = \{b\in\mathcal{K}: \mathcal{A}(b)\in\mathfrak B\}.\) Define \(R_{\max}:=\max_{k\in[m],\, b\in\mathcal{K}} R_{kb}.\)

We define the following variables: \(\theta_{kb}\) is a binary indicator for whether segment \(k\) chooses bundle index \(b\); \(p_b\) is the price of bundle index \(b\); \(P_{kb}\) is the effective price paid by segment \(k\); \(s_k\) is the surplus of segment \(k\); \(Z_{kb}\) is the profit from assigning bundle index \(b\) to segment \(k\); and \(\mathcal{A}(b)\) is the set of components which define bundle index \(b\). We denote \(\boldsymbol{\Theta}:=[\theta_{kb}]\), \(\mathbf{p}:=(p_b)\), \(\mathbf{P}^{\mathrm{pay}}:=[P_{kb}]\), \(\mathbf{s}:=(s_k)\), and \(\mathbf{Z}:=[Z_{kb}]\). The mixed bundle pricing problem over \(\mathfrak B\) can be written as follows: \[\begin{align} \small \Xi(\mathfrak B):=\max \quad & \sum_{k=1}^m \sum_{b \in \mathcal{I}_{\mathfrak B}} \alpha_k \cdot Z_{kb} \\ \text{s.t.} \quad & \sum_{b \in \mathcal{I}_{\mathfrak B}} \theta_{kb} = 1, && \forall k \\ & s_k \ge R_{kb} - p_b, && \forall k,b \in \mathcal{I}_{\mathfrak B} \\ & p_b - R_{\max}(1 - \theta_{kb}) \leq P_{kb}, && \forall k, b \in \mathcal{I}_{\mathfrak B} \\ & P_{kb} \leq p_b, && \forall k, b \in \mathcal{I}_{\mathfrak B} \\ & s_k = \sum_{b \in \mathcal{I}_{\mathfrak B}} (R_{kb}\theta_{kb} - P_{kb}), && \forall k \\ & R_{kb}\theta_{kb} - P_{kb} \ge 0, && \forall k, b \in \mathcal{I}_{\mathfrak B} \\ & s_k \geq \sum_{b \in \mathcal{I}_{\mathfrak B}} (R_{kb}\theta_{\gamma b} - P_{\gamma b}), && \forall k,\gamma \\ & Z_{kb} = P_{kb} - c_{kb}\theta_{kb}, && \forall k, b \in \mathcal{I}_{\mathfrak B} \\ & p_b \le \sum_{c\in\mathcal{C}} p_c, && \forall b\in\mathcal{I}_{\mathfrak B},\; \forall \mathcal{C}\subseteq \mathcal{I}_{\mathfrak B} \text{ with } \mathcal{A}(b)\subseteq \bigcup_{c\in\mathcal{C}} \mathcal{A}(c) \label{eq:hm95cover95subadditivity} \\ & p_b, P_{kb},s_k \geq 0, R_{k0}\theta_{k0} = P_{k0}, \theta_{kb} \in \{0,1\} \end{align}\tag{1}\]

The full mixed-bundling formulation is obtained as the special case \(\Xi(\mathfrak F), \mathfrak F=2^{[n]},\) where all product subsets are offered as candidate bundles. A restricted candidate family should be interpreted as a new pricing problem over the retained menu, not simply as the set of bundles purchased in a full-menu optimum. Appendix 9.3.2 provides an example illustrating this distinction.

In \(\Xi(\mathfrak B)\), price subadditivity is imposed in the cover form 1 : a retained collection \(\mathcal{C}\subseteq\mathcal{I}_{\mathfrak B}\) covers a target bundle \(b\in\mathcal{I}_{\mathfrak B}\) if \(\mathcal{A}(b)\subseteq\bigcup_{c\in\mathcal{C}}\mathcal{A}(c)\), and the model requires \(p_b\le\sum_{c\in\mathcal{C}}p_c\). This form is used for restricted candidate families because they are generally not closed under partitions; see Appendix 9.3 for its connection to full-universe partition subadditivity.

We note that the main computational challenge of the above formulation lies in the exponential number of variables corresponding to all possible bundles. In the following, we demonstrate how the GNN framework can be leveraged to prune the search space and help solve large-scale bundle pricing problems efficiently.

4 GNN-Based Strategies↩︎

In this section, we design GNN-based strategies to solve the optimal bundle pricing problem. Specifically, we employ a GNN architecture characterized by an edge-centric update mechanism. Unlike standard graph convolutional networks (GCNs) that primarily focus on learning node embeddings over static edge attributes, our framework executes a dual update process: it iteratively refines both node and edge embeddings. Our model predicts a segment–product inclusion probability matrix \(\mathbf{P} \in \mathbb{R}^{m\times n}\), where \(\mathbf{P}_{kj}\) is the predicted probability that customer segment \(k\) selects product \(j\), derived directly from the final learned edge embeddings. These probabilities serve as a compact representation of heterogeneous preferences and provide the basis for pruning the exponentially large bundle space. We describe our strategies in detail in the following sections.

4.1 Motivation for GNN↩︎

We motivate our choice of architecture by analyzing the structural limitations of standard Multilayer Perceptrons (MLPs) in the context of bundle pricing, and explaining how GNNs effectively address these challenges:

4.1.0.1 Limitations of MLPs.

Standard feed-forward networks face two fundamental bottlenecks when applied to combinatorial problems of this nature:

1. Inability to Generalize Across Problem Sizes: MLPs are constrained by fixed input and output dimensions. Consequently, instances with different numbers of products or customer segments typically require training separate models. While zero-padding can accommodate smaller instances, this approach is inefficient and fails to generalize to instances larger than the training configuration. In contrast, GNNs operate on graph topologies of arbitrary size using shared weights, allowing for seamless scalability.

2. Lack of Permutation Equivariance: MLPs treat inputs as flattened feature vectors and lack the inductive bias to capture the relational structure among agents. In our context, reordering product or customer indices should not alter the optimal product assignment (i.e., the target mapping is permutation equivariant). However, because MLPs assign weights based on specific input positions, permuting the input order arbitrarily changes the output. GNNs naturally respect this symmetry by applying the same local operations across all nodes and edges, ensuring consistent predictions regardless of indexing.

4.1.0.2 Advantages of GNNs.

In contrast, GNNs operate directly on graph-based representations, where nodes correspond to entities (products and customer segments) and edges encode their interactions (utilities).

1. Scalability through Weight Sharing: A key feature of GNNs is the weight-sharing mechanism, where the same learnable weight matrices are applied across all nodes and edges. This decouples the number of parameters from the graph size, allowing a model trained on small-scale instances to generalize effectively to large-scale instances without retraining.

2. Context-Aware Representations: Through message passing, GNNs explicitly model the bipartite interactions. Each node iteratively aggregates information from its neighbors, building a representation that reflects both its local attributes (e.g., costs) and its structural context (e.g., demand from connected segments). This enables the model to capture the latent patterns of optimal bundles in a manner that is both scalable and permutation invariant.

4.2 Graph Representation↩︎

Figure 2 illustrates the graph representation of the optimal bundle pricing problem, which serves as the input to our Graph Neural Network (GNN). The graph consists of two types of nodes: product nodes (\(\square\)) and customer nodes (\(\bigcirc\)), as shown in Figure 2.

Let \(\mathcal{V}\) and \(\mathcal{E}\) denote the sets of nodes and edges, respectively. We define the raw features for the problem instance as follows:

• Product Features: Let \(\mathbf{y}_j^P \in \mathbb{R}^4\) denote the initial feature vector for product \(j\) (\(1 \le j \le n\)). It is constructed as \(\mathbf{y}_j^P=\bigl[c_j^{\mathit{u}},\;\frac{1}{m}\sum_{k=1}^m u_{kj},\;0,\;0\bigr]\).

• Customer Features: Let \(\mathbf{y}_k^S \in \mathbb{R}^4\) denote the initial feature vector for customer segment \(k\) (\(1 \le k \le m\)). It is constructed as \(\mathbf{y}_k^S=\bigl[0,\;0,\;\alpha_k,\;c_k^s\bigr]\).

• Edge Features: Each edge connecting product \(j\) and customer \(k\) is characterized by the scalar utility value \(u_{kj}\).

4.3 GNN Structure↩︎

We implement our model following the GNN framework proposed by [40], utilizing the Generalized Graph Convolution [41] for effective message aggregation. Specifically, we construct an edge-output graph neural network consisting of \(\omega\) graph network blocks followed by a shared scalar edge head. Each graph network block follows the same message-passing template: node embeddings are first updated by aggregating information from neighboring nodes and incident edges, and edge embeddings are then updated using the newly updated endpoint node embeddings. After \(\omega\) graph network blocks, the scalar edge head maps each final product–segment edge embedding to a logit, which is subsequently transformed into a segment-product inclusion probability by a sigmoid function.

4.3.0.1 General Layer Notation.

Let \(G^{(l)} = (\mathbf{V}^{(l)}, \mathbf{Z}^{(l)})\) denote the graph embeddings at layer \(l\). The node feature matrix \(\mathbf{V}^{(l)} \in \mathbb{R}^{(n+m) \times d_v^{(l)}}\) is defined as the vertical concatenation of the node feature vectors. The matrix is constructed as: \[\mathbf{V}^{(l)} = \begin{bmatrix} \mathbf{v}_1^{(l)} \\ \vdots \\ \mathbf{v}_{n+m}^{(l)} \end{bmatrix}.\] Here, \(\mathbf{v}_w^{(l)} \in \mathbb{R}^{1 \times d_v^{(l)}}\) represents the embedding vector of the \(w\)-th node, and \(d_v^{(l)}\) denotes the embedding dimensionality at layer \(l\).

Furthermore, \(\mathbf{V}^{(l)}\) can be explicitly partitioned into product and customer components. We denote \(\mathbf{Y}^{P,(l)} \in \mathbb{R}^{n \times d_v^{(l)}}\) as the embedding matrix for the \(n\) product nodes and \(\mathbf{Y}^{S,(l)} \in \mathbb{R}^{m \times d_v^{(l)}}\) as the embedding matrix for the \(m\) customer nodes. Thus, \(\mathbf{V}^{(l)} = \begin{bmatrix} \mathbf{Y}^{P,(l)} \\ \mathbf{Y}^{S,(l)} \end{bmatrix}.\)

Similarly, the edge embedding tensor \(\mathbf{Z}^{(l)}\in\mathbb{R}^{n\times m\times d_e^{(l)}}\) collects the embeddings for all product–segment pairs. The relationship between the tensor \(\mathbf{Z}^{(l)}\) and an individual edge embedding vector \(\mathbf{e}_{jk}^{(l)}\in\mathbb{R}^{d_e^{(l)}}\) is \(\mathbf{Z}^{(l)}[j,k,:]=\mathbf{e}_{jk}^{(l)}\), for \(j\in[n],\;k\in[m].\) Here \(d_e^{(l)}\) denotes the edge embedding dimension at layer \(l\). The input edge feature is one-dimensional \((d_e^{(0)}=1)\), while all graph network blocks output hidden edge embeddings of dimension \(d_e^{(l)}=d_{\mathrm{mid}}\). Therefore, after the \(\omega\)-th graph network block, the final hidden edge embedding tensor is \(\mathbf{Z}^{(\omega)} \in \mathbb{R}^{n\times m\times d_{\mathrm{mid}}}.\) A shared scalar edge head, introduced below, maps each vector \(\mathbf{e}_{jk}^{(\omega)}\in\mathbb{R}^{d_{\mathrm{mid}}}\) to a scalar logit.

4.3.0.2 Initialization (\(l=0\)).

The input features for the GNN are initialized directly from the raw feature vectors defined in Section 4.2. For node features, we map the product and customer features to the global node index \(w\): \[\mathbf{v}_w^{(0)} = \begin{cases} \mathbf{y}_w^P & \text{if } 1 \le w \le n \quad (\text{Product Node}) \\ \mathbf{y}_{w-n}^S & \text{if } n < w \le n+m \quad (\text{Customer Node}) \end{cases}\] For edge features, we initialize the tensor with the utility values \(\mathbf{e}_{jk}^{(0)} = [u_{kj}]\).

Each GNN layer executes two update functions sequentially: a node update followed by an edge update. We provide a schematic illustration of the GNN architecture in Figure 3. Below, we describe the computations performed in each layer.

Node Update (\(\phi^v\)). In the first stage of layer \(l\), we update the embeddings of all nodes to capture the structural information from their neighbors. We employ the Generalized Graph Convolution of [41] as the update function. For a target node \(w\) and its neighbor set \(\mathcal{N}(w)\), the update consists of five operations:

4.3.0.3 1. Message Construction.

For a directed message from neighbor \(z\) to target node \(w\), we first project the source-node embedding, the target-node embedding, and the edge embedding into a common hidden dimension \(d_{\mathrm{mid}}\): \(\bar{\mathbf{v}}_{z\to w}^{(l-1)} = \mathbf{v}_z^{(l-1)}\mathbf{W}_{\mathrm{s}}^{(l)}, \bar{\mathbf{v}}_{w}^{(l-1)} = \mathbf{v}_w^{(l-1)}\mathbf{W}_{\mathrm{t}}^{(l)}, \bar{\mathbf{e}}_{zw}^{(l-1)} = \mathbf{e}_{zw}^{(l-1)}\mathbf{W}_{\mathrm{e}}^{(l)}.\) Here \(\bar{\mathbf{v}}_{z\to w}^{(l-1)},\; \bar{\mathbf{v}}_{w}^{(l-1)},\; \bar{\mathbf{e}}_{zw}^{(l-1)} \in\mathbb{R}^{1\times d_{\mathrm{mid}}}.\) The matrices \(\mathbf{W}_{\mathrm{s}}^{(l)}, \mathbf{W}_{\mathrm{t}}^{(l)}, \mathbf{W}_{\mathrm{e}}^{(l)}\) are learnable parameters shared across all nodes and edges in layer \(l\). If an input dimension already matches \(d_{\mathrm{mid}}\), the corresponding projection can be interpreted as the identity map.

The message from \(z\) to \(w\) is constructed as \[\mathbf{m}_{z\to w}^{(l)} = \operatorname{ReLU} \left( \bar{\mathbf{v}}_{z\to w}^{(l-1)} + \bar{\mathbf{e}}_{zw}^{(l-1)} \right) + \epsilon_{\mathrm{mp}}\mathbf{1}, \qquad \forall z\in\mathcal{N}(w).\] Here, \(\operatorname{ReLU}(x)=\max\{x,0\}\) is applied element-wise when \(x\) is a vector. The constant \(\epsilon_{\mathrm{mp}}>0\) is a fixed numerical-stability constant and is not learned.

4.3.0.4 2. Softmax Aggregation.

We aggregate incoming messages using a Softmax mechanism, which assigns data-dependent weights to different neighbors: \[\hat{\mathbf{m}}_w^{(l)} = \sum_{z\in\mathcal{N}(w)} \left( \frac{ \exp\!\left(\beta_{\mathrm{mp}}\mathbf{m}_{z\to w}^{(l)}\right) }{ \sum_{\gamma\in\mathcal{N}(w)} \exp\!\left(\beta_{\mathrm{mp}}\mathbf{m}_{\gamma\to w}^{(l)}\right) } \right) \odot \mathbf{m}_{z\to w}^{(l)}.\] Here, \(\exp(\cdot)\), division, and multiplication are applied coordinate-wise, and \(\odot\) denotes element-wise multiplication. The scalar \(\beta_{\mathrm{mp}}>0\) is the inverse-temperature hyperparameter of the Softmax aggregator.

4.3.0.5 3. Residual Update.

The aggregated message is combined with the target node’s own projected embedding through a residual connection and passed through an MLP: \(\mathbf{h}_w^{(l)} = \operatorname{MLP}_{v}^{(l)} \left( \bar{\mathbf{v}}_{w}^{(l-1)} + \hat{\mathbf{m}}_w^{(l)} \right), \quad \mathbf{h}_w^{(l)}\in\mathbb{R}^{1\times d_{\mathrm{mid}}}.\) An MLP is a feed-forward neural network obtained by composing affine transformations and nonlinear activation functions. The same \(\operatorname{MLP}_{v}^{(l)}\) is shared across all nodes in layer \(l\).

4.3.0.6 4. Activation.

We then apply a ReLU activation to the node-update output \(\mathbf{v}_w^{(l)} = \operatorname{ReLU} \left( \mathbf{h}_w^{(l)} \right).\)

4.3.0.7 5. Node Dropout.

To reduce overfitting during training, we apply dropout to the updated node embeddings before they are used in the edge update. Let \(p_D\in[0,1)\) denote the dropout probability. For each node \(w\) and feature coordinate \(r\), let \(v_{wr}^{(l)}\) be the \(r\)-th entry of the updated node embedding \(\mathbf{v}_w^{(l)}\). We define the post-dropout feature coordinate \(\widetilde{v}_{wr}^{(l)}\) by \[\widetilde{v}_{wr}^{(l)} = \begin{cases} \dfrac{\xi_{wr}^{(l)}}{1-p_D}v_{wr}^{(l)}, & \text{during training},\\[3mm] v_{wr}^{(l)}, & \text{during validation and testing}, \end{cases} \qquad \xi_{wr}^{(l)}\sim \operatorname{Bernoulli}(1-p_D).\]

Edge Update (\(\phi^e\)). In the second stage of layer \(l\), we explicitly update the embedding of each product–segment edge. Unlike standard GNNs where edge features may remain static, our framework refines edge features based on the newly updated endpoint node embeddings and the edge’s own history.

For the edge connecting product \(j\) and customer segment \(k\), whose customer node has global index \(n+k\), we construct the concatenated input vector \(\mathbf{z}_{jk}^{(l)} = \widetilde{\mathbf{v}}_{j}^{(l)} \Vert \widetilde{\mathbf{v}}_{n+k}^{(l)} \Vert \mathbf{e}_{jk}^{(l-1)},\) where \(\Vert\) denotes concatenation, and \(\widetilde{\mathbf{v}}_{j}^{(l)}\) and \(\widetilde{\mathbf{v}}_{n+k}^{(l)}\) are the post-dropout endpoint node embeddings. This vector is processed by a two-layer edge MLP followed by layer normalization: \[\mathbf{e}_{jk}^{(l)} = \operatorname{LayerNorm}_{e}^{(l)} \left( \operatorname{ReLU} \left( \mathbf{z}_{jk}^{(l)}\mathbf{W}_{e,1}^{(l)} + \mathbf{b}_{e,1}^{(l)} \right) \mathbf{W}_{e,2}^{(l)} + \mathbf{b}_{e,2}^{(l)} \right), \qquad \mathbf{e}_{jk}^{(l)}\in\mathbb{R}^{1\times d_{\mathrm{mid}}}.\] Here, \(\mathbf{W}_{e,1}^{(l)} \in \mathbb{R}^{(2d_{\mathrm{mid}}+d_e^{(l-1)})\times d_{\mathrm{mid}}}, \mathbf{W}_{e,2}^{(l)} \in \mathbb{R}^{d_{\mathrm{mid}}\times d_{\mathrm{mid}}}, \mathbf{b}_{e,1}^{(l)},\mathbf{b}_{e,2}^{(l)} \in \mathbb{R}^{1\times d_{\mathrm{mid}}}\) are learnable parameters of the edge update network. The parameters of \(\operatorname{LayerNorm}_{e}^{(l)}\), defined below, are also learnable. The same edge update network is shared across all product–segment edges in layer \(l\).

Layer normalization, denoted by \(\operatorname{LayerNorm}_{e}^{(l)}\), normalizes the coordinates of a feature vector and then applies learnable scale and shift parameters. Specifically, for a vector \(\mathbf{x}\in\mathbb{R}^{1\times d}\), \[\operatorname{LayerNorm}_{e}^{(l)}(\mathbf{x}) = \boldsymbol{\gamma}_{e}^{(l)} \odot \frac{\mathbf{x}-\mu(\mathbf{x})\mathbf{1}}{\sqrt{\sigma^2(\mathbf{x})+\varepsilon_{\mathrm{LN}}}} + \boldsymbol{\delta}_{e}^{(l)}.\] Here, \(\mu(\mathbf{x})\) and \(\sigma^2(\mathbf{x})\) are the mean and variance of the coordinates of \(\mathbf{x}\), \(\varepsilon_{\mathrm{LN}}>0\) is a fixed numerical-stability constant, and \(\boldsymbol{\gamma}_{e}^{(l)},\boldsymbol{\delta}_{e}^{(l)} \in\mathbb{R}^{1\times d}\) are learnable parameters.

Output. After \(\omega\) graph network blocks, each product–segment edge has a final hidden embedding \(\mathbf{e}_{jk}^{(\omega)}\in\mathbb{R}^{1\times d_{\mathrm{mid}}}, j\in[n],\;k\in[m].\) Equivalently, \(\mathbf{Z}^{(\omega)} \in \mathbb{R}^{n\times m\times d_{\mathrm{mid}}}\) collects all final hidden edge embeddings.

We apply dropout to the terminal edge embeddings during training. For each edge coordinate \(r\in[d_{\mathrm{mid}}]\), let \(e_{jk,r}^{(\omega)}\) be the \(r\)-th entry of \(\mathbf{e}_{jk}^{(\omega)}\). Define \[\widetilde{e}_{jk,r}^{(\omega)} = \begin{cases} \dfrac{\zeta_{jk,r}^{(\omega)}}{1-p_D}e_{jk,r}^{(\omega)}, & \text{during training},\\[3mm] e_{jk,r}^{(\omega)}, & \text{during validation and testing}, \end{cases} \qquad \zeta_{jk,r}^{(\omega)} \sim \operatorname{Bernoulli}(1-p_D).\] Let \(\widetilde{\mathbf{e}}_{jk}^{(\omega)}\) denote the resulting post-dropout terminal edge embedding. Dropout introduces no learnable parameters; \(p_D\) is a hyperparameter.

A shared scalar edge head \(\rho_e:\mathbb{R}^{1\times d_{\mathrm{mid}}}\to\mathbb{R}\) then maps each terminal edge embedding to a scalar logit: \[\Lambda_{jk} = \rho_e\!\left(\widetilde{\mathbf{e}}_{jk}^{(\omega)}\right) = \widetilde{\mathbf{e}}_{jk}^{(\omega)}\mathbf{w}_o+b_o, \qquad j\in[n],\;k\in[m],\] where \(\mathbf{w}_o\in\mathbb{R}^{d_{\mathrm{mid}}\times 1}, b_o\in\mathbb{R}\) are learnable parameters shared across all product–segment edges. Collecting these logits gives the product-by-segment logit matrix \(\boldsymbol{\Lambda} = [\Lambda_{jk}]_{j\in[n],\,k\in[m]} \in \mathbb{R}^{n\times m}.\)

We apply the sigmoid function entry-wise and transpose the result to obtain the segment-product probability matrix \[\mathbf{P} = \sigma(\boldsymbol{\Lambda})^\top \in (0,1)^{m\times n},\] where \(\sigma(x)=1/(1+\exp(-x))\). Thus, \(\mathbf{P}_{kj} = \sigma(\Lambda_{jk})\) is the predicted probability that product \(j\) is included in the bundle purchased by segment \(k\).

4.4 Training Process↩︎

To train the GNN, we first collect historical problem instances, and then use the Mixed Bundling policy proposed by [1] to obtain the corresponding optimal bundle assignments \(\theta_{kb}^{\star}\)’s. We then project these bundle assignments to segment-product binary labels where \(q_{kj}^{\star}=1\) if product \(j\) is included in the optimal bundle assigned to segment \(k\), and \(q_{kj}^{\star}=0\) otherwise, i.e., \(q_{kj}^{\star}:=\sum_{\substack{b: j\in \mathcal{A}(b)}} \theta_{kb}^{\star}\). Specifically, our training procedure adopts a supervised learning paradigm, where the GNN is trained to predict selection probabilities for each customer-product pair, serving as a probabilistic predictor of the underlying binary product-selection decisions. Importantly, the model does not learn bundle prices directly; instead, it learns to approximate the optimal bundle assignment structure implied by the Mixed Bundling solution, which is leveraged to conduct high-quality search space pruning before solving the MILP, thereby accelerating the solution process by reducing computational complexity while maintaining near-optimal solution quality.

Thus, each instance \(\ell\) is characterized by \((\mathbf{c}_{(\ell)}^u, \mathbf{c}_{(\ell)}^s, \mathbf{U}_{(\ell)}, \boldsymbol{\alpha}_{(\ell)}, f(\cdot), \mathbf{Q}^{\star}_{(\ell)})\), where \(\mathbf{U}_{(\ell)} = [u_{kj,(\ell)}]_{m \times n}\) and \(\mathbf{Q}^{\star}_{(\ell)} = [Q^{\star}_{kj,(\ell)}]_{m \times n}\). The instances were split into a training set (80%) and a validation set (20%). We then convert each problem instance into the corresponding graph representation and construct the input tensors \(\mathbf{Y}^P\) and \(\mathbf{Y}^S\). We formulate the edge prediction task as a binary classification problem. During each training epoch, we first pass the graphs in the training set to the Graph Neural Network to obtain the predicted product selection probability \(\mathbf{P}_{kj}\)’s. We compute the weighted binary cross-entropy loss \(\mathcal{L}\). Let \(S_T\) denote the number of instances in the training set. To address the inherent class imbalance (as optimal bundles typically contain only a small subset of available products), we introduce a positive weight parameter \(w_{pos} = (1 - \kappa) / \kappa\), where \(\kappa = \frac{\sum_{\ell=1}^{S_T}\sum_{k=1}^{m_{(\ell)}}\sum_{j=1}^{n_{(\ell)}} q_{kj,(\ell)}^{\star}}{\sum_{\ell=1}^{S_T} m_{(\ell)}n_{(\ell)}}\) is the proportion of positive edges across the entire training set. The loss function is defined as: \[\mathcal{L} = -\frac{1}{ \sum_{\ell=1}^{S_T}m_{(\ell)}n_{(\ell)} } \sum_{\ell=1}^{S_T} \sum_{k=1}^{m_{(\ell)}} \sum_{j=1}^{n_{(\ell)}} \left[ w_{\mathrm{pos}}q_{kj,(\ell)}^{\star} \log \mathbf{P}_{kj,(\ell)} + (1-q_{kj,(\ell)}^{\star}) \log(1-\mathbf{P}_{kj,(\ell)}) \right].\] The AdamW optimizer is used to update the model weights to minimize the training loss. The resulting model is then applied to the inference policies described in Section 4.5.

4.5 Pruning-based Inference Policies↩︎

In this subsection, we leverage the probability matrix \(\mathbf{P}\) predicted by the GNN model to guide two pruning strategies that construct a compact yet high-quality bundle space. The resulting reduced problem is then solved using the MILP formulation of [1], but restricted to the pruned bundle set rather than the full exponential space.

4.5.0.1 Fixed Cutoff Pruning (FCP)

Our first policy generates one candidate bundle for each customer segment. For each segment \(k \in [m] = \{1,\dots,m\}\) and \(j \in [n] = \{1,\dots,n\}\), we define its candidate bundle by a fixed cutoff value. In our default method, we set the fixed cutoff to 0.5. (Although the choice of the cutoff threshold can be tuned based on performance, we provide a theoretical justification for this thresholding rule in Section 5. In addition, we show a sensitivity analysis over a wide range of cutoff values, showing that 0.5 offers a robust profit–runtime trade-off.) That is, \[B_k^{\mathrm{FCP}} = \{j\in[n]\mid \mathbf{P}_{kj}\ge 0.5\}, \qquad k\in[m],\] and the overall candidate bundle family is \(\mathfrak B^{\mathrm{FCP}} = \{B_k^{\mathrm{FCP}}:k\in[m]\}\cup\{\emptyset\} \subseteq\mathfrak F.\) Here, the overall candidate family is shared across all segments, meaning that any segment is free to select any \(B_k^{\mathrm{FCP}}\), not only its own. This reduces the feasible space from \(2^n\) to \(O(m)\), making the subsequent MILP much more tractable. We then solve \(\Xi(\mathfrak B^{\mathrm{FCP}})\).

In the case where all probabilities of a segment fall below the cutoff, we retain the single product with the highest probability to avoid producing an empty bundle. The detailed procedure is given in Algorithm [alg:fcp] in Appendix 14.

4.5.0.2 Progressive Cutoff Pruning (PCP).

The PCP policy is a more conservative pruning rule. For each segment \(k\), we first retain products whose predicted probabilities exceed the cutoff \(U_k=\{j\in[n]:\mathbf{P}_{kj}\ge 0.5\},\) and sort them in descending order of \(\mathbf{P}_{kj}\): \(U_k=(j_{k,(1)},j_{k,(2)},\dots,j_{k,(|U_k|)}).\) We then construct a sequence of prefix bundles \[B_{k,i}^{\mathrm{PCP}} = \{j_{k,(1)},j_{k,(2)},\dots,j_{k,(i)}\}, \qquad i=1,\dots,|U_k|.\] The resulting candidate family is \(\mathfrak B^{\mathrm{PCP}} = \bigcup_{k=1}^m \{B_{k,i}^{\mathrm{PCP}}:i=1,\dots,|U_k|\} \cup\{\emptyset\} \subseteq \mathfrak F.\) We then solve \(\Xi(\mathfrak B^{\mathrm{PCP}})\) over the associated index set \(\mathcal{I}_{\mathfrak B^{\mathrm{PCP}}}\).

For PCP, the retained prefix bundles lead to a richer candidate family and therefore more cover-form price-subadditivity constraints than FCP. We handle these constraints using a cutting-plane implementation, with details provided in Appendix 15.

Figure 4 illustrates the construction of candidate bundles under the FCP and PCP policies for a scenario with \(m_{\mathrm{test}}=1\) and \(n_{\mathrm{test}}=5\). First, products 2 and 4 are excluded from consideration because their predicted probabilities fall below the cutoff \(\tau=0.5\). Under the FCP policy, the remaining high-probability products form a single candidate bundle, \(B_1^{\mathrm{FCP}}=\{1,3,5\}.\) Thus, the resulting candidate bundle family is \(\mathfrak B^{\mathrm{FCP}} = \{\emptyset,\{1,3,5\}\}.\) In contrast, under the PCP policy, given the probability order \(\mathbf{P}_{11}>\mathbf{P}_{15}>\mathbf{P}_{13},\) the products are added sequentially to form a nested chain. Consequently, the segment-specific PCP candidate family is \(\mathfrak B_1^{\mathrm{PCP}} = \bigl\{\{1\},\{1,5\},\{1,5,3\}\bigr\},\) and the global candidate bundle family is \(\mathfrak B^{\mathrm{PCP}} = \{\emptyset\}\cup\mathfrak B_1^{\mathrm{PCP}}.\)

4.5.0.3 Pruning Effect on Model Size.

The full mixed-bundling formulation optimizes over the entire bundle family \(\mathfrak F=2^{[n]}\), leading to \(O(m2^n)\) assignment-related variables and self-selection constraints, together with exponentially many price-subadditivity constraints. FCP provides the strongest compression by retaining at most one candidate bundle per segment, so \(|\mathfrak B^{\mathrm{FCP}}|\le m+1\). PCP is more conservative: it retains a nested prefix chain for each segment, so \(|\mathfrak B^{\mathrm{PCP}}|\le mn+1\), and uses cutting planes to avoid pre-enumerating all cover-subadditivity constraints. Thus, both policies replace the exponential bundle universe with polynomial-size restricted families; FCP prioritizes efficiency, while PCP retains more pricing flexibility.

4.6 GNN-Guided Local Search Policy↩︎

In order to further improve our solution, we propose an inference policy based on GNN-guided local search. The basic idea is to iteratively modify the bundle assignment by either adding an unselected product or dropping a selected product, and accept modifications if the profit increases.

To maximize the efficiency of the local search, we develop a global top-k local search strategy guided by the probability matrix \(\mathbf{P}\) predicted by the GNN model. Instead of evaluating neighbors for all segments, we focus only on the most promising modifications globally and prioritize candidates based on a confidence score. Specifically, for any segment \(k\) and product \(j\), we define the confidence score for a potential modification as: \[\text{Score}_{kj} = \begin{cases} \mathbf{P}_{kj} & \text{if product } j \text{ is currently unselected (Candidate for Addition)} \\ 1 - \mathbf{P}_{kj} & \text{if product } j \text{ is currently selected (Candidate for Removal)} \end{cases}\] Intuitively, a higher \(\mathbf{P}_{kj}\) indicates a strong suggestion to include the product, while a lower \(\mathbf{P}_{kj}\) implies a strong suggestion to exclude it.

In each iteration, we select the top \(K\) candidate additions and the top \(K\) candidate removals according to the GNN confidence scores. The selected moves are pooled and evaluated in descending order of confidence; the first improving move is accepted, and the search restarts from the updated assignment. We use the sublinear budget \(K=\lceil\sqrt m\rceil\), with supporting sensitivity analysis reported in Appendix 12.

The local search process terminates when no improvement is found among the top candidates in a full cycle, or when the predefined iteration limit is reached. Furthermore, we rely on the LP relaxation to efficiently evaluate each neighbor, as it provides a valid lower bound on the restricted MILP optimum and significantly reduces computational cost. We demonstrate a high degree of consistency between LP and MILP objective improvements during the local search process, justifying the use of LP as a reliable proxy in Appendix 13. We use FCPLS to denote this algorithm.

4.7 Iterative Self-Improvement↩︎

The preceding policies exploit a key advantage of GNNs: a model trained on small instances can be applied to larger instances because its parameters are shared across nodes and edges. This provides useful scalability and enables small-to-large transfer. Nevertheless, this scalability can be limited. When the product dimension increases substantially, the structure of the optimal bundle family becomes richer, and a model trained only on small exact labels may be less well calibrated for larger product spaces. As a result, its predicted segment-product probabilities may become less accurate, leading to lower-quality candidate families and weaker inference-policy performance.

A natural solution would be to train directly on large-scale optimal labels. However, it is computationally infeasible to generate one exact label by solving the full mixed-bundling formulation \(\Xi(\mathfrak F)\), whose bundle space grows exponentially in \(n\). We therefore propose an iterative self-improvement strategy that improves scalability stage by stage. In particular, starting from a base GNN trained on small exact labels, we use its predictions to run PCP on larger unlabeled instances. The proposed policies are fast and provide high-quality solutions for retraining. The retrained GNN can then be used to generate solutions for even larger instances, so the procedure can in principle be repeated over multiple stages. In this sense, the method is a bootstrapping scheme in which the current learned policy generates improved large-scale supervision for the next model. An illustration of this method is shown in Figure 5, while the detailed algorithm is given in Appendix 14. In the numerical experiments, we report the effect of one self-improvement round to show the effect of this step. The numerical results in Section 6.2 show that this self-improved model substantially improves transfer to larger product spaces while preserving the computational efficiency of the pruning framework.

5 Theoretical Foundation↩︎

In this section, we provide a theoretical justification for the edge-level prediction target used by our GNN-guided pruning framework. Recall that the implemented network in Section 4.3 outputs a segment-product probability matrix \(\mathbf{P} \in[0,1]^{m\times n}\), where \(\mathbf{P}_{kj}\) is interpreted as the probability that product \(j\) is included in the bundle selected by segment \(k\). The inference policies in Section 4.5 use this matrix to construct restricted candidate bundle families before solving the mixed-bundling MILP.

Our GNN does not attempt to represent the full mixed-bundling solution, which contains exponentially many bundle-level prices and choice variables. Instead, we focus on the product-assignment projection of an optimal mixed-bundling solution. Let \(W^\star=(\boldsymbol{\Theta}^\star,\mathbf{p}^\star, \mathbf{P}^{\mathrm{pay},\star},\mathbf{s}^\star,\mathbf{Z}^\star)\) be an optimal solution of the full formulation \(\Xi(\mathfrak F)\), where \(\mathfrak F=2^{[n]}\). Let \(\mathbf{M}\in\{0,1\}^{(2^n)\times n}\) denote the bundle-product incidence matrix with \(m_{bj}:=\mathbb{1}_{\{j\in \mathcal{A}(b)\}}.\) The product-assignment projection is \[\mathbf{Q}^{\star}=\boldsymbol{\Theta}^{\star} \mathbf{M}\in\{0,1\}^{m\times n}, \qquad Q^{\star}_{kj} = \sum_{b:j\in \mathcal{A}(b)}\theta^{\star}_{kb}.\] Thus, \(Q^{\star}_{kj}=1\) if and only if product \(j\) is contained in the bundle selected by segment \(k\) in the optimal mixed-bundling solution. This projection is naturally an edge-level object on the compact customer-product graph since each coordinate corresponds to one customer-product pair.

In this section, for ease of theoretical analysis, we generalize the network in Section 4.3 and consider a class \(\mathcal{F}_{\mathrm{edge}}^{W,\sigma}\) of finite-layer GNN. We first consider a class \(\mathcal{F}_{\mathrm{edge}}^{W}\) of finite-layer GNN without the last sigmoid layer.

Let \(X=(G,H,E)\in\mathcal{X}_{m,n}\) denote the customer-product graph input, where the formal input space \(\mathcal{X}_{m,n}\) and the permutation action are defined in Appendix ¿sec:sec:app:preliminaries?. Let \(\mathbf{y}_k^S\), \(\mathbf{y}_j^P\), and \(u_{kj}\) denote the raw customer, product, and customer-product edge features defined in Section 4.2. We choose learnable continuous maps \(\iota_S,\iota_P,\iota_E, \{\phi_S^{(\ell)},\phi_P^{(\ell)},\Gamma_S^{(\ell)}, \Gamma_P^{(\ell)},\Gamma_E^{(\ell)}\}_{\ell=1}^{\omega}, R_E, R_1,\) which are shared across customers, products, and edges of the same type.

The logit-output GNN class first initializes \[\mathbf{v}_k^{S,(0)}=\iota_S(\mathbf{y}_k^S), \qquad \mathbf{v}_j^{P,(0)}=\iota_P(\mathbf{y}_j^P), \qquad \mathbf{e}_{jk}^{(0)} = \iota_E(\mathbf{y}_j^P,\mathbf{y}_k^S,u_{kj}).\] At each layer \(\ell=1,\ldots,\omega\), the customer and product embeddings are updated by shared message-passing maps: \[\mathbf{v}_k^{S,(\ell)} = \Gamma_S^{(\ell)} \left( \mathbf{v}_k^{S,(\ell-1)}, \sum_{j=1}^n u_{kj}\, \phi_P^{(\ell)}(\mathbf{v}_j^{P,(\ell-1)}) \right), \qquad \mathbf{v}_j^{P,(\ell)} = \Gamma_P^{(\ell)} \left( \mathbf{v}_j^{P,(\ell-1)}, \sum_{k=1}^m u_{kj}\, \phi_S^{(\ell)}(\mathbf{v}_k^{S,(\ell-1)}) \right).\]

The edge embedding is then updated using the current endpoint embeddings and the previous edge embedding: \[\mathbf{z}_{jk}^{(\ell)} := \mathbf{v}_j^{P,(\ell)} \Vert \mathbf{v}_k^{S,(\ell)} \Vert \mathbf{e}_{jk}^{(\ell-1)}, \qquad \mathbf{e}_{jk}^{(\ell)} = \Gamma_E^{(\ell)}\bigl(\mathbf{z}_{jk}^{(\ell)}\bigr),\] where \(\Vert\) denotes concatenation. Let \(\mathcal{U}^{(\omega)}(X)\) collect all terminal customer, product, and edge embeddings. The logit output is defined by \(F(X)_{kj} = R_E\bigl(\mathbf{e}_{jk}^{(\omega)},\mathcal{U}^{(\omega)}(X)\bigr), k\in[m],\;j\in[n],\) where the edge readout \(R_E\) is invariant in \(\mathcal{U}^{(\omega)}(X)\). Each \(F\in\mathcal{F}_{\mathrm{edge}}^{W}\) maps an instance \(X\) to a real-valued matrix \(F(X)\in\mathbb{R}^{m\times n}.\) The corresponding scalar-output class is denoted by \(\mathcal{F}_{\mathrm{edge}}^1\), where a scalar output has the form \(f(X)=R_1(\mathcal{U}^{(\omega)}(X))\) with a permutation-invariant scalar readout \(R_1\).

Then, we define the corresponding class \(\mathcal{F}_{\mathrm{edge}}^{W,\sigma}\) of finite layer GNNs with a sigmoid layer as \[\mathcal{F}_{\mathrm{edge}}^{W,\sigma} = \{\, \sigma \circ \Phi : \Phi \in \mathcal{F}_{\mathrm{edge}}^{W} \,\},\] where the sigmoid function \(\sigma(\cdot)\) is applied entry-wise.

The main representation result below will be stated for a canonical single-valued optimal assignment mapping \(\Phi_Q\), which is constructed in Appendix ¿sec:sec:subsec:opt95assign?. Before stating that result, we first introduce the WL refinement used to characterize the expressive power of \(\mathcal{F}_{\mathrm{edge}}^{W}\).

5.1 Customer-Product WL Refinement and Unfoldability↩︎

In this subsection, following the idea of [31], given a finite set of instances, we establish that GNN can approximate the optimal product assignment to full mixed-bundling problem \(\Xi(\mathfrak F)\) under some conditions.

We now describe the color-refinement procedure used to analyze the expressive power of \(\mathcal{F}_{\mathrm{edge}}^{W}\). The Weisfeiler–Lehman (WL) test [42] iteratively assigns colors to nodes based on their own colors and the multiset of neighbor colors. Different from the node-output WL test in [31], our network maintains edge embeddings and outputs an edge-level matrix. We therefore use a customer-product WL refinement that augments node-side colors with edge-side colors. This modification allows us to track which customer-product edge coordinates are distinguishable by the edge-output architecture.

We next introduce the notion of foldability. Intuitively, a foldable instance contains two customer nodes or two product nodes that remain indistinguishable under all rounds of WL refinement. Such nodes cannot be assigned different representations by a message-passing GNN that respects the graph symmetry.

For an unfoldable instance, the node-side WL colors become discrete after at most \(m+n\) rounds of refinement, where all customer-side colors are pairwise distinct and all product-side colors are pairwise distinct. Since the edge hash in Algorithm 6 is injective in the current endpoint colors, unfoldability also implies that distinct edge coordinates are distinguishable after finitely many rounds.

Finally, we define the cross-instance WL equivalence relation used in the separation theorem. This relation compares two customer-product graphs up to a simultaneous relabeling of customer and product indices.

The next subsection uses this relation to connect the separation power (ability to distinguish non-isomorphic graphs) of \(\mathcal{F}_{\mathrm{edge}}^1\) and \(\mathcal{F}_{\mathrm{edge}}^W\) to \(\mathrm{WL}_{\mathrm{CP}}\), and then states the representation theorem for the canonical optimal assignment mapping.

5.2 Representation result↩︎

We now define the target mapping represented by \(\mathcal{F}_{\mathrm{edge}}^{W,\sigma}\). Throughout this subsection, the target is defined with respect to the full mixed-bundling formulation \(\Xi(\mathfrak F)\), where \(\mathfrak F=2^{[n]}\) is the full bundle family. For an instance \(X\in\mathcal{X}_{m,n}\), let \(\Phi_{\mathrm{feas}}(X) := \mathbb{1}\{\Xi(\mathfrak F;X)\text{ is feasible}\}, \mathcal{D}_{\mathrm{solu}} := \mathcal{X}_{m,n}\cap \Phi_{\mathrm{feas}}^{-1}(1).\) Here \(\Xi(\mathfrak F;X)\) denotes the full mixed-bundling formulation induced by instance \(X\). By the boundedness argument in Lemma 2, every instance in \(\mathcal{D}_{\mathrm{solu}}\) admits an optimal solution.

For a feasible solution \(W=(\boldsymbol{\Theta},\mathbf{p}, \mathbf{P}^{\mathrm{pay}},\mathbf{s},\mathbf{Z})\) of \(\Xi(\mathfrak F;X)\), define its product-assignment projection by \(\mathbf{Q}(W):=\boldsymbol{\Theta}(W)\mathbf{M},\) where \(\mathbf{M}\in\{0,1\}^{|\mathcal{K}|\times n}\) is the bundle-product incidence matrix with entries \(m_{bj}=\mathbb{1}_{\{j\in \mathcal{A}(b)\}}, b\in\mathcal{K},\;j\in[n].\)

Since \(\Xi(\mathfrak F;X)\) may have multiple optimal solutions, the optimal product-assignment projection may not be unique. We therefore define a single-valued canonical optimal assignment mapping \(\Phi_Q: \mathcal{D}_{\mathrm{solu}}\setminus\mathcal{D}_{\mathrm{foldable}} \to \{0,1\}^{m\times n}.\) Formally, Appendix ¿sec:sec:subsec:opt95assign? constructs \(\Phi_Q\) by first applying a WL-based canonical sorting map \(\Phi_{\mathrm{sort}}\) on non-foldable instances and then selecting the lexicographically smallest matrix among all optimal product-assignment projections. This construction makes \(\Phi_Q\) well defined and permutation equivariant.

We state the separation result connecting this edge-output GNN class to the WL refinement on the customer-product graph.

Theorem 1 shows that \(\mathcal{F}_{\mathrm{edge}}^{W}\) has the same separation power as the WL refinement.

A direct implication of Theorem 2 is that the empirical binary cross-entropy loss \(\mathcal{L}_D(F_\varepsilon^\sigma)\) can be made arbitrarily small on the finite dataset \(D\): \(\mathcal{L}_D(F_\varepsilon^\sigma) \le |D|mn[-\log(1-\varepsilon)]\to0 \quad\text{as }\varepsilon\downarrow0.\)

6 Numerical Experiments↩︎

We now present comprehensive numerical experiments designed to evaluate the effectiveness and scalability of our proposed GNN-guided bundle pricing strategies. The goals of this section are twofold: (i) to demonstrate that our approaches consistently achieve high profit ratios while maintaining significant computational efficiency, and (ii) to validate their robustness across heterogeneous problem settings with varying numbers of customer segments and products.

To thoroughly test scalability and robustness, we evaluate our strategies under varying numbers of customer segment \(m\) and product \(n\). For each scenario, results are averaged over 100 instances. The GNN model training was conducted with a single NVIDIA A100 GPU, while the GNN model inference and the inference policies were performed on a local machine featuring an Apple M1 Pro CPU (3.2 GHz) and 16GB RAM, utilizing up to 8 threads. The implementation uses Gurobi 12.0, Python 3.10, PyTorch 2.8, and CUDA 12.8.

To evaluate any method, we adopt two normalized metrics: profit ratio (PR) and time ratio (TR) \[PR_{A,B} = \tfrac{\text{Profit of policy A}}{\text{Profit of policy B}},\quad TR_{A,B} = \tfrac{\text{Runtime of policy A}}{\text{Runtime of policy B}},\] where PR measures solution quality and TR captures efficiency. These two metrics allow direct comparison of solution quality and speed across different algorithms.

6.1 Base Simulations↩︎

We compare our proposed approach against the two most fundamental and widely recognized policies in the literature: the exact Mixed Bundling (MB) policy and the approximation-based Bundle Size Pricing (BSP) policy:

1. Mixed Bundling (MB) policy. For MB, we follow [1]’s MILP formulation. The MB formulation is a MILP that assigns a separate price to each candidate bundle. Detailed formulation is provided in Section 3.1. We solve the MB policy with a MIPGap threshold of 0.1%.

2. Bundle Size Pricing (BSP) policy. The BSP baseline is proposed by [2]. The BSP approximates MB by assigning the same price to all bundles of equal size. Detailed formulation is provided in Appendix 9.1. We solve the BSP policy with an MIPGap threshold of 0.1%.

It is worth noting that we exclude certain heuristic algorithms and general neural-network-based MILP methods from our baselines due to their inapplicability to the general bundle pricing problem, the detailed discussion is provided in Appendix 16.1.

The instance-generation procedure and GNN training details are provided in Appendix 16.2. In the main experiments, we train the base GNN on 3000 small instances with \((m_{\mathrm{train}},n_{\mathrm{train}})=(10,10)\), where labels are obtained by solving the exact MB formulation and projecting the optimal assignments to segment-product labels. The dataset is split into training and validation sets, and the best checkpoint is selected by validation loss. To mitigate random initialization effects, we train 10 GNN models with seeds 1 to 10 (see Appendix 16.3 for details) and report the average performance of the resulting inference policies. All MILP-based policies are solved with a MIPGap threshold of \(0.1\%\).

In Table 2, we report the numerical results of our algorithms compared with both baselines for problems with \(n=10\) and varying numbers of customer segments. All three GNN-based policies recover more than \(98\%\) of the optimal MB profit while requiring only a small fraction of the MB runtime. In contrast, BSP is computationally efficient but incurs a substantial profit loss, achieving only about \(85\%\)–\(88\%\) of the MB benchmark. Among the proposed policies, FCP provides the fastest policy, PCP achieves the highest profit by retaining a richer candidate family, and FCPLS offers an intermediate trade-off by improving upon FCP through local search while remaining faster than PCP.

For problems with more than 10 products, solving the full MB formulation becomes computationally challenging. Therefore, we use BSP as the baseline for evaluating larger instances. The results are reported in Table 3.

Overall, FCP achieves strong performance with high computational efficiency, consistently outperforming BSP while using only a fraction of its runtime. PCP further improves the profit ratio by retaining a richer candidate bundle family, but its runtime grows substantially; hence it is not evaluated for \((m_{\mathrm{test}}=20,n_{\mathrm{test}}=40)\). FCPLS provides a middle ground by improving upon FCP through GNN-guided local search while remaining more scalable than PCP on larger instances.

6.1.0.1 Scalability across the number of customer segments \(m_{\mathrm{test}}\).

Under a 600-second time limit, Figure 7 illustrates the scalability of the FCP policy as the number of customer segments \(m_{\mathrm{test}}\) varies while fixing \(n_{\mathrm{test}}=20\). The FCP policy consistently outperforms BSP in terms of profit across all tested values of \(m_{\mathrm{test}}\), achieving profit ratios above 1.16. In terms of runtime, FCP remains substantially faster than BSP for all sizes of \(m_{\mathrm{test}}\). Overall, these results show that FCP scales effectively with the number of customer segments while maintaining strong solution quality.

6.1.0.2 Scalability across the number of products \(n_{\mathrm{test}}\).

Under a 600-second time limit, Figure 8 illustrates the scalability of FCP and PCP as the number of products \(n_{\mathrm{test}}\) varies while fixing \(m_{\mathrm{test}}=10\). PCP demonstrates strong solution quality, consistently achieving profit ratios above 1.16 relative to BSP and improving as \(n_{\mathrm{test}}\) increases. This benefit, however, comes with higher computational cost because PCP retains a larger candidate bundle family. By contrast, FCP remains extremely fast across all tested product sizes, but its solution quality decreases as \(n_{\mathrm{test}}\) grows. In particular, when \(n_{\mathrm{test}}=50\), FCP becomes comparable to BSP in terms of profit.

These results reveal a useful but incomplete form of size generalization. The GNN trained on \(n_{\mathrm{train}}=10\) captures transferable customer-product patterns, but a large increase in the product dimension changes the combinatorial structure of the optimal bundle family. Rather than treating this as a limitation of the GNN framework, we use it as an opportunity for self-improvement: the small-scale model can still generate high-quality large-scale solutions through PCP, and these solutions can in turn serve as near-optimal labels for retraining at the larger scale.

6.2 Iterative Self-Improvement↩︎

To evaluate the efficacy of the iterative self-improvement strategy proposed in Section 4.7, we conduct a specific experiment to assess cross-scale generalization.

In the numerical experiments, we first train a GNN model using instances of size \(m_{\mathrm{train}}=10, n_{\mathrm{train}}=10\) following the setting in Section 6.1. Then, we generate 3000 near-optimal instances with size \(m_{\mathrm{train}}=10\) and \(n_{\mathrm{train}}=50\) generated by PCP policy, and we train a separate GNN model using these 3000 samples with \(m_{\mathrm{train}}=10\) and \(n_{\mathrm{train}}=50\). The corresponding inference policies using the new model are denoted by FCP-I, PCP-I, and FCPLS-I.

Figure 9 illustrates the effect of the iterative self-improvement strategy. The self-improved model trained on \(n_{\mathrm{train}}=50\) improves size generalization, especially for problems with larger product size. For FCP, FCP-I is slightly worse than the base model at \(n_{\mathrm{test}}=20\), but outperforms it for all \(n_{\mathrm{test}}\ge30\), with larger gains as \(n_{\mathrm{test}}\) increases. For PCP, PCP-I is worse than the base PCP model at \(n_{\mathrm{test}}=20\), but becomes comparable at \(n_{\mathrm{test}}=30\) and outperforms the base model for all \(n_{\mathrm{test}}\ge 40\), with the largest gain observed at \(n_{\mathrm{test}}=50\).

Runtime remains comparable between FCP and FCP-I, while PCP-I is generally faster than the base PCP model despite achieving stronger solution quality. These results show that iterative self-improvement improves transfer to larger product spaces while preserving the computational efficiency of the pruning framework.

Table 4 reports the performance of the self-improved policies on larger instances up to \(n_{\mathrm{test}}=100\). FCP-I is the most scalable policy: it consistently outperforms BSP in profit, with profit ratios between about \(1.13\) and \(1.21\), while requiring only a small fraction of BSP’s runtime. FCPLS-I provides a moderate but consistent quality improvement over FCP-I, at the cost of additional local-search time. PCP-I achieves the highest profit ratios in the reported \(m_{\mathrm{test}}=10\) settings, but its runtime grows quickly with \(n_{\mathrm{test}}\), so it is not reported for the larger \(m_{\mathrm{test}}=20\) settings.

7 Extension: Additive Random Valuations↩︎

The main experiments focus on the deterministic segment-level mixed-bundling setting studied in Sections 3–6. In that setting, each instance contains a fixed set of customer segments with known bundle valuations, and the objective of \(\Xi(\mathfrak F)\) is the exact profit for that deterministic population. In this section, we conduct an external validation experiment under additive random valuations. Under additive valuations, the value of bundle \(b\) for customer \(i\) is \(R_{ib} = \sum_{j\in \mathcal{A}(b)} v_{ij}.\)

In a random-valuation problem, customer valuations are drawn from an underlying distribution. The seller’s objective is therefore an expected profit, rather than the profit from a fixed deterministic set of segments. A common approach for solving such stochastic pricing problems is sample average approximation (SAA): draw a finite sample of customers, solve the deterministic optimization problem induced by this sample, and then evaluate the resulting policy on fresh out-of-sample valuation draws. We adopt this SAA approach here. Thus, the sampled customers play a role analogous to the deterministic customer segments in the main formulation, but the evaluation criterion is different: a good policy must perform well not only on the optimization sample but also out of sample.

This experiment serves two purposes. First, it tests whether the proposed GNN-guided pruning framework remains effective in a random-valuation setting when combined with an SAA formulation. Second, it evaluates whether the restricted mixed-bundling formulation solved by FCP can produce prices that generalize to unseen random valuations. For a sampled instance, FCP constructs a candidate family \(\mathfrak B^{\mathrm{FCP}}\) from the GNN predictions and solves the restricted formulation \(\Xi(\mathfrak B^{\mathrm{FCP}})\). We then compare its in-sample and out-of-sample performance with two structured benchmarks: bundle-size pricing (BSP) and CPBSD-A, an additive approximation of component pricing with bundle-size discounts based on [15] and described in Appendix 16.4.

In this section, we use \(N\) to denote the number of products, \(K_{\mathrm{in}}\) to denote the number of in-sample customers used in the SAA optimization problem, and \(K_{\mathrm{oos}}\) to denote the number of fresh out-of-sample valuation draws used for evaluation. We consider two product-scale blocks. The first uses \(N=10\) and \(K_{\mathrm{in}}=50\); the second uses \(N=30\) and \(K_{\mathrm{in}}=50\). In both cases, the resulting price menu is evaluated on \(K_{\mathrm{oos}}\) out-of-sample draws. For each scale, we test three product cost settings: zero costs, independent random costs, and costs correlated with valuation means. For FCP, we reuse the same GNN checkpoint trained once on smaller additive random-valuation instances, without any scale-specific or cost-specific retraining. We report in-sample profit (InS), out-of-sample profit (OOS), and runtime. Details on valuation generation, cost definitions, the CPBSD-A formulation, implementation settings, and the FCP evaluation procedure under SAA are provided in Appendix 16.4; the procedure is summarized in Algorithm [alg:saa95fcp95random95valuation].

Table 5 shows that FCP has the best overall performance across these additive random-valuation problems. It achieves the lowest runtime in every cost setup at both product scales and obtains the highest OOS profit in five out of six settings. The only exception is the zero-cost case at \(N=30,K_\mathrm{in}=50\), where BSP is marginally higher in OOS profit, while the difference is negligible relative to the runtime gap. Under the two random-cost settings, FCP attains the highest InS and OOS profit in all four rows, indicating that its advantage appears both in the SAA objective and in OOS generalization. In the simpler zero-cost cases, all three methods have close profit values: CPBSD-A achieves the highest InS profit at \(N=10,K_{\mathrm{in}}=50\), whereas FCP achieves the highest InS profit at \(N=30,K_{\mathrm{in}}=50\). The scalability contrast is especially clear at \(N=30,K_{\mathrm{in}}=50\): CPBSD-A reaches the 300-second time limit in all three cost setups, whereas FCP remains within a few seconds. These results support the main message that learned pruning works well under additive random-valuation problems while preserving both profitability and computational efficiency.

8 Conclusion↩︎

This paper develops a GNN-guided pruning framework for large-scale mixed bundle pricing. Instead of constructing and optimizing over the full exponential bundle space, we learn segment-product inclusion probabilities on a compact customer-product graph. These predictions are used to prune the search space and construct a much smaller candidate bundle family, over which we solve the mixed-bundling formulation.

Algorithmically, we propose two pruning policies, FCP and PCP, that trade off candidate-family size and solution quality. We further introduce a GNN-guided local-search procedure that improves solution quality from a given initial solution. To enhance scalability, we propose an iterative self-improvement scheme that uses GNN-generated high-quality large-scale solutions as near-optimal labels to improve subsequent GNN-based policies. Numerical experiments show that the proposed framework achieves near-optimal performance on tractable instances, scales to much larger product spaces, and outperforms structured heuristics. Additional validation under additive random valuations shows that the same learning-and-pruning idea can be embedded in an SAA pipeline and can produce prices that generalize well out of sample while maintaining a clear computational advantage.

Theoretically, we show that the optimal product-assignment projection is the right edge-level object for learning since it avoids the exponential bundle-level MILP graph while preserving the assignment structure needed for pruning. We further establish, via a customer-product WL refinement with edge colors, that a specific edge-output GNN can represent this canonical assignment target on finite solvable and non-foldable datasets and recover the binary assignment by thresholding.

References↩︎

---

1. https://www.ordergroove.com/blog/product-bundling/↩︎