DIRECT 算法：算法原理与 Julia 代码实现

本文主要信息来源于 1993 年提出 DIRECT 算法的原文 ¹ .

[!tldr] 本文导读
DIRECT 算法是一种全局的优化算法，并且可以优化无噪声的黑盒函数.

1. 它在维度较高时的性质也相比以前的优化算法（Shubert 算法等）更稳定；
2. 找到最优值的必要条件：函数连续，或者至少在全局最优点的邻域附近连续.
3. 适用范围：搜索空间为闭方体、Lipschitz 条件函数.

一维的全局优化：从 Shubert 算法开始

本文的开头先介绍了一种一维的 Lipschitz 优化算法，它能够很好地解决一维 Lipschitz 连续函数的全局优化 ² . 这个算法称为 Shubert 算法. 它的思路很简单，就是利用了 Lipschitz 函数一个很重要的性质，可以用 V 形的一个分段线性函数作为下界：

设 Lipschitz 函数 $f$ 满足

$$ |f(x) - f(x')| \leqslant K |x -x'|, \quad \forall x,x'\in [a,b] $$

这样的话就可以知道

$$ \begin{aligned} f(x) & \geqslant f(a) - K(x-a) \\ f(x) & \geqslant f(b) + K(x-b) \end{aligned} $$

其实这就是上面图片的 V 形，图中对应的量为：

$$ \begin{aligned} X(a,b,f,K) & = \frac{a+b}{2} + \frac{f(a)-f(b)}{2K} \\ B(a,b,f,K) & = \frac{f(a)+f(b)}{2} - K(b-a) \end{aligned} $$

Shubert 算法的想法是：如果我们不断改进这个下界，那么根据 Lipschitz 函数的性质，总可以得到真正的全局最优解：

改进的方法很简单，就是不断寻找 V 形的折点，然后通过获得折点上的函数值，更新 Lipschitz 下界，不断操作就可以逼近真正的下界.

需要注意的是：

• Shubert 算法要求对函数的 Lipschitz 常数 $K$ 有一定的认知，如果 $K$ 选取过小，则算法收敛不到最优值（这个是显然的，因为这将导致 V 形不再是真正的下界），但是选择过大将导致收敛较慢.
• Shubert 算法可以扩展到高维，也就是用超平面，但是 $n$ 维就需要计算 $2^{n}$ 个边界点的函数值，这在极端高维情况下是几乎不可能实现的.

算法代码见文末，可让 AI 解释其中的细节.

DIRECT 算法就是在这样的背景下诞生的，下面我们先引入一维的 DIRECT 算法.

一维的 DIRECT 算法

Shubert 算法的需要改进点在于高维情况下的边界点过多，计算成本过大，因此 DIRECT 算法在这个上面进行了修改：考虑的不再是边界点，而是中点，一维情形下是区间中点，高维情形下是各坐标中点，也就是超矩形的中心.

我们用一维的图像进行理解，实际上就是先计算 $c = \dfrac{b-a}{2}$ 处的函数值 $f(c)$ ，此时根据 Lipschitz 条件有

$$ \begin{aligned} & f(x) \geqslant f(c) + K(x-c) ,\quad \forall x \leqslant c\\ & f(x) \geqslant f(c) - K(x-c) ,\quad \forall x \geqslant c\\ \end{aligned} $$

在图上就是反 V 形，下界就是端点处，可以知道是

$$ \text{lower bound} = f(c) - \frac{(b-a)K}{2} $$

我们令 $d$ 表示区间的半长，那么 $f(c)-Kd$ 其实就代表了一个区间的“潜在最小”的可能性，我们如果将这个作为下界去不断逼近真正的下界，是不是就能真正达到最优解呢？基于这种想法，作者给出了“潜在最优区间”的概念.

[!note] 定义：潜在最优区间
假定现在已经将区间 $[a,b]$ 分割为子区间 $[a_{i}, b_{i}]$ ，中点分别为 $c_{i}$ ，其中 $i=1,2,\cdots, m$ ，令 $\varepsilon>0$ 为常数，$f_{\min}$ 是已经计算的目前最优函数值. 如果对于区间 $j$ ，存在常数 $\widetilde{K}>0$ ，使得

$$ f(c_{j}) - \widetilde{K} d_{j} \leqslant f(c_{i}) - \widetilde{K} d_{i}, \quad \forall i=1,2,\cdots, m $$

且

$$ f(c_{j}) - \widetilde{K}d_{j} \leqslant f_{\min} - \varepsilon|f_{\min}| $$

则称 $j$ 是潜在最优区间.

这里是 DIRECT 算法最精华的部分之一，因为潜在最优此时无需预先给定 Lipschitz 常数，这就使得其有较高的适应性，只要保证 Lipschitz 常数存在即可.

第二个不等式的来源，原文当中是这样说的：

The second condition insists that the lower bound for the interval, based on the rate-of-change constant $\widetilde{K}$, exceed the current best solution by a nontrivial amount. For example, if $\varepsilon=0.01$ , then the lower bound for the interval would have to exceed the current best solution by more than 1%. This second condition is needed to prevent the algorithm from becoming too local in its orientation, wasting precious function evaluations in pursuit of extremely small improvements

也就是避免算法过于局部化，造成浪费计算资源，如果不能超过目前最优的 $\varepsilon$ 比例，就不考虑继续在这个局部迭代下去. 另外，如果 $K$ 事实上已知，则算法在实际实现时需要令 $\widetilde{K} \leqslant K$ .

为了找寻潜在最优区间，其实对于一个区间只需要记录 $(d, f(c))$ ，对于相同的 $d$ ，中点函数值 $f(c)$ 较小的区间才是潜在最优区间，否则就不是，因此当我们有很多子区间的时候，就可以画出如下的图像：

这里面的每个点都是一个区间对应的 $(d,f(c))$ ，当用斜率 $K$ 的直线从下往上移动，最先触碰到的点就是 $\widetilde{K}$ 对应的潜在最优区间. 那么，实质上这个过程就是找右下凸包的过程：

在 DIRECT 算法的实际操作当中，不断进行三等分获得子区间，然后选择潜在最优进行进一步三等分. 算法伪代码见下图.

为加深理解，这里我们给出一个一维的算法实例.

[!faq] 例：一维函数全局最优化
对 $x^{2}$ 寻找 $[-1, 2]$ 上的最优解，适用 DIRECT 算法迭代两轮.

区间 $[-1,2]$ ，那么中心点为 $c_{1}=\dfrac{1}{2}$ ，半长为 $d_{1}=\dfrac{3}{2}$ ，此时的当前最优为 $f(c_{1}) = \dfrac{1}{4}$ ，令 $\varepsilon$ 初始化为 $10^{-4}$ .

由于初始仅有一个区间 $[-1,2]$ ，因此 $S = \left\lbrace [-1,2] \right\rbrace$ . 进入第一次迭代：$\delta = 1$ ，即区间的三分之一长，三等分为：

• 左区间：$[-1,0]$ ，中心为 $- \dfrac{1}{2}$ ，$f=0.25$ .
• 中区间：$[0,1]$ ，中心为 $\dfrac{1}{2}$ ，$f=0.25$ .
• 右区间：$[1,2]$ ，中心为 $\dfrac{3}{2}$ ，$f=2.25$ .

将 $[-1,2]$ 从 $S$ 当中移除，此时 $S$ 变空，第一轮结束，第二轮开始，考虑计算潜在最优区间，对于三个区间：

$$ \begin{cases} d_{1} = d_{2}=d_{3} = 0.5 \\ f(c_{1}) = f(c_{2}) = 0.25 \\ f(c_{3})=2.25 \end{cases} $$

因此左区间和中区间都是潜在最优，因此 $S=\left\lbrace [-1,0], [0,1] \right\rbrace$ ，下一步就是继续分割：

• 对 $[-1,0]$ ：

  • 左区间：$\left[-1, - \dfrac{2}{3}\right]$ ，中心为 $- \dfrac{5}{6}$ ，$f= \dfrac{25}{36}$ .
  • 中区间：$\left[- \dfrac{2}{3}, - \dfrac{1}{3}\right]$ ，中心为 $-\dfrac{1}{2}$ ，$f=0.25$ .
  • 右区间：$\left[- \dfrac{1}{3},0\right]$ ，中心为 $-\dfrac{1}{6}$ ，$f= \dfrac{1}{36}$ .
• 对 $[0,1]$ 是对称的.

由此可以又把 $S$ 移空，然后进行下一步的迭代. 第三步开始的潜在最优区间就是 $\left[- \dfrac{1}{3},0\right]$ 和 $\left[0, \dfrac{1}{3}\right]$ . $\square$

多维 DIRECT 算法

将一维的 DIRECT 算法推广到高维，就是将区间变为超矩形，其中的一些定义发生了改变. 多维的中心点记为 $\boldsymbol{c}$ ，到各个顶点的距离定义为

$$ d = \frac{1}{2}\sqrt{\sum\limits_{i=1}^{n} s_{i}^{2}} $$

这里 $s_{i}$ 是第 $i$ 维的边长，$d$ 此时就是半对角线. 多维情况下，划分规则就要考虑最大边长，记 $\delta$ 为 $\dfrac{1}{3}$ 的最大边长，计算每个维度上的：

$$ w_{i} = \min \left\lbrace f(\boldsymbol{c}+ \delta \boldsymbol{e}_{i}), f(\boldsymbol{c}- \delta \boldsymbol{e}_{i}) \right\rbrace $$

这里 $i\in I$ ，$I$ 是具有最大边长的所有维度. 这个过程见下图：

然后整个算法的流程见下：

Julia 代码实现

我最近学了一下 Julia 代码，感觉很快，强类型，适合写数学算法类的代码. 因此本文算是一个练手，用于写真正的优化代码.

Shubert 算法

"""
shubert.jl

Shubert 算法的简单实现

@author: xzqbear
@time: 2026-03-05
"""

function shubert(
    f::Function,
    lb::Float64,
    ub::Float64,
    L::Float64,
    tol::Float64=1e-4,
    max_iter::Int64=1000
)::Tuple{Float64,Float64}
    """
    shubert(f, lb, ub, L, tol, max_iter)

    Description
        f: objective function
        lb: lower bound for varibale
        ub: lower bound for varibale
        L: Lipshitz constant
        tol: tolerance
        max_iter: maximum iterations
    """
    function shubert
        # function body
    end
    # 初始化
    X = Float64[lb, ub]
    F = Float64[f(lb), f(ub)]

    # V 型下界
    U(x::Float64, a::Float64, b::Float64) = min(f(a) - L * (x - a), f(b) + L * (x - b))
    current_best = minimum(F)

    # 在 for 之前创建数组
    # 如果在 for 循环内创建数组将会在 heap 上分配，速度较慢
    Z = Float64[]
    Uz = Float64[]

    for i in 1:max_iter
        l = length(X)

        empty!(Z)
        empty!(Uz)

        for j in 1:(l-1)
            z = (X[j] + X[j+1]) / 2 + (f(X[j]) - f(X[j+1])) / (2L)
            push!(Z, z)
            push!(Uz, U(z, X[j], X[j+1]))
        end

        # 采样点
        z_star = Z[argmin(Uz)]

        # 采样与更新
        insert_idx = searchsortedfirst(X, z_star)
        insert!(X, insert_idx, z_star)
        insert!(F, insert_idx, f(z_star))
        current_best = minimum(F)

        # 终止判断
        if abs(minimum(Uz) - current_best) < tol
            println("提前退出，共迭代 $i 次")
            break
        end
    end

    x_min = X[argmin(F)]
    f_min = minimum(F)

    return x_min, f_min
end

下面再给出个测试用例：

# 测试用例

ε = 1e-4
L = 4.0

f(x::Float64)::Float64 = (x - 2) * x * (x - 4)
a = 0.0
b = 4.0

@time shubert(f, a, b, L, ε)

Julia 很快就给出了结果：

提前退出，共迭代 456 次
  0.001535 seconds (162 allocations: 290.551 KiB, 7377.44% compilation time)
(3.1546970083418735, -3.0792014356348374)

REPL 也显示最终的极小值点和最小值为 (3.1546970083418735, -3.0792014356348374) ，和实际结果几乎一致.

需要注意这里的一些实现细节：

1. Julia 当中，for 循环内不要创建数组，否则会拖慢速度，Julia 的循环内创建数组是在 heap 上进行的，如果不是算法必须，请不要在循环里面创建数组.
2. 尽量多用函数而少使用全局变量，Julia 由于其动态编译语言的特性，类型安全很重要，但是全局变量无法推断其类型 ³，同时全局变量会在内存上运算，局部变量则是 Cache 上，速度有显著的差别，如果将上述实现用全局变量再跑一遍，几乎有 10 倍的差距.

DIRECT 算法

"""
direct.jl

DIRECT 算法的实现（基于论文《DIRECT算法及其实现》）
作者：xzqbear
时间：2026-03-06

主要特点：
- 使用无穷范数距离（最长边的一半）代替欧氏距离，简化计算
- 用分割次数指数存储边长，避免重复计算
- 分割时按采样点的最小函数值排序维度，将好点放入较大矩形
- 潜在最优矩形判定采用论文中的条件（结合凸包思想）
- 支持归一化搜索空间
"""

using Base: @kwdef

# 高维超矩形
@kwdef mutable struct Rectangle
    center::Vector{Float64}      # 归一化中心点坐标
    exponents::Vector{Int}        # 各维度的分割次数指数（初始为0，边长=3^{-exponent}）
    f_center::Float64             # 中心点函数值
    d::Float64                     # 无穷范数距离：0.5 * 3^{-max(exponents)}
    max_exp::Int                   # 最大指数（即最长边的指数）
end

# 多重派发
function Rectangle(
    center::Vector{Float64},
    exponents::Vector{Int},
    f_center::Float64
)::Rectangle
    """
    Rectangle

    给定中心与维度分割次数，返回矩形对象
    center: 归一化中心坐标
    exponents: 各维度分割次数
    f_center: 中心点函数值
    """
    max_exp = maximum(exponents)
    d = 0.5 * 3.0^(-max_exp)
    return Rectangle(center, exponents, f_center, d, max_exp)
end

# 归一化和逆归一化
normalize(x, lb, ub) = (x .- lb) ./ (ub .- lb)
denormalize(x_norm, lb, ub) = lb .+ x_norm .* (ub .- lb)


function direct(f, bounds;
    max_evals::Int=1000,
    max_iter::Int=500,
    ϵ::Float64=1e-4,
    tol::Float64=1e-6,
    verbose::Bool=false
)::Tuple{Array,Float64}
    """
    direct 优化方法

    max_evals: 最大函数 evaluation 次数
    max_iter: 最大迭代次数
    ϵ: Jones 因子
    tol: 容忍限
    verbose: 是否打印进度
    """

    # 获得维度和对应上下界
    n = length(bounds)
    lb = [b[1] for b in bounds]
    ub = [b[2] for b in bounds]

    # 函数式，包装标准化后的函数
    f_norm(x_norm) = f(denormalize(x_norm, lb, ub))

    # 初始化：整个单位超立方体作为一个矩形
    c0 = fill(0.5, n)                 # 中心点
    exp0 = zeros(Int, n)               # 初始分割次数为0
    f0 = f_norm(c0)
    rect0 = Rectangle(c0, exp0, f0)

    rectangles = [rect0]                # 所有矩形列表
    f_min = f0
    x_min_norm = copy(c0)

    n_evals = 1                         # 已进行的函数评估次数

    # 主循环
    for iter in 1:max_iter
        n_rect = length(rectangles)

        # 步骤1：找出所有潜在最优矩形的索引
        S = find_potentially_optimal(rectangles, f_min, ϵ)

        if isempty(S)
            verbose && println("迭代 $iter: 无潜在最优矩形，停止")
            break
        end

        # 步骤2：处理每个潜在最优矩形
        new_rects = Rectangle[]          # 临时存储本次划分产生的新矩形

        for idx in S
            rect = rectangles[idx]

            # 找出当前矩形的最大指数（即最长边）
            max_exp = rect.max_exp
            # 所有具有最大指数的维度
            I = findall(==(max_exp), rect.exponents)
            δ = 0.5 * 3.0^(-max_exp) * 2 / 3   # 采样步长：δ = (边长)/3，边长 = 2*d = 3^{-max_exp}
            # 更直接：边长 = 3^{-max_exp}，所以 δ = 边长/3 = 3^{-max_exp-1}
            δ = 3.0^(-max_exp - 1)

            # 采样：在每个最长边维度上取左右两点，记录 w_i = min(f_left, f_right)
            sampled = Dict{Int,Tuple{Vector{Float64},Float64,Vector{Float64},Float64}}()
            w = Dict{Int,Float64}()

            for i in I
                left = copy(rect.center)
                left[i] -= δ
                right = copy(rect.center)
                right[i] += δ
                f_left = f_norm(left)
                f_right = f_norm(right)
                n_evals += 2

                # 更新全局最优
                if f_left < f_min
                    f_min = f_left
                    x_min_norm = left
                end
                if f_right < f_min
                    f_min = f_right
                    x_min_norm = right
                end

                sampled[i] = (left, f_left, right, f_right)
                w[i] = min(f_left, f_right)
            end

            # 按 w 升序对维度排序（将好点优先放入大矩形）
            sorted_I = sort(collect(I), by=i -> w[i])

            # 开始划分：当前矩形将逐步缩小为中间矩形
            current_center = copy(rect.center)
            current_exps = copy(rect.exponents)

            for i in sorted_I
                left, f_left, right, f_right = sampled[i]

                # 左子矩形的指数：该维度指数+1，其他不变
                left_exps = copy(current_exps)
                left_exps[i] += 1
                left_rect = Rectangle(left, left_exps, f_left)
                push!(new_rects, left_rect)

                # 右子矩形
                right_exps = copy(current_exps)
                right_exps[i] += 1
                right_rect = Rectangle(right, right_exps, f_right)
                push!(new_rects, right_rect)

                # 更新当前矩形的指数（中间矩形）：该维度指数+1
                current_exps[i] += 1
                # 中心不变
            end

            # 更新原矩形为最终的中间矩形
            rect.center = current_center
            rect.exponents = current_exps
            # 重新计算 max_exp 和 d（通过构造函数或手动更新）
            # 这里我们手动更新字段，避免重新分配
            rect.max_exp = maximum(current_exps)
            rect.d = 0.5 * 3.0^(-rect.max_exp)
            # f_center 不变
        end

        # 将本次划分产生的新矩形添加到全局列表
        append!(rectangles, new_rects)

        # 检查终止条件
        # 条件1：函数评估次数超过上限
        if n_evals >= max_evals
            verbose && println("迭代 $iter: 达到最大函数评估次数 $max_evals")
            break
        end

        # 条件2：当前最优矩形（即 f_min 对应的矩形）的最长边小于容忍度
        # 找到 f_min 对应的矩形
        best_rect = nothing
        for rect in rectangles
            if rect.f_center == f_min  # 注意可能有多个，取第一个即可
                best_rect = rect
                break
            end
        end
        if best_rect !== nothing
            max_side = 3.0^(-best_rect.max_exp)  # 最长边长度
            if max_side < tol
                verbose && println("迭代 $iter: 最优矩形最长边 $max_side < $tol")
                break
            end
        end

        verbose && println("迭代 $iter: f_min=$f_min, 矩形数=$(length(rectangles)), 评估次数=$n_evals")
    end

    x_opt = denormalize(x_min_norm, lb, ub)
    return x_opt, f_min
end


function find_potentially_optimal(rectangles::Rectangle, f_min::Float64, ϵ::Float64)
    """
    find_potentially_optimal(rectangles, f_min, ϵ)

    查找潜在最优矩形的索引数组
    """

    n = length(rectangles)
    # 构建点列表 (d, f, idx)
    points = Vector{Tuple{Float64,Float64,Int}}(undef, n)
    for i in 1:n
        rect = rectangles[i]
        points[i] = (rect.d, rect.f_center, i)
    end

    # 按 d 排序
    sort!(points, by=x -> x[1])

    # 分组：每个 d 只保留最小 f 的代表点，同时记录该组所有索引
    groups = Vector{Tuple{Float64,Float64,Vector{Int}}}()
    i = 1
    while i ≤ n
        d = points[i][1]
        f_min_group = points[i][2]
        idxs = [points[i][3]]
        j = i + 1
        while j ≤ n && abs(points[j][1] - d) ≤ 1e-12
            if points[j][2] < f_min_group - 1e-12
                f_min_group = points[j][2]
                idxs = [points[j][3]]
            elseif abs(points[j][2] - f_min_group) ≤ 1e-12
                push!(idxs, points[j][3])
            end
            j += 1
        end
        push!(groups, (d, f_min_group, idxs))
        i = j
    end

    # 代表点列表 (d, f)
    rep = [(g[1], g[2]) for g in groups]

    # 计算右下凸包（使用斜率单调递增判断）
    stack = Tuple{Float64,Float64}[]
    for (d, f) in rep
        while length(stack) ≥ 2
            d1, f1 = stack[end-1]
            d2, f2 = stack[end]
            # 斜率 k1 = (f2-f1)/(d2-d1), k2 = (f-f2)/(d-d2)
            # 如果 k1 ≥ k2，则中间点 (d2,f2) 不在凸包上
            if (f2 - f1) / (d2 - d1) ≥ (f - f2) / (d - d2) - 1e-12
                pop!(stack)
            else
                break
            end
        end
        push!(stack, (d, f))
    end

    # 筛选满足第二条件的代表点
    S_indices = Int[]
    l = length(stack)
    for i in 1:l
        d, f = stack[i]
        # 找到对应的组
        group = nothing
        for g in groups
            if abs(g[1] - d) ≤ 1e-12 && abs(g[2] - f) ≤ 1e-12
                group = g
                break
            end
        end
        group === nothing && continue
        idxs = group[3]

        if i == length(stack)
            # 最后一个点，直接满足
            append!(S_indices, idxs)
        else
            d_next, f_next = stack[i+1]
            K = (f_next - f) / (d_next - d)
            lower_bound = f - K * d
            if lower_bound ≤ f_min - ϵ * abs(f_min) + 1e-12
                append!(S_indices, idxs)
            end
        end
    end

    return S_indices
end

测试用例可以考虑：

function test_goldstein_price()
    # Goldstein-Price 函数
    function gp(x)
        x1, x2 = x[1], x[2]
        term1 = 1 + (x1 + x2 + 1)^2 * (19 - 14x1 + 3x1^2 - 14x2 + 6x1 * x2 + 3x2^2)
        term2 = 30 + (2x1 - 3x2)^2 * (18 - 32x1 + 12x1^2 + 48x2 - 36x1 * x2 + 27x2^2)
        return term1 * term2
    end

    bounds = [[-2.0, 2.0], [-2.0, 2.0]]
    x_opt, f_opt = direct(gp, bounds, max_evals=1000, max_iter=100, ϵ=1e-4, tol=1e-6, verbose=false)
    println("Goldstein-Price 优化结果:")
    println("最优解: ", x_opt)
    println("最优值: ", f_opt)
    println("理论最优值 ≈ 3.0")
end

迭代出来还是速度比较快的.