Multi-Dimensional Knapsack, in depth: d-constraint DP, LP bounds, dual-price greedy and MILP solvers on a GPU packing example

The problem and what changes with d constraints

There are n items. Item i has value v_i and, for each of d resources k, a non-negative weight w_ik. Resource k has capacity C_k. Choose x_i in {0, 1} to maximise the sum of v_ix_i subject to, for every k, the sum of w_ikx_i being at most C_k. With d = 1 this is ordinary 0/1 knapsack. Written as a matrix, it is the integer program: maximise v^Tx subject to Wx ≤ C and x binary, where W has one row per resource.

Three facts shape every algorithmic choice:

• It is NP-hard for every d, because d = 1 already is. For d = 1 the hardness is only pseudo-polynomial: a DP over capacity runs in O(nC). For MKP the equivalent DP runs in O(n times the product of all C_k), which explodes as d grows.
• No FPTAS for d ≥ 2 unless P = NP. A PTAS exists for fixed d, but it is impractical.
• The LP relaxation is tight in structure. If you allow 0 ≤ x_i ≤ 1, a basic optimal solution of the LP has at most d fractional variables. Rounding those down gives a feasible integer solution that loses at most the value of d items. When items are small compared with capacities, that loss is small, which is why LP-based methods dominate in practice.

Worked example: packing jobs on a GPU node

A scheduler has one node with 8 GPUs, 256 GB of host memory and 48 CPU cores, and seven queued jobs. Each job has a priority value, and the scheduler wants the subset with the highest total priority that fits all three dimensions at once. Splitting a job is not allowed.

Brute force over all 128 subsets gives one optimum: {B, C, F, G} with value 26, using 7 GPUs, 240 GB and 36 cores. Notice that the two most valuable jobs, A and D, are both left out. Together they need 8 GPUs and 256 GB, so they fit those two dimensions exactly, but they leave no room for anything else and score only 19. Each dimension pushes the answer in a different direction, and the binding resource changes depending on which other jobs are in the set.

The LP relaxation (solved with HiGHS through SciPy) reaches 28.2, taking B, C and G fully, A at 0.3, E at 0.46 and F at 0.67, with all three resources used to the limit. Three fractional variables for three constraints is the most the structural result allows. The bound tells us no integer solution can beat 28, so 26 is within about 8 percent of the bound before we have even proved optimality.

Exact DP over several capacities

The exact DP generalises directly. Let best[c₁, ..., c_d] be the best value using at most c₁ to c_d of each resource. Process items one at a time and, as in one-dimensional 0/1 knapsack, extend only states from before the current item so each item is used at most once. A dictionary keyed by tuples of used amounts keeps only reachable states, which is often far fewer than the full grid.

def prune_dominated(states):
    # O(s^2) for clarity; sort by value and sweep for large state sets
    items = sorted(states.items(), key=lambda kv: -kv[1][0])
    kept = {}
    for used, entry in items:
        if not any(all(a <= b for a, b in zip(k, used)) for k in kept):
            kept[used] = entry      # nothing kept so far is cheaper everywhere and at least as good
    return kept

def mkp_dp(values, weights, caps):
    # weights[i] is a tuple of d resource amounts for item i
    # state: tuple of used amounts -> (best value, chosen items)
    states = {tuple(0 for _ in caps): (0, ())}
    for i, (v, w) in enumerate(zip(values, weights)):
        new_states = dict(states)
        for used, (val, items) in states.items():
            nxt = tuple(u + wk for u, wk in zip(used, w))
            if all(n <= cap for n, cap in zip(nxt, caps)):
                if nxt not in new_states or new_states[nxt][0] < val + v:
                    new_states[nxt] = (val + v, items + (i,))
        # dominance pruning: drop states that use more of every resource for no more value
        states = prune_dominated(new_states)
    return max(states.values())

The full grid for the example has 9 × 257 × 49 = 113,337 cells. Every RAM figure is a multiple of 16 GB, so measuring RAM in 16 GB units cuts it to 9 × 17 × 49 = 7,497. Dividing each dimension by the greatest common divisor of its weights is exact; rounding to coarser units is not. Dominance pruning helps further: a state using at least as much of every resource for no more value can never win, so keeping only the Pareto frontier can shrink the state set by orders of magnitude.

Greedy scores: collapsing d weights into one

With one constraint, sorting by v_i/w_i is the natural greedy. With d constraints you must first collapse d weights into one number, and the way you collapse them decides the answer. Here are four choices, each followed by first-fit in score order, on the worked example:

Ignoring memory and cores does worst. The best score uses dual prices: the LP's shadow price y_k says what one more unit of resource k is worth. HiGHS reports 1.6 per GPU, 0.0375 per GB and 0 per core, so cores are free at the margin, and greedy ranked this way finds the optimum of 26.

Do not over-read that last row. Jobs A, E and F all score exactly 1.0 under the dual prices, because items with fractional LP values sit on the price boundary, so a different tie-break could change the result. Treat dual-price greedy as a strong starting solution, improve it with swap-based local search, and compare it with the LP bound. Folding constraints into the objective with multipliers is Lagrangian or surrogate relaxation; subgradient updates refine the multipliers when the LP itself is too large.

Solving it with a MILP solver

For most production instances, the right answer is a mixed-integer solver. MKP is a textbook integer program, and modern branch-and-cut solvers combine the LP bound, cutting planes and primal heuristics far better than hand-written code. SciPy ships an interface to HiGHS, so you can solve the example without new dependencies:

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds

v = np.array([10, 7, 6, 9, 4, 8, 5])
W = np.array([
    [4,  2,  1,   4,  1,   2,  2],     # GPUs
    [96, 64, 32, 160, 64, 128, 16],    # RAM GB
    [12, 16, 8,  16,  24,  8,  4],     # CPU cores
])
C = np.array([8, 256, 48])

res = milp(
    c=-v,                                        # milp minimises
    constraints=LinearConstraint(W, -np.inf, C),
    integrality=np.ones(len(v)),                 # all variables integer
    bounds=Bounds(0, 1),                         # ... and binary
    options={"time_limit": 5.0, "mip_rel_gap": 0.0},
)
chosen = np.flatnonzero(res.x > 0.5)            # [1, 2, 5, 6] -> B, C, F, G
print(-res.fun, chosen)                          # 26.0 (up to float rounding)

Read the solution with a threshold such as 0.5, because solvers return values like 0.9999999, and always set a time limit and log the status and gap: a run stopped by the limit returns a feasible but possibly non-optimal answer.

Branch and bound for MKP follows the general method in Integer Programming and LP relaxation: solve the relaxation, branch on a fractional variable, and prune any node whose bound cannot beat the incumbent. Unlike single-constraint knapsack, the relaxation is a real LP rather than the fractional greedy from Fractional Knapsack, which only works with one constraint.

Choosing a method by instance size

Choosing a method is mostly about n, d and the size of the capacities after scaling:

Online admission is how schedulers usually use MKP: solve the LP over a recent window for prices y_k, then admit each arriving job that fits and whose value exceeds the sum of y_kw_ik. Scarce resources get high prices, so the rule protects the bottleneck.

Failure modes

Most MKP bugs come from modelling, not from the algorithm:

• A forgotten dimension. The first greedy above ignores two resources and loses value. In production, the same mistake overcommits a node, which looks like a stability bug rather than an optimisation bug. List every resource that can be exhausted, including those that are cheap now.
• Rounding that breaks feasibility. Converting 12.5 GB into 13 units, or 12, changes the problem. Round weights up and capacities down if you must round, so every solution stays feasible.
• Trusting a heuristic without a bound. A greedy answer can look excellent and still be far from optimal on some instances. Log the LP bound next to every heuristic answer and alert on the gap.
• Equality where you meant a threshold. Reading x_i == 1 from a solver drops items returned as 0.9999999.
• Ignoring time limits. A MILP that solves in 50 ms on the test set can take minutes on an adversarial instance. Without a limit, the scheduler stalls.

Trade-offs

The central trade-off is optimality against time and predictability. The DP gives exact answers with memory you can compute in advance, but only for small d and capacities. A solver gives exact or near-exact answers with a reported gap, but its running time varies widely between instances and it adds a dependency. Greedy with dual prices is fast, explainable (each item has a score you can show) and easy to run online, but carries no guarantee. The LP bound connects them, because it costs one LP solve and tells you how much any method is leaving on the table.

MKP also assumes values add up and items are independent. Pairing or exclusion rules between jobs are one extra constraint each in a solver, which is another reason to prefer one. If an item can be taken several times, see Bounded Knapsack.

What to do next

1. Write down every resource your items consume and the capacity of each, then check whether any dimension never binds; drop it only after measuring.
2. Scale each dimension by the greatest common divisor of its weights, and compute the product of the scaled capacities to see whether an exact DP is affordable.
3. Solve the LP relaxation and record the bound and the dual prices for your typical instance.
4. Implement dual-price greedy plus one-swap local search and compare its value with the LP bound on a week of real instances.
5. If the gap matters, model the problem with a MILP solver, set a time limit, warm-start from greedy, and log the solver's status and final gap.
6. Add tests that brute-force small random instances and check that every method returns a feasible set and that exact methods match the brute force.