diff --git a/Project.toml b/Project.toml
index bdec9be..868b4bc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,12 +8,15 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [weakdeps]
+AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
 
 [extensions]
+AbstractGPsDisjunctiveProgramming = "AbstractGPs"
 InfiniteDisjunctiveProgramming = "InfiniteOpt"
 
 [compat]
+AbstractGPs = "0.5"
 Aqua = "0.8"
 JuMP = "1.18"
 Reexport = "1"
@@ -30,4 +33,4 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
 
 [targets]
-test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt"]
+test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt", "AbstractGPs"]
diff --git a/README.md b/README.md
index 543d1f7..5350c97 100644
--- a/README.md
+++ b/README.md
@@ -176,6 +176,7 @@ The following reformulation methods are currently supported:
 
     - `optimizer`: Optimizer to use when solving subproblems to determine M values. This is a required value.
     - `default_M`: Default big-M value to use if no big-M is specified for a logical variable (1e9).
+    - `sampler`: M-value sampler for infinite models. Default: `ExhaustiveSampler()`, which solves an M subproblem at every support. Pass a `GPSampler` to instead solve a subset of the supports and fill the rest with a conservative Gaussian-process estimate. Ignored for finite models.
 
 5. [P-Split](https://arxiv.org/abs/2202.05198): This method reformulates each disjunct constraint into P constraints, each with a partitioned group defined by the user. This method requires that terms in the constraint be convex additively seperable with respect to each variable. The `PSplit` struct is created with the following required arguments:
 
@@ -223,6 +224,8 @@ optimize!(model, gdp_method = Hull())
 value(W)
 ```
 
+When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) enables an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. To opt in, pass a `GPSampler` via the `sampler` keyword, e.g. `MBM(optimizer, sampler = GPSampler())` (squared exponential kernel with the lengthscale selected by marginal likelihood), `GPSampler(Matern52Kernel())` (a custom kernel, lengthscale selected the same way), or `GPSampler(Matern52Kernel(), lengthscales = [0.2])` (lengthscale pinned). The filled values are heuristic upper estimates rather than certificates; see the `GPSampler` docstring for the tuning keywords (`kappa`, `budget`, `detect_uniform_M`, `lengthscales`, `jitter`, `seeds`).
+
 ## Release Notes
 
 Prior to `v0.4.0`, the package did not leverage the JuMP extension capabilities and was not as robust. For these earlier releases, refer to [Perez, Joshi, and Grossmann, 2023](https://arxiv.org/abs/2304.10492v1) and the following [JuliaCon 2022 Talk](https://www.youtube.com/watch?v=AMIrgTTfUkI).
diff --git a/ext/AbstractGPsDisjunctiveProgramming.jl b/ext/AbstractGPsDisjunctiveProgramming.jl
new file mode 100644
index 0000000..6efdf75
--- /dev/null
+++ b/ext/AbstractGPsDisjunctiveProgramming.jl
@@ -0,0 +1,173 @@
+module AbstractGPsDisjunctiveProgramming
+
+import AbstractGPs
+import AbstractGPs.KernelFunctions
+import DisjunctiveProgramming as DP
+
+################################################################################
+#                               SAMPLER CONFIG
+################################################################################
+# The concrete sampler behind DP.GPSampler; the base package only
+# carries the function stub, so constructing one requires AbstractGPs.
+struct _GPSampler{F} <: DP.AbstractMBMSampler
+    f::F
+    std_dev_margin::Float64
+    frac_supports::Float64
+    detect_uniform_M::Bool
+    initial_supports::Union{Int, Vector{Float64}}
+
+    function _GPSampler(
+        f::F;
+        std_dev_margin::Real = 2.5,
+        frac_supports::Real = 0.25,
+        detect_uniform_M::Bool = true,
+        initial_supports = 4
+        ) where {F}
+        f isa Union{Nothing, AbstractGPs.GP} || error(
+            "`f` must be an `AbstractGPs.GP` prior, e.g. " *
+            "`GP(Matern52Kernel())`.")
+        std_dev_margin >= 0 || error("`std_dev_margin` must be nonnegative.")
+        0 < frac_supports <= 1 || error("`frac_supports` must be in `(0, 1]`.")
+        if initial_supports isa Int
+            initial_supports >= 2 ||
+                error("`initial_supports` must be at least 2.")
+        else
+            initial_supports = collect(Float64, initial_supports)
+            (!isempty(initial_supports) &&
+                all(frac -> 0 <= frac <= 1, initial_supports)) ||
+                error("`initial_supports` must be fractions in `[0, 1]`.")
+        end
+        new{F}(f, Float64(std_dev_margin), Float64(frac_supports),
+            detect_uniform_M, initial_supports)
+    end
+end
+
+DP.GPSampler(f = nothing; kwargs...) = _GPSampler(f; kwargs...)
+
+################################################################################
+#                                 GP FITTING
+################################################################################
+# Normalized to [0, 1]^d so one lengthscale works across dimensions.
+# An independent parameter contributes one coordinate; a dependent
+# group contributes its joint-support column.
+function _support_coords(grids::Tuple)::Vector{Vector{Float64}}
+    axis_coords = map(grids) do g
+        g isa AbstractMatrix ? [g[:, j] for j in axes(g, 2)] :
+            [[v] for v in g]
+    end
+    coords = [reduce(vcat, getindex.(axis_coords, Tuple(I)))
+              for I in vec(CartesianIndices(length.(axis_coords)))]
+    dims = eachindex(first(coords))
+    mins = [minimum(c[d] for c in coords) for d in dims]
+    ranges = [max(maximum(c[d] for c in coords) - mins[d], eps())
+              for d in dims]
+    return [[(c[d] - mins[d]) / ranges[d] for d in dims]
+            for c in coords]
+end
+
+# Defaults behind GPSampler(): candidate lengthscales for the
+# marginal-likelihood fit (relative to the unit box) and the
+# observation-noise nugget.
+const _LENGTHSCALES = [0.05, 0.1, 0.2, 0.4, 0.8]
+const _JITTER = 1e-8
+
+# A user prior is used as given; the default squared exponential has
+# its lengthscale selected by marginal likelihood
+function _fit_posterior(
+    sampler::_GPSampler,
+    X::Vector{Vector{Float64}},
+    y::Vector{Float64}
+    )
+    sampler.f === nothing ||
+        return AbstractGPs.posterior(sampler.f(X, _JITTER), y)
+    best_posterior, best_log_prob = nothing, -Inf
+    for lengthscale in _LENGTHSCALES
+        kernel = KernelFunctions.with_lengthscale(
+            KernelFunctions.SqExponentialKernel(), lengthscale)
+        finite_gp = AbstractGPs.GP(kernel)(X, _JITTER)
+        log_prob = AbstractGPs.logpdf(finite_gp, y)
+        if log_prob > best_log_prob
+            best_posterior = AbstractGPs.posterior(finite_gp, y)
+            best_log_prob = log_prob
+        end
+    end
+    return best_posterior
+end
+
+function _mean_sd(
+    sampler::_GPSampler,
+    X::Vector{Vector{Float64}},
+    solved::Dict{Int, Float64}
+    )
+    solved_indices = collect(keys(solved))
+    y = [solved[i] for i in solved_indices]
+    y_mean = sum(y) / length(y)
+    # floored so near-equal solved values still cushion the filled ones
+    y_scale = max(sqrt(sum(abs2, y .- y_mean) / max(length(y) - 1, 1)),
+        1e-2 * abs(y_mean), 1e-8)
+    posterior = _fit_posterior(
+        sampler, X[solved_indices], (y .- y_mean) ./ y_scale)
+    posterior_mean = AbstractGPs.mean(posterior, X)
+    posterior_var = max.(AbstractGPs.var(posterior, X), 0.0)
+    return posterior_mean .* y_scale .+ y_mean,
+        sqrt.(posterior_var) .* y_scale
+end
+
+################################################################################
+#                              M VALUE SAMPLING
+################################################################################
+# Solve M at max-UCB selected supports, fill the rest with the bound
+function DP.sample_M_values(
+    sampler::_GPSampler,
+    objectives::AbstractArray,
+    sub::DP.GDPSubmodel,
+    method::DP._MBM,
+    support_grids::Tuple
+    )
+    indices = collect(CartesianIndices(objectives))
+    n = length(indices)
+    solved = Dict{Int, Float64}()
+    solve_at(index::Int) = begin
+        M_val = DP.raw_M(sub, objectives[indices[index]], method)
+        M_val === nothing && return false
+        solved[index] = M_val
+        return true
+    end
+    # an evenly spaced initial count, or user-given fractions
+    fractions = sampler.initial_supports isa Int ?
+        range(0, 1, length = sampler.initial_supports) :
+        sampler.initial_supports
+    for index in unique(1 .+ round.(Int, fractions .* (n - 1)))
+        solve_at(index) || return nothing
+    end
+    if sampler.detect_uniform_M
+        # a uniform M needs no fit
+        probes = collect(values(solved))
+        all(==(first(probes)), probes) && return first(probes)
+    end
+    solve_target = min(ceil(Int, sampler.frac_supports * n), n)
+    X = _support_coords(support_grids)
+    while length(solved) < solve_target
+        means, sds = _mean_sd(sampler, X, solved)
+        acquisition = means .+ sampler.std_dev_margin .* sds
+        for index in keys(solved)
+            acquisition[index] = -Inf
+        end
+        solve_at(argmax(acquisition)) || return nothing
+    end
+    M_vals = Array{Float64}(undef, size(objectives))
+    if length(solved) == n # nothing left to estimate
+        for (index, I) in enumerate(indices)
+            M_vals[I] = solved[index]
+        end
+        return M_vals
+    end
+    means, sds = _mean_sd(sampler, X, solved)
+    for (index, I) in enumerate(indices) # exact M values are nonnegative
+        M_vals[I] = get(solved, index,
+            max(means[index] + sampler.std_dev_margin * sds[index], 0.0))
+    end
+    return M_vals
+end
+
+end
diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl
index 5f77dce..2cc5c0a 100644
--- a/ext/InfiniteDisjunctiveProgramming.jl
+++ b/ext/InfiniteDisjunctiveProgramming.jl
@@ -238,56 +238,84 @@ function DP.prepare_max_M_objective(
     return obj.set.lower - obj_func
 end
 
-# Constant interpolation
-function _interpolate(
-    grids::NTuple{N, AbstractVector{<:Real}},
-    values::AbstractArray{<:Real, N}
-    ) where {N}
+# Candidate indices along one axis of the M array: the corners of the
+# grid cell bracketing a scalar query (independent parameter), or the
+# column matching a joint-support query (dependent group; every column
+# when the query is off-support, so the estimate stays conservative)
+function _axis_candidates(grid::AbstractVector{<:Real}, arg::Real)
+    lo = clamp(searchsortedlast(grid, arg), 1, length(grid) - 1)
+    return lo:(lo + 1)
+end
+function _axis_candidates(grid::AbstractMatrix{<:Real}, arg)
+    j = findfirst(k -> isapprox(view(grid, :, k), arg, atol = 1e-10),
+                  axes(grid, 2))
+    return isnothing(j) ? axes(grid, 2) : (j:j)
+end
+
+# Constant interpolation: max of `values` over the candidate indices
+function _interpolate(grids::Tuple, values::AbstractArray{<:Real})
     # mimic the call form of Interpolations.jl's interpolation
-    return (args...) -> _interpolate_at(grids, values, args)
+    return (args...) -> maximum(
+        values[I...] for I in Iterators.product(
+            map(_axis_candidates, grids, args)...))
+end
+
+# The infinite parameters of `mini_expr` and their supports, in the
+# ascending order of `parameter_refs`. An independent parameter gives
+# its sorted support vector; a dependent group gives the matrix whose
+# columns are its joint supports. Grids are read off the mini model so
+# their column order matches the transcription axes; the returned
+# prefs are the main-model parameters.
+function _support_grids(
+    sub::DP.GDPSubmodel, mini_expr::JuMP.AbstractJuMPScalar)
+    reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map)
+    mini_prefs = InfiniteOpt.parameter_refs(mini_expr)
+    prefs = Tuple(getindex.(Ref(reverse_map), p) for p in mini_prefs)
+    return prefs, Tuple(InfiniteOpt.supports(p) for p in mini_prefs)
 end
 
-function _interpolate_at(
-    grids::NTuple{N, AbstractVector{<:Real}},
-    values::AbstractArray{<:Real, N},
-    args::NTuple{N, <:Real}
-    ) where {N}
-    # lower-corner cell index per dimension
-    idx_lo = ntuple(d -> 
-        clamp(searchsortedlast(grids[d], args[d]),1, length(grids[d]) - 1), N
+# Solve the M subproblem exactly at every support
+function DP.sample_M_values(
+    sampler::DP.ExhaustiveSampler,
+    objectives::AbstractArray,
+    sub::DP.GDPSubmodel,
+    method::DP._MBM,
+    support_grids::Tuple
     )
-    # max over the 2^N corners; bit d of k picks lower or upper
-    return maximum(
-        values[ntuple(d -> idx_lo[d] +((k >> (d - 1)) & 1), N)...]
-        for k in 0:(2^N - 1)
-        )
+    M_vals = Array{Float64}(undef, size(objectives))
+    for I in eachindex(objectives)
+        m = DP.raw_M(sub, objectives[I], method)
+        m === nothing && return nothing
+        M_vals[I] = m
+    end
+    return M_vals
 end
 
-# Transcribe mini_expr, solve per support on the transcribed JuMP
-# model, and aggregate to a scalar if uniform, else to a parameter
-# function on main.
+# Transcribe mini_expr, compute the per-support M values with the
+# method's sampler, and aggregate to a scalar if uniform, else to a
+# parameter function on main.
 function DP.raw_M(
     sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel},
     mini_expr::JuMP.AbstractJuMPScalar,
     method::DP._MBM
     )
     objectives = InfiniteOpt.transformation_expression(mini_expr)
-    transcribed = InfiniteOpt.transformation_model(sub.model)
-    inner_sub = DP.GDPSubmodel(transcribed,JuMP.VariableRef[],
-        Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}()
-        )
-    M_vals = Array{typeof(method.default_M)}(undef, size(objectives))
-    for I in eachindex(objectives)
-        m = DP.raw_M(inner_sub, objectives[I], method)
-        m === nothing && return nothing
-        M_vals[I] = m
+    # transcription orders the dimensions by parameter group, which is
+    # not the ascending order `parameter_refs` gives the grids below
+    group_idxs = InfiniteOpt.parameter_group_int_indices(mini_expr)
+    if length(group_idxs) > 1 && ndims(objectives) == length(group_idxs)
+        objectives = permutedims(objectives, sortperm(group_idxs))
     end
+    transcribed = InfiniteOpt.transformation_model(sub.model)
+    inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[],
+        Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}())
+    prefs, grids = _support_grids(sub, mini_expr)
+    M_vals = DP.sample_M_values(method.sampler, objectives,
+        inner_sub, method, grids)
+    M_vals === nothing && return nothing
+    M_vals isa Number && return M_vals
     all(==(first(M_vals)), M_vals) && return first(M_vals)
-    mini_prefs = InfiniteOpt.parameter_refs(mini_expr)
-    reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map)
-    prefs = Tuple(reverse_map[p] for p in mini_prefs)
-    main = JuMP.owner_model(first(prefs))
-    grids = Tuple(InfiniteOpt.supports(p) for p in prefs)
+    main = JuMP.owner_model(first(keys(sub.fwd_map)))
     param_func = InfiniteOpt.build_parameter_function(
         error, _interpolate(grids, M_vals), prefs)
     return InfiniteOpt.add_parameter_function(main, param_func)
diff --git a/src/datatypes.jl b/src/datatypes.jl
index bdcd4c2..434cd1d 100644
--- a/src/datatypes.jl
+++ b/src/datatypes.jl
@@ -368,38 +368,72 @@ struct BigM{T} <: AbstractReformulationMethod
 end
 
 """
-    MBM{O, T, L <: LogicalVariableRef} <: AbstractReformulationMethod
+    AbstractMBMSampler
+
+Abstract type for the M-value samplers used by [`MBM`](@ref) on
+infinite models. Concrete samplers implement
+[`sample_M_values`](@ref).
+"""
+abstract type AbstractMBMSampler end
+
+"""
+    ExhaustiveSampler <: AbstractMBMSampler
+
+The default M-value sampler for [`MBM`](@ref): solve an M subproblem
+at every support of the infinite model.
+"""
+struct ExhaustiveSampler <: AbstractMBMSampler end
+
+"""
+    MBM{O, T, S} <: AbstractReformulationMethod
 
 A type for using the multiple big-M reformulation approach for disjunctive constraints.
 
 **Fields**
 - `optimizer::O`: Optimizer to use when solving mini-models (required).
 - `default_M::T`: Default big-M value to use if no big-M is specified for a logical variable (1e9).
-"""
-mutable struct MBM{O, T} <: AbstractReformulationMethod
+- `sampler::S`: M-value sampler for infinite models
+  (`ExhaustiveSampler()`). An [`ExhaustiveSampler`](@ref) solves an M
+  subproblem at every support; a [`GPSampler`](@ref) solves a subset
+  of the supports and fills the rest with a conservative
+  Gaussian-process estimate (see [`sample_M_values`](@ref)). Ignored
+  for finite models.
+"""
+mutable struct MBM{O, T, S <: AbstractMBMSampler} <:
+    AbstractReformulationMethod
     optimizer::O
     default_M::T
-    
+    sampler::S
+
     # Constructor with optimizer (required) and optional default_M
-    function MBM(optimizer::O, default_M::T = 1e9) where {O, T}
-        new{O, T}(optimizer, default_M)
+    # (kwargs cannot bind static parameters, hence typeof(sampler))
+    function MBM(
+        optimizer::O, default_M::T = 1e9;
+        sampler = ExhaustiveSampler()
+        ) where {O, T}
+        new{O, T, typeof(sampler)}(optimizer, default_M, sampler)
     end
 end
 
-mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMethod
+mutable struct _MBM{O, T, S <: AbstractMBMSampler,
+    M <: JuMP.AbstractModel} <: AbstractReformulationMethod
     optimizer::O
     M::Dict{LogicalVariableRef{M}, Any}
     default_M::T
+    sampler::S
     subproblem_indicators::Vector{LogicalVariableRef{M}}
     # Cached submodels: indicator => GDPSubmodel.
     # Typed Any so extensions can store different types.
     model_cache::Dict{LogicalVariableRef{M}, Any}
 
-    function _MBM(method::MBM{O, T}, model::M) where {O, T, M <: JuMP.AbstractModel}
-        new{O, T, M}(
+    function _MBM(
+        method::MBM{O, T, S}, model::M
+        ) where {O, T, S, M <: JuMP.AbstractModel}
+        new{O, T, S, M}(
             method.optimizer,
             Dict{LogicalVariableRef{M}, Any}(),
             method.default_M,
+            method.sampler,
             Vector{LogicalVariableRef{M}}(),
             Dict{LogicalVariableRef{M}, Any}()
         )
diff --git a/src/extension_api.jl b/src/extension_api.jl
index ad3566f..02efd8f 100644
--- a/src/extension_api.jl
+++ b/src/extension_api.jl
@@ -38,3 +38,82 @@ Y(t, x)
 ```
 """
 function InfiniteLogical end
+
+"""
+    GPSampler(
+        f = nothing;
+        std_dev_margin::Real = 2.5,
+        frac_supports::Real = 0.25,
+        detect_uniform_M::Bool = true,
+        initial_supports = 4
+        )
+
+A Gaussian-process M sampler for the `sampler` field of [`MBM`](@ref)
+on infinite models. Instead of solving an M subproblem at every
+support, it solves the initial supports, then the supports selected
+by an upper-confidence-bound acquisition until `frac_supports` of
+them are solved, and fills the remaining supports with the posterior
+upper confidence bound `mean + std_dev_margin * sd`. The filled
+values are heuristic upper estimates of the exact M values, not
+certificates. This requires that AbstractGPs be imported first.
+
+**Arguments**
+- `f`: An `AbstractGPs.GP` prior, used as given, e.g.
+  `GP(Matern52Kernel())`; `nothing` (the default) uses a squared
+  exponential kernel with its lengthscale selected by marginal
+  likelihood. The prior is fit to support coordinates normalized to
+  `[0, 1]` per dimension and to M values standardized to zero mean
+  and unit scale, so any lengthscale baked into `f` is relative to
+  the unit box.
+- `std_dev_margin::Real`: Standard deviations added above the
+  posterior mean, both to select the next support to solve and to
+  fill the unsolved supports (2.5).
+- `frac_supports::Real`: Fraction of the supports to solve exactly,
+  in `(0, 1]` (0.25). The initial supports are always solved;
+  `frac_supports = 1.0` solves every support.
+- `detect_uniform_M::Bool`: If `true` (the default), M values that
+  agree at the initial supports are taken to be uniform and used for
+  every support. Set it to `false` to always fit the GP, which
+  leaves the usual `std_dev_margin * sd` cushion on the unsolved
+  supports at the cost of the extra solves.
+- `initial_supports`: Number of evenly spaced supports solved before
+  the first GP fit (4), or a vector of fractions in `[0, 1]` giving
+  their positions along the supports. With `detect_uniform_M`,
+  evenly spaced initial supports can read a periodic M as uniform;
+  pass unevenly spaced fractions to guard against that.
+
+**Example**
+```julia
+julia> using DisjunctiveProgramming, InfiniteOpt, AbstractGPs, HiGHS
+
+julia> method = MBM(HiGHS.Optimizer,
+           sampler = GPSampler(GP(Matern52Kernel()), std_dev_margin = 4.0))
+```
+"""
+function GPSampler end
+
+"""
+    sample_M_values(sampler, objectives, sub, method, support_grids)
+
+Compute the MBM M values at the transcription supports of an infinite
+model. `sampler` is the `sampler` field of [`MBM`](@ref),
+`objectives` is the array of per-support objective expressions, `sub`
+is the transcribed submodel wrapped as a `GDPSubmodel`, `method` is
+the `_MBM` data, and `support_grids` is a tuple of the supports of
+the infinite parameters (a sorted vector per independent parameter,
+a matrix of joint supports per dependent group).
+Returns an array of M values shaped like `objectives`, a
+scalar when M is uniform across the supports, or `nothing` if an M
+subproblem is infeasible. Extensions implement methods that dispatch
+on `sampler`: an [`ExhaustiveSampler`](@ref) solves an M subproblem
+at every support, while a [`GPSampler`](@ref) solves a subset of the
+supports and
+fills the rest with a Gaussian-process upper confidence bound.
+"""
+function sample_M_values(sampler, objectives, sub, method, support_grids)
+    error("Unrecognized `sampler` value `$(repr(sampler))` for MBM " *
+          "on an infinite model. Use `ExhaustiveSampler()` to solve " *
+          "an M subproblem at every support, or a `GPSampler` (with " *
+          "AbstractGPs loaded) to estimate M values with a " *
+          "Gaussian process.")
+end
diff --git a/src/mbm.jl b/src/mbm.jl
index 23c5096..8f74269 100644
--- a/src/mbm.jl
+++ b/src/mbm.jl
@@ -87,7 +87,9 @@ function reformulate_disjunct_constraint(
         },
     method::_MBM
     )
-    ref_cons = reformulate_disjunction(model, con, MBM(method.optimizer))
+    ref_cons = reformulate_disjunction(model, con, MBM(
+        method.optimizer, method.default_M,
+        sampler = method.sampler))
     new_ref_cons = Vector{JuMP.AbstractConstraint}()
     for ref_con in ref_cons
         append!(new_ref_cons,
diff --git a/test/extensions/AbstractGPsDisjunctiveProgramming.jl b/test/extensions/AbstractGPsDisjunctiveProgramming.jl
new file mode 100644
index 0000000..4fd11cb
--- /dev/null
+++ b/test/extensions/AbstractGPsDisjunctiveProgramming.jl
@@ -0,0 +1,310 @@
+using InfiniteOpt, HiGHS, AbstractGPs
+import DisjunctiveProgramming as DP
+
+# subtype without a sample_M_values method, for the fallback error
+struct _UnimplementedSampler <: DP.AbstractMBMSampler end
+
+function test_gp_sampler_kwargs()
+    @test MBM(HiGHS.Optimizer).sampler === ExhaustiveSampler()
+    sampler = GPSampler()
+    @test sampler.f === nothing
+    @test sampler.std_dev_margin == 2.5
+    @test sampler.frac_supports == 0.25
+    @test sampler.detect_uniform_M
+    @test sampler.initial_supports == 4
+    f = GP(Matern52Kernel())
+    sampler = GPSampler(f, std_dev_margin = 4.0, frac_supports = 0.1,
+        detect_uniform_M = false, initial_supports = [0.0, 0.3, 1.0])
+    @test sampler.f === f
+    @test sampler.std_dev_margin == 4.0
+    @test sampler.frac_supports == 0.1
+    @test !sampler.detect_uniform_M
+    @test sampler.initial_supports == [0.0, 0.3, 1.0]
+    @test GPSampler(initial_supports = 6).initial_supports == 6
+    @test MBM(HiGHS.Optimizer, sampler = sampler).sampler === sampler
+    # a bare kernel is not a prior
+    @test_throws ErrorException GPSampler(SqExponentialKernel())
+    @test_throws ErrorException GPSampler(std_dev_margin = -1)
+    @test_throws ErrorException GPSampler(frac_supports = 0)
+    @test_throws ErrorException GPSampler(frac_supports = 1.5)
+    @test_throws ErrorException GPSampler(initial_supports = 1)
+    @test_throws ErrorException GPSampler(initial_supports = [1.5])
+    @test_throws ErrorException GPSampler(initial_supports = Float64[])
+end
+
+# Mirror of test_raw_M_infinite_scalar: uniform seed M values collapse
+# to the exactly-solved scalar under the GP sampler
+function test_gp_raw_M_scalar()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @constraint(model, con, x >= 5, Disjunct(Y[1]))
+    @constraint(model, con2, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(
+        MBM(HiGHS.Optimizer, sampler = GPSampler()), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    @test DP.raw_M(sub, obj, mbm) == 5.0
+end
+
+# With frac_supports = 1.0 every support is solved exactly, so the GP
+# sampler must reproduce the exact grid parameter function
+function test_gp_raw_M_matches_exact()
+    function pfunc_values(sampler, supports)
+        model = InfiniteGDPModel()
+        @infinite_parameter(model, t ∈ [0, 1], supports = supports)
+        @variable(model, 0 <= x <= 10, Infinite(t))
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @parameter_function(model, f == t -> 2*t)
+        @constraint(model, con, x <= f, Disjunct(Y[1]))
+        @constraint(model, con2, x >= 0.5, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        mbm = DP._MBM(
+            MBM(HiGHS.Optimizer, sampler = sampler), model)
+        sub = DP.copy_model_with_constraints(
+            model, DP.DisjunctConstraintRef[con2], mbm)
+        obj = DP.prepare_max_M_objective(
+            model, JuMP.constraint_object(con), sub)
+        M = DP.raw_M(sub, obj, mbm)
+        @test M isa InfiniteOpt.GeneralVariableRef
+        return [InfiniteOpt.raw_function(M)(t_val) for t_val in supports]
+    end
+    supports = [0.0, 0.25, 0.5, 0.75, 1.0]
+    exact_vals = pfunc_values(ExhaustiveSampler(), supports)
+    # the default prior, a user prior, and a pinned lengthscale all
+    # solve the same supports here, so the M values match exactly
+    @test pfunc_values(GPSampler(frac_supports = 1.0), supports) ==
+        exact_vals
+    @test pfunc_values(
+        GPSampler(GP(Matern52Kernel()), frac_supports = 1.0), supports) ==
+        exact_vals
+    @test pfunc_values(GPSampler(GP(with_lengthscale(
+        SqExponentialKernel(), 0.2)), frac_supports = 1.0), supports) ==
+        exact_vals
+end
+
+# Two independent parameters: the GP path builds 2-D coordinates and
+# fits a multivariate GP; with frac_supports = 1.0 every support is solved
+# exactly, so the parameter function matches the exhaustive one. Setup
+# as in test_raw_M_infinite_two_params: M(t, s) = 10 - t - s.
+function test_gp_raw_M_two_params()
+    function pfunc_values(sampler)
+        model = InfiniteGDPModel()
+        @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+        @infinite_parameter(model, s ∈ [0, 1], supports = [0.0, 1.0])
+        @variable(model, 0 <= x <= 10, Infinite(t, s))
+        @variable(model, Y[1:2], InfiniteLogical(t, s))
+        @constraint(model, con, x <= t + s, Disjunct(Y[1]))
+        @constraint(model, con2, x >= 0.5, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        mbm = DP._MBM(
+            MBM(HiGHS.Optimizer, sampler = sampler), model)
+        sub = DP.copy_model_with_constraints(
+            model, DP.DisjunctConstraintRef[con2], mbm)
+        obj = DP.prepare_max_M_objective(
+            model, JuMP.constraint_object(con), sub)
+        M = DP.raw_M(sub, obj, mbm)
+        @test M isa InfiniteOpt.GeneralVariableRef
+        raw_fn = InfiniteOpt.raw_function(M)
+        return [raw_fn(t_val, s_val)
+                for t_val in [0.0, 0.5, 1.0], s_val in [0.0, 1.0]]
+    end
+    @test pfunc_values(GPSampler(frac_supports = 1.0)) ==
+        pfunc_values(ExhaustiveSampler())
+end
+
+# Dependent parameters: the joint supports become the GP coordinates
+# directly; with frac_supports = 1.0 every support is solved exactly, so the
+# M values match the exhaustive ones. Setup as in
+# test_raw_M_infinite_dependent_varying: M(ξ) = 10 - ξ[1] - ξ[2].
+function test_gp_raw_M_dependent()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 6)
+    @variable(model, 0 <= x <= 10, Infinite(ξ))
+    @variable(model, Y[1:2], InfiniteLogical(ξ))
+    @constraint(model, con, x <= ξ[1] + ξ[2], Disjunct(Y[1]))
+    @constraint(model, con2, x >= 0.5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(
+        MBM(HiGHS.Optimizer, sampler = GPSampler(frac_supports = 1.0)), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    M = DP.raw_M(sub, obj, mbm)
+    @test M isa InfiniteOpt.GeneralVariableRef
+    raw_fn = InfiniteOpt.raw_function(M)
+    S = InfiniteOpt.supports(ξ)
+    for j in axes(S, 2)
+        @test raw_fn(S[:, j]) ≈ 10.0 - S[1, j] - S[2, j] atol = 1e-6
+    end
+end
+
+# an empty disjunct region makes the M subproblems infeasible; both
+# samplers propagate that up to the reformulation error
+function test_gp_infeasible_disjunct()
+    function build()
+        model = InfiniteGDPModel(HiGHS.Optimizer)
+        set_silent(model)
+        @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+        @variable(model, 0 <= x <= 10, Infinite(t))
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @parameter_function(model, f == t -> 2*t)
+        @constraint(model, x <= f, Disjunct(Y[1]))
+        @constraint(model, x >= 8, Disjunct(Y[2]))
+        @constraint(model, x <= 3, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        @objective(model, Max, 𝔼(x, t))
+        return model
+    end
+    for sampler in (ExhaustiveSampler(), GPSampler())
+        model = build()
+        @test_throws ErrorException optimize!(model,
+            gdp_method = MBM(HiGHS.Optimizer, sampler = sampler))
+    end
+end
+
+# optimum (10) needs M(t) >= 10 - 2t pointwise; the GP fill is heuristic
+function test_gp_mbm_solve_equivalence()
+    function solve_with(sampler)
+        model = InfiniteGDPModel(HiGHS.Optimizer)
+        set_silent(model)
+        @infinite_parameter(model, t ∈ [0, 1], num_supports = 20)
+        @variable(model, 0 <= x <= 10, Infinite(t))
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @parameter_function(model, f == t -> 2*t)
+        @constraint(model, x <= f, Disjunct(Y[1]))
+        @constraint(model, x >= 0.5, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        @objective(model, Max, 𝔼(x, t))
+        optimize!(model,
+            gdp_method = MBM(HiGHS.Optimizer, sampler = sampler))
+        @test termination_status(model) == MOI.OPTIMAL
+        return objective_value(model)
+    end
+    obj_exact = solve_with(ExhaustiveSampler())
+    obj_gp = solve_with(GPSampler())
+    obj_tuned = solve_with(GPSampler(std_dev_margin = 4.0, frac_supports = 0.2))
+    @test obj_exact ≈ 10.0 atol = 1e-4
+    # over-M can't raise the optimum, under-M can only shave it a bit
+    @test obj_gp <= obj_exact + 1e-6
+    @test obj_gp ≈ obj_exact atol = 1e-2
+    @test obj_tuned <= obj_exact + 1e-6
+    @test obj_tuned ≈ obj_exact atol = 1e-2
+end
+
+# Seed placement vs a periodic M. With f(t) = 2|cos(2*pi*t)| on these
+# supports, M = 10 - f is 8 at supports 1, 3, 5 and 10 at supports
+# 2, 4. Seeds that only hit the M = 8 supports alias the periodic M
+# to uniform 8, which caps x at 8 and cuts the optimum from 10 down
+# to 9; the default and denser seed grids see both values and stay
+# exact.
+function test_gp_periodic_M_seeds()
+    supports = [0.0, 0.25, 0.5, 0.75, 1.0]
+    function solve_with(sampler)
+        model = InfiniteGDPModel(HiGHS.Optimizer)
+        set_silent(model)
+        @infinite_parameter(model, t ∈ [0, 1], supports = supports)
+        @variable(model, 0 <= x <= 10, Infinite(t))
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @parameter_function(model, f == t -> 2 * abs(cos(2 * pi * t)))
+        @constraint(model, x <= f, Disjunct(Y[1]))
+        @constraint(model, x >= 0.5, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        @objective(model, Max, 𝔼(x, t))
+        optimize!(model,
+            gdp_method = MBM(HiGHS.Optimizer, sampler = sampler))
+        return objective_value(model)
+    end
+    @test solve_with(ExhaustiveSampler()) ≈ 10.0 atol = 1e-6
+    @test solve_with(GPSampler()) ≈ 10.0 atol = 1e-6
+    @test solve_with(GPSampler(initial_supports = 5)) ≈ 10.0 atol = 1e-6
+    @test solve_with(GPSampler(initial_supports = [0.0, 0.5, 1.0])) ≈
+        9.0 atol = 1e-6
+end
+
+# With detection off the uniform M is not collapsed to a scalar: the
+# GP is fit and the unsolved supports keep their std_dev_margin * sd cushion,
+# which must sit above the M that detection would have returned.
+function test_gp_detect_uniform_M_off()
+    function raw_M_with(detect)
+        model = InfiniteGDPModel()
+        @infinite_parameter(model, t ∈ [0, 1], num_supports = 20)
+        @variable(model, 0 <= x <= 10, Infinite(t))
+        @variable(model, Y[1:2], InfiniteLogical(t))
+        @constraint(model, con, x >= 5, Disjunct(Y[1]))
+        @constraint(model, con2, x <= 3, Disjunct(Y[2]))
+        @disjunction(model, Y)
+        mbm = DP._MBM(MBM(HiGHS.Optimizer,
+            sampler = GPSampler(detect_uniform_M = detect)), model)
+        sub = DP.copy_model_with_constraints(
+            model, DP.DisjunctConstraintRef[con2], mbm)
+        obj = DP.prepare_max_M_objective(
+            model, JuMP.constraint_object(con), sub)
+        return DP.raw_M(sub, obj, mbm)
+    end
+    @test raw_M_with(true) == 5.0
+    M = raw_M_with(false)
+    @test M isa InfiniteOpt.GeneralVariableRef
+    raw_fn = InfiniteOpt.raw_function(M)
+    vals = [raw_fn(t) for t in range(0, 1, length = 20)]
+    @test all(vals .>= 5.0 - 1e-6)
+    @test maximum(vals) > 5.0
+end
+
+# Dependent parameters have no support grid, so turning detection off
+# leaves the GP with nothing to fit over
+# detect_uniform_M = false forces coordinate construction from the
+# joint supports; the uniform M still comes back exact
+function test_gp_detect_uniform_M_off_dependent()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4)
+    @variable(model, 0 <= x <= 10, Infinite(ξ))
+    @variable(model, Y[1:2], InfiniteLogical(ξ))
+    @constraint(model, con, x >= 5, Disjunct(Y[1]))
+    @constraint(model, con2, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer,
+        sampler = GPSampler(detect_uniform_M = false)), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    @test DP.raw_M(sub, obj, mbm) == 5.0
+end
+
+function test_gp_unknown_sampler_error()
+    # a non-AbstractMBMSampler is rejected at construction
+    @test_throws TypeError MBM(HiGHS.Optimizer, sampler = :grid)
+    model = InfiniteGDPModel(HiGHS.Optimizer)
+    set_silent(model)
+    @infinite_parameter(model, t ∈ [0, 1], num_supports = 5)
+    @variable(model, 0 <= x <= 10, Infinite(t))
+    @variable(model, Y[1:2], InfiniteLogical(t))
+    @parameter_function(model, f == t -> 2*t)
+    @constraint(model, x <= f, Disjunct(Y[1]))
+    @constraint(model, x >= 0.5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    @objective(model, Max, 𝔼(x, t))
+    @test_throws ErrorException optimize!(model,
+        gdp_method = MBM(HiGHS.Optimizer,
+            sampler = _UnimplementedSampler()))
+end
+
+@testset "AbstractGPsDisjunctiveProgramming" begin
+    test_gp_sampler_kwargs()
+    test_gp_raw_M_scalar()
+    test_gp_raw_M_matches_exact()
+    test_gp_raw_M_two_params()
+    test_gp_raw_M_dependent()
+    test_gp_mbm_solve_equivalence()
+    test_gp_periodic_M_seeds()
+    test_gp_detect_uniform_M_off()
+    test_gp_detect_uniform_M_off_dependent()
+    test_gp_unknown_sampler_error()
+    test_gp_infeasible_disjunct()
+end
diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl
index 47bb884..3f05a4d 100644
--- a/test/extensions/InfiniteDisjunctiveProgramming.jl
+++ b/test/extensions/InfiniteDisjunctiveProgramming.jl
@@ -413,6 +413,82 @@ function test_raw_M_infinite_param_function()
     end
 end
 
+# raw_M over two infinite parameters with different support counts.
+# Transcription orders the objective dimensions by parameter group,
+# which need not be the ascending order of the grids, so the M values
+# must be permuted to line up. Setup: x(t, s) in [0, 10],
+# disj1: x <= t + s, disj2: x >= 0.5. Slack r(x) = x - t - s
+# maximized over x in [0.5, 10]: 10 - t - s.
+function test_raw_M_infinite_two_params()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0])
+    @infinite_parameter(model, s ∈ [0, 1], supports = [0.0, 1.0])
+    @variable(model, 0 <= x <= 10, Infinite(t, s))
+    @variable(model, Y[1:2], InfiniteLogical(t, s))
+    @constraint(model, con, x <= t + s, Disjunct(Y[1]))
+    @constraint(model, con2, x >= 0.5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    M = DP.raw_M(sub, obj, mbm)
+    @test M isa InfiniteOpt.GeneralVariableRef
+    raw_fn = InfiniteOpt.raw_function(M)
+    for t_val in [0.0, 0.5, 1.0], s_val in [0.0, 1.0]
+        @test raw_fn(t_val, s_val) >= 10.0 - t_val - s_val - 1e-6
+    end
+end
+
+# Dependent parameters with a uniform M short-circuit to a scalar
+# before any support grid is needed. Setup as in
+# test_raw_M_infinite_scalar, over a dependent parameter array.
+function test_raw_M_infinite_dependent_params()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4)
+    @variable(model, 0 <= x <= 10, Infinite(ξ))
+    @variable(model, Y[1:2], InfiniteLogical(ξ))
+    @constraint(model, con, x >= 5, Disjunct(Y[1]))
+    @constraint(model, con2, x <= 3, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    @test DP.raw_M(sub, obj, mbm) == 5.0
+end
+
+# Dependent parameters with M varying over the joint supports: the
+# parameter function looks M up at each joint support and falls back
+# to the conservative max off-support. Setup: x(ξ) in [0, 10],
+# disj1: x <= ξ[1] + ξ[2], disj2: x >= 0.5, so M(ξ) = 10 - ξ1 - ξ2.
+function test_raw_M_infinite_dependent_varying()
+    model = InfiniteGDPModel()
+    @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4)
+    @variable(model, 0 <= x <= 10, Infinite(ξ))
+    @variable(model, Y[1:2], InfiniteLogical(ξ))
+    @constraint(model, con, x <= ξ[1] + ξ[2], Disjunct(Y[1]))
+    @constraint(model, con2, x >= 0.5, Disjunct(Y[2]))
+    @disjunction(model, Y)
+    mbm = DP._MBM(MBM(HiGHS.Optimizer), model)
+    sub = DP.copy_model_with_constraints(
+        model, DP.DisjunctConstraintRef[con2], mbm)
+    obj = DP.prepare_max_M_objective(
+        model, JuMP.constraint_object(con), sub)
+    M = DP.raw_M(sub, obj, mbm)
+    @test M isa InfiniteOpt.GeneralVariableRef
+    raw_fn = InfiniteOpt.raw_function(M)
+    S = InfiniteOpt.supports(ξ)
+    expected = [10.0 - S[1, j] - S[2, j] for j in axes(S, 2)]
+    for j in axes(S, 2)
+        @test raw_fn(S[:, j]) ≈ expected[j] atol = 1e-6
+    end
+    # off-support queries fall back to the maximum over all supports
+    @test raw_fn([0.1234, 0.5678]) ≈ maximum(expected) atol = 1e-6
+end
+
 # Piecewise-constant max-of-corners: returns the maximum value over
 # the 2^n corners of the cell containing the query.
 function test_interpolate()
@@ -834,6 +910,9 @@ end
         test_interpolate()
         test_raw_M_infinite_scalar()
         test_raw_M_infinite_param_function()
+        test_raw_M_infinite_two_params()
+        test_raw_M_infinite_dependent_params()
+        test_raw_M_infinite_dependent_varying()
         test_mbm_finite_and_integer_var()
         test_mbm_infinite_simple()
         test_mbm_infinite_param_dependent()
diff --git a/test/runtests.jl b/test/runtests.jl
index 06e8813..8034f6e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,4 +24,5 @@ include("constraints/disjunction.jl")
 include("print.jl")
 include("solve.jl")
 include("extensions/InfiniteDisjunctiveProgramming.jl")
+include("extensions/AbstractGPsDisjunctiveProgramming.jl")