组合优化与凸优化相关算法 Julia实现

线性规划

方法有单纯形法（简单，非多项式），椭圆法（复杂，多项式，仅有理论价值），内点法（非多项式，实际效率高）。

以例子说明，目标函数
$$\min z=0x_1+0.1x_2+0.2x_3+0.3x_4+0.8x_5$$
约束条件
$$s.t.\{\begin{array}{l}{x}_1+2x_2+x_4=100\\ 2x_3+2x_4+x_5=100\\ 3x_1+x_2+2x_3+3x_5=100\\ x_1,x_2,x_3,x_4,x_5\ge 0\end{array}$$

代码：

using JuMP, GLPK
m = Model(GLPK.Optimizer)
@variable(m, x[1:5] >= 0, integer = true)
@objective(m, Min, 0.1x[2] + 0.2x[3] + 0.3x[4] + 0.8x[5])
@constraint(m, constraint1, x[1] + 2x[2] + x[4] == 100)
@constraint(m, constraint2, 2x[3] + 2x[4] + x[5] == 100)
@constraint(m, constraint3, 3x[1] + x[2] + 2x[3] + 3x[5] == 100)
#print(m)
JuMP.optimize!(m)
#println("Optimal Soutions:")
for i in 1:5
    println("x$i = ", JuMP.value(x[i]))
end
println("z = ", JuMP.objective_value(m))

结果：

x1 = 30.0
x2 = 10.0
x3 = 0.0
x4 = 50.0
x5 = 0.0
z = 16.0

一维搜索

黄金分割和 Fibonacci 法

取 $\alpha <\lambda <\mu <\beta$ 和 $x\in [\alpha ,\beta ]$ 上单谷函数 $\Phi (x)$，求其最值 $x^{\ast }$，则有：
一，若 $\Phi (\lambda )\ge \Phi (\mu )$，则 $x^{\ast }\in [\lambda ,\beta ]$；
二，若 $\Phi (\lambda )<\Phi (\mu )$，则 $x^{\ast }\in [\alpha ,\mu ]$；
如此迭代缩减区间长度，进行精确一维搜索。

对称性推导

通过选择不同的 $\lambda ,\mu$ 有不同策略，而一个好的策略保证对称性 $\lambda -\alpha =\beta -\mu$，下面推导之。

设归一化区间参数 $p=\frac{\lambda -\alpha }{\beta -\alpha },q=\frac{\beta -\mu }{\beta -\alpha }$，则任一个策略即可表示为点对 $(p,q)$，且有 $p,q\in (0,1)$ 和 $p+q<1$，不妨设 $x^{\ast }$ 在 $[\alpha ,\beta ]$ 上均匀分布，则分以下三种情况：
一，$x^{\ast }\in [\alpha ,\lambda ]$，概率为 $p$，此时必然有 $\Phi (\lambda )<\Phi (\mu )$，下轮区间缩减为 $[\alpha ,\mu ]$；
二，$x^{\ast }\in [\mu ,\beta ]$，概率为 $q$，此时必然有 $\Phi (\lambda )>\Phi (\mu )$，下轮区间缩减为 $[\lambda ,\beta ]$；
三，$x^{\ast }\in [\lambda ,\mu ]$，概率为 $1-p-q$，此时不妨设 $\Phi (\lambda )<\Phi (\mu )$ 和 $\Phi (\lambda )>\Phi (\mu )$ 发生的概率各占一半。

用缩减比 $t(p,q)$ 来表示策略的好坏，可由下轮区间长度与本轮区间长度的比值得到。综上，可分为两种情况：
情况一，下轮区间缩减为 $[\alpha ,\mu ]$ 的概率为 $\frac{1+p-q}{2}$，缩减比为 $1-q$；
情况二，下轮区间缩减为 $[\lambda ,\beta ]$ 的概率为 $\frac{1+q-p}{2}$，缩减比为 $1-p$；
则有
$$t(p,q)=(1-q)\frac{1+p-q}{2}+(1-p)\frac{1+q-p}{2}=\frac{2-(p+q)+(p-q{)}^2}{2}$$
而好的策略意味着低的缩减比。

可得在 $p+q$ 一定的条件下，好的策略保证 $p=q$，即对称性 $\lambda -\alpha =\beta -\mu$。

黄金缩减比推导

经过以上讨论，若只考虑缩减比最优，那 $p=q=\frac{1}{2}-ϵ$ 无疑是最优策略。可实际上，计算 $\Phi (\lambda )$ 和 $\Phi (\mu )$ 在大多数情况下都是该算法流程中的主要时间消耗，而区间缩减为 $[\alpha ,\mu ]$ 的情况下 $\Phi (\lambda )$ 是可以作为下轮计算的已知量的，同理区间缩减为 $[\lambda ,\beta ]$ 的情况下 $\Phi (\mu )$ 也是可以作为下轮计算的已知量。

换言之，结果区间精度与初始区间精度的比值为 $r$，在最优缩减比 $t^{\ast }=\frac{1}{2}+ϵ$ 下每轮计算两次 $\Phi (x)$，时间消耗为 $T(t^{\ast })=2{\log }_{t^{\ast }}r$，若存在缩减比 $t^{\prime }$ 下每轮计算一次 $\Phi (x)$ 的策略，其时间消耗为 $T(t^{\prime })={\log }_{t^{\prime }}r$，只需要满足条件 $t^{\prime }<\sqrt{t^{\ast }}=\frac{1}{\sqrt{2}}+ϵ$，即 $t^{\prime }\le \frac{1}{\sqrt{2}}\approx 0.707$ 即可保证该策略比最优缩减比策略更好。

对于本轮区间参数 $(\alpha ,\lambda ,\mu ,\beta )$ 和下轮区间参数 $({\alpha }^{\prime },{\lambda }^{\prime },{\mu }^{\prime },{\beta }^{\prime })$，保证缩减比
$$t=\frac{\mu -\alpha }{\beta -\alpha }=\frac{\beta -\lambda }{\beta -\alpha }=\frac{{\mu }^{\prime }-{\alpha }^{\prime }}{{\beta }^{\prime }-{\alpha }^{\prime }}=\frac{{\beta }^{\prime }-{\lambda }^{\prime }}{{\beta }^{\prime }-{\alpha }^{\prime }}$$
的情况下，希望有如下性质：
当 ${\alpha }^{\prime }=\alpha ,{\beta }^{\prime }=\mu$ 时，有 ${\mu }^{\prime }=\lambda$；（${\lambda }^{\prime }=\lambda$ 是无法满足的）
当 ${\alpha }^{\prime }=\lambda ,{\beta }^{\prime }=\beta$ 时，有 ${\lambda }^{\prime }=\mu$；（${\mu }^{\prime }=\mu$ 是无法满足的）
由于左右两种情况对称，故以 ${\alpha }^{\prime }=\alpha ,{\mu }^{\prime }=\lambda ,{\beta }^{\prime }=\mu$ 为例，设此时缩减比为 $t_{\varphi }$，则有
$$\begin{gathered} \beta -\lambda =\mu -\alpha =t_{\varphi }(\beta -\alpha ) \\ \lambda -\alpha ={\mu }^{\prime }-{\alpha }^{\prime }=t_{\varphi }({\beta }^{\prime }-{\alpha }^{\prime })=t_{\varphi }(\mu -\alpha )=t_{\varphi }^2(\beta -\alpha ) \end{gathered}$$
得方程
$$\beta -\alpha =(\beta -\lambda )+(\lambda -\alpha )=t_{\varphi }(\beta -\alpha )+t_{\varphi }^2(\beta -\alpha )$$
化简之，即
$$t_{\varphi }^2+t_{\varphi }-1=0$$
由 $t_{\varphi }>0$ 舍去负根，得 $t_{\varphi }=\frac{\sqrt{5}-1}{2}\approx 0.618$ 恰为黄金分割比，满足上文条件为最优策略。此即一维搜索的黄金分割法。

Fibonacci 法

设 Fibonacci 数列 $F_0=0,F_1=1,F_n=F_{n-1}+F_{n-2}$，由于 ${\lim }_{n\to \infty }\frac{F_n}{F_{n+1}}=\frac{\sqrt{5}-1}{2}$，故可用反向 Fibonacci 数列指导每轮缩减，是为 Fibonacci 法，代码如下：

function FibonacciSection(f, l, u, tol)
    fib = [1, 1]
    while fib[end] * tol < u - l
        push!(fib, fib[end] + fib[end-1])
    end
    m, n = (fib[end-1] * l + fib[end-2] * u) / fib[end], (fib[end-2] * l + fib[end-1] * u) / fib[end]
    for i = length(fib)-1 : -1 : 3
        if f(m) > f(n)
            l, m = m, n
            n = (fib[i-2] * l + fib[i-1] * u) / fib[i]
        else
            u, n = n, m
            m = (fib[i-1] * l + fib[i-2] * u) / fib[i]
        end
    end
    return m
end

二分法

取当前区间 $[\alpha ,\beta ]$ 中点 $\lambda =\frac{\alpha +\beta }{2}$，求导数值 ${\Phi }^{\prime }(\lambda )$，并判断：
若 ${\Phi }^{\prime }(\lambda )>0$，区间更新为 $[\alpha ,\lambda ]$；
若 ${\Phi }^{\prime }(\lambda )<0$，区间更新为 $[\lambda ,\beta ]$；
若 ${\Phi }^{\prime }(\lambda )=0$，$\lambda$ 为最小点，算法结束。

代码：

function BinarySection(f, grad_f, l, u, tol)
    m = (l + u) / 2
    while u - l > tol
        g = grad_f(m)
        if g > 0.
            u = m
        elseif g < 0.
            l = m
        else
            break
        end
        m = (l + u) / 2
    end
    return m
end

Shubert-Piyavskii 法

若有 $\left|\Phi (\lambda )-\Phi (\mu )\right|\le l\left|\lambda -\mu \right|,\forall \lambda ,\mu \in [a,b]$，则称函数 $\Phi (x)$ 在 $[a,b]$ 上是 $l$-Lipschitz 连续的，即导数的幅度存在上界。

Shubert-Piyavskii 法用斜率为 $\pm l$ 的线段，迭代构造 $\Phi (x)$ 的锯齿状下界，代码：

function Shubert_Piyavskii(f, L, l, u, tol)
    intersection(p, q) = ((p[1] + q[1]) / 2 + (p[2] - q[2]) / (2 * L), (p[1] - q[1]) * L / 2 + (p[2] + q[2]) / 2)
    m = (l + u) / 2
    Up = [(l, f(l)), (m, f(m)), (u, f(u))]
    Lo = [intersection(Up[1], Up[2]), intersection(Up[2], Up[3])]
    f_min = minimum(map(last, Up))
    while true
        m = argmin(map(last, Lo))
        y = f(Lo[m][1])
        if f_min > y
            f_min = y
            if (f_min - Lo[m][2]) * 2 < L * tol
                break
            end
        end
        Up = vcat(Up[1 : m], [(Lo[m][1], y)], Up[m+1 : end])
        Lo = vcat(Lo[1 : m-1], [intersection(Up[m], Up[m+1]), intersection(Up[m+1], Up[m+2])], Lo[m+1 : end])
    end
    return Lo[m][1]
end

不精确搜索

考虑从 $x^{(k)}$ 出发，沿方向 $d^{(k)}$ 寻找新迭代点 $x^{(k+1)}=x^{(k)}+s^{(k)}=x^{(k)}+{\lambda }_k{d}^{(k)}$ 使得 $\Phi (x^{(k+1)})$ 尽可能小。选取步长增大系数 $\alpha >1$ 和步长缩短系数 $\beta <1$ 迭代求步长 ${\lambda }_k$，直到满足给定规则。

代码：

macro LineSearch(cond2)
    quote
        alpha, beta, lambda = 1.5, 0.5, 1.
        while true
            g, s = grad_f(x), lambda * d
            if f(x + s) - f(x) > rho * g' * s
                lambda = beta * lambda
            elseif $cond2
                lambda = alpha * lambda
            else
                break
            end
        end
        return lambda
    end |> esc
end

Goldstein(f, grad_f, x, d, rho)         = @LineSearch(f(x + s) - f(x) < (1. - rho) * g' * s)
Armijo(f, grad_f, x, d, rho, mu)        = @LineSearch(f(x + s) - f(x) < (mu * rho) * g' * s)
Goldstein2(f, grad_f, x, d, rho, sigma) = @LineSearch(f(x + s) - f(x) < sigma * g' * s)
WP(f, grad_f, x, d, rho, sigma)         = @LineSearch(grad_f(x + s)' * s < sigma * g' * s)
WP2(f, grad_f, x, d, rho, eta)          = @LineSearch(abs(grad_f(x + s)' * d) > -eta * g' * d)

Goldstein 规则

(1) $\Phi (x^{(k+1)})-\Phi (x^{(k)})\le \rho \nabla {\Phi }^T(x^{(k)})s^{(k)}$
(2) $\Phi (x^{(k+1)})-\Phi (x^{(k)})\ge (1-\rho )\nabla {\Phi }^T(x^{(k)})s^{(k)}$
其中 $\rho \in (0,\frac{1}{2})$，实际取 0.1 或更小。

Armijo 规则

(1) 同 Goldstein(1)
(2) $\Phi (x^{(k+1)})-\Phi (x^{(k)})\ge \mu \nabla {\Phi }^T(x^{(k)})s^{(k)}$
其中 $\mu \in (5,10)$。

Goldstein 改进规则

(1) 同 Goldstein(1)
(2) $\Phi (x^{(k+1)})-\Phi (x^{(k)})\ge \sigma \nabla {\Phi }^T(x^{(k)})s^{(k)}$
其中 $\sigma \in (\rho ,1)$。

Wolfe-Powell 规则

(1) 同 Goldstein(1)
(2) $\nabla {\Phi }^T(x^{(k+1)})s^{(k)}\ge \sigma \nabla {\Phi }^T(x^{(k)})s^{(k)}$
其中 $\sigma \in (\rho ,1)$。

Wolfe-Powell 改进规则

(1) 同 Goldstein(1)
(2) $\left|\nabla {\Phi }^T(x^{(k+1)})d^{(k)}\right|\le -\eta \nabla {\Phi }^T(x^{(k)})d^{(k)}$
其中 $\eta \in (0,1)$，当 $\eta =0$ 即为精确一维搜索，当 $\eta =0.1$ 可看作近似的精确一维搜索。

无约束最优化

无约束最优化问题 $\min f(x),f:{\mathbb{R}}^n\to \mathbb{R}$ 往往采用迭代下降的方法求解，即 $x^{(k+1)}=x^{(k)}+{\lambda }_k{d}^{(k)}$，不同算法间的区别在于求得 $d^{(k)}$ 的方式和得出的值不同。

当 $f(x)$ 为二次函数时，写成二次型 $f(x)=\frac{1}{2}{x}^{T}Ax+b^{T}x+c,\nabla f(x)=Ax+b$ 带入方程
$$\nabla f^T(x+\lambda d)d=0$$
得
$$(A(x+\lambda d)+b{)}^{T}d=0$$
化简得最佳步长 $\lambda$ 为
$$\lambda =-\frac{d^T(Ax+b)}{d^{T}Ad}$$
若 $f(x)$ 非二次函数，则需要线搜索得到最佳步长。

有约束最优化问题可通过罚因子方法或增广 Lagrange 方法转化为无约束最优化问题求解。

测试选用 Rosenbrock Banana 函数
$$f(x)=(a-x_1{)}^2+b(x_2-x_1^2{)}^2$$
其有全局极小点 $(a,a^2)$，取 $a=1,b=5$，初始解 $x_0=(1,-1{)}^T$。代码：

using LinearAlgebra, Optim
using Plots

macro iterDescent(part0, part1, part2)
    quote
        x, r = x0, [x0]
        $part0
        for i = 1 : max_iter
            $part1
            if norm(d) < tol
                break
            end
            a = Optim.minimizer(optimize(a -> f.f(x + a * d), -1., 1.))
            s = a * d
            if norm(s) < tol
                break
            end
            x = x + s
            $part2
            push!(r, x)
        end
        return r
    end |> esc
end

function draw(p)
    for i = 2 : length(p)
        plot!([p[i - 1][1], p[i][1]], [p[i - 1][2], p[i][2]], color = :black, legend = false, w = 1)
    end
end

function RosenbrockBanana(a, b)
    (
        f = x::Vector{Float64} -> (a - x[1])^2 + b * (x[2] - x[1]^2)^2,
        g = x::Vector{Float64} -> [-2 * (a - x[1]) - 4 * b * x[1] * (x[2] - x[1]^2); 2 * b * (x[2] - x[1]^2)],
        h = x::Vector{Float64} -> [(2 - 4 * b * (x[2] - 3 * x[1]^2)) (-4 * b * x[1]); (-4 * b * x[1]) (2 * b)]
    )
end

rb = RosenbrockBanana(1, 5)
contour(range(-2, 2, length = 1000), range(-2, 2, length = 1000), (x, y) -> rb.f([x, y]),
    aspect_ratio = 1, levels = 100, color=:turbo, xlims = (-2, 2), ylims = (-2, 2))
x0 = [1., -1.]

坐标轮换法

迭代方向为坐标轴，即单位基向量，每 $n$ 轮一循环。

代码：

CyclicCoordinate(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent(
    n = length(x0),
    begin
        d = zeros(n); d[(i - 1) % n + 1] = 1
    end,
    nothing
)
draw(CyclicCoordinate(rb, x0))

结果：

Hooke-Jeeves 法

在坐标轮换法基础上加上第 $n+1$ 轮的加速步，方向为 $d^{(k)}=x^{(k)}-x^{(k-n-1)}$，每 $n+1$ 轮一循环。

代码：

HookeJeeves(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent(
    begin
        n = length(x0); l = zeros(n)
    end,
    begin
        if i % (n + 1) > 0
            d = zeros(n); d[i % (n + 1)] = 1
        else
            d = l
        end
    end,
    begin
        if i % (n + 1) > 0
            l[i % (n + 1)] += a
        else
            l = s
        end
    end
)
draw(HookeJeeves(rb, x0))

结果：

Rosenbrock 法

第 $1$ 到 $n$ 轮同坐标轮换法，然后将前 $n$ 轮的位移 $d^{(k)}=x^{(k)}-x^{(k-n)}$ 作为第 $n+1$ 轮的方向，并进行 Gram-Schmidt 正交化形成新的基，每 $n$ 轮一循环。
$$a_j=\{\begin{array}{ll}{d}_j,& {\lambda }_j=0\\ {\Sigma }_{i=j}^n{\lambda }_i{d}_i,& {\lambda }_j\ne 0\end{array},b_j=\{\begin{array}{ll}{a}_j,& j=1\\ a_j-{\Sigma }_{i=1}^{j-1}(a_j^T{\overline{d}}_i){\overline{d}}_i,& j\ge 2\end{array},{\overline{d}}_j=\frac{b_j}{\Vert b_j\Vert }$$

代码：

Rosenbrock(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent(
    begin
        n = length(x0); l = zeros(n)
        b = []
        for i = 1 : n
            push!(b, zeros(n))
            b[end][i] = 1
        end
    end,
    d = b[(i - 1) % n + 1],
    begin
        l[(i - 1) % n + 1] = a
        if i % n == 0
            t = zeros(n)
            for j = n : -1 : 1
                if l[j] != 0
                    t += l[j] * b[j]
                    b[j] = t
                end
            end
            for j = 1 : n
                t = zeros(n)
                for k = 1 : j - 1
                    t += (b[j]' * b[k]) * b[k]
                end
                b[j] -= t
                b[j] /= norm(b[j])
            end
            l = zeros(n)
        end
    end
)
draw(Rosenbrock(rb, x0))

结果：

梯度下降法

迭代方向为梯度下降方向。（二维和坐标轮换法类似，更高维情况有不同）

代码：

GradientDescent(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent nothing d = -f.g(x) nothing
draw(GradientDescent(rb, x0))

结果：

Newton 法

在当前解局部用二次函数近似，求其梯度方向 $g_k$ 和 Hessian 矩阵 $G_k$ 的逆 $H_k=G_k^{-1}$，将 $-H_k{g}_k$ 作为迭代方向。

代码：

Newton(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent nothing d = -f.h(x) \ f.g(x) nothing
draw(Newton(rb, x0))

结果：

DFP 法

考虑到 Hessian 矩阵及其逆不易求得，故用公式近似 Hessian 矩阵的逆 $H_k$，这类方法为拟 Newton 法，有 DFP 法和 BFGS 法，区别在于近似公式不同。

$$H^{k+1}=H^k-\frac{H^k{y}^k(H^k{y}^k{)}^T}{(y^k{)}^T{H}^k{y}^k}+\frac{s^k(s^k{)}^T}{(y^k{)}^T{s}^k}$$

代码：

DFP(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent(
    begin
        g, H = f.g(x0), Matrix{Float64}(I, length(x0), length(x0))
    end,
    d = -H * g,
    begin
        g_new = f.g(x)
        y, g = g_new - g, g_new
        H = H - ((H * y) * (H * y)') / (y' * H * y) + (s * s') / (y' * s)
    end
)
draw(DFP(rb, x0))

结果：

BFGS 法

$$H^{k+1}=(I-{\rho }_k{y}^k(s^k{)}^T{)}^T{H}^k(I-{\rho }_k{y}^k(s^k{)}^T)+{\rho }_k{s}^k(s^k{)}^T$$
其中 ${\rho }_k=\frac{1}{(s^k{)}^T{y}^k}$

代码：

BFGS(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent(
    begin
        g, H, E = f.g(x0), Matrix{Float64}(I, length(x0), length(x0)), Matrix{Float64}(I, length(x0), length(x0))
    end,
    d = -H * g,
    begin
        g_new = f.g(x)
        y, g = g_new - g, g_new
        rho = 1 / (s' * y)
        H = (E - rho * y * s')' * H * (E - rho * y * s') + rho * s * s'
    end
)
draw(BFGS(rb, x0))

结果：

Fletcher-Reeves 共轭梯度法

若 $d_{k+1}^T{G}_k{d}_k=0$，则称 $d_{k+1}$ 是共轭方向。

代码：

ConjugateGradientFR(f, x0, tol = 1e-6, max_iter = 1000) = @iterDescent(
    begin
        g = f.g(x0); d = -g
    end,
    nothing,
    begin
        g_new = f.g(x)
        b, g = (g_new' * g_new) / (g' * g), g_new
        d = - g_new + b * d
    end
)
draw(ConjugateGradientFR(rb, x0))

结果：

References

北大李东风 Julia语言入门 教材：https://www.math.pku.edu.cn/teachers/lidf/docs/Julia/html/_book/index.html
快速入门: 一份简单而粗略的语言概览 ：https://cheatsheet.juliadocs.org/zh-cn/
最优化：建模、算法与理论 (刘浩洋等)