Frank-Wolfe算法MATLAB实现：大规模约束优化的高效解法

简介：Frank-Wolfe算法是一种经典的约束凸优化方法，由J. Frank和D. Wolfe于1956年提出，在处理大型稀疏数据集时尤为高效，适合需要求解带约束目标函数最小化问题的场景。这份Matlab实现资源，面向正在学习优化算法原理、需要将理论转化为可调试代码的学生和算法工程师，问题定位明确。包内共3个文件，以txt和html为主，压缩后仅1KB；核心代码位于Frank-Wolfe (matlab).txt，完整实现了初始化、梯度计算、约束集边界上的最优方向搜寻、迭代更新以及终止判断等关键流程，另有html和txt各一份，可配合理解算法来源与Matlab编程细节。已有1026人学习/下载，说明其具备一定的参考价值；通过阅读并运行代码，读者既能对照理论步骤逐行理解Frank-Wolfe算法的工作机制，也能在此基础上扩展改造，适用于课程作业、算法对比实验或工程中的优化模块预研，是一份简洁实用的入门级Matlab优化程序。

1. Frank-Wolfe算法的MATLAB实现，为什么你在求解大规模约束优化时应该先看它

做优化的人都有这种体验：一个带线性约束和非光滑正则项的模型，手里已经有梯度，但投影算子算起来太贵。梯度下降每步都要把变量拉回可行域，投影一次可能就是一个二次规划。这时候Frank-Wolfe算法反而是最省事的——它不投影，只求解一个线性规划子问题，本质上是在做「沿着可行域顶点方向前进」。这个算法在机器学习、交通分配、张量分解和矩阵补全里都在用，MATLAB里实现一个完整版本不过几十行代码，比你想的简单得多。如果你在做一个目标函数可微、约束是凸多面体的优化问题，Frank-Wolfe应该是你考虑的第一个候选解法。本文将从一个工程视角带你从原理走到可运行的MATLAB程序，包含参数怎么设、问题怎么组织、结果怎么验证，拿到标题里的那个zip文件，你自己的版本也能照着跑通。

2. Frank-Wolfe的迭代格式与线搜索：先理解方向子问题和步长策略

2.1 Frank-Wolfe算法的核心迭代框架

Frank-Wolfe算法解决的是下面这个问题：

$$ \min_{x \in C} f(x) $$

其中$C$是一个凸多面体（单纯形、多面体、范数球都可以），$f$是连续可微的。算法每次迭代做三件事：计算梯度、在可行域上最小化线性函数、沿方向做线搜索。三件事都不涉及投影。

迭代格式写成伪代码就是：

初始化 x0 ∈ C for k = 0, 1, 2, ...: s_k = argmin_{s ∈ C} <∇f(x_k), s> # 方向子问题 d_k = s_k - x_k # Frank-Wolfe方向 γ_k = 最优步长 # 线搜索或衰减步长 x_{k+1} = x_k + γ_k * d_k

在MATLAB里的骨架对应如下。注意这里刻意把方向子问题单独写成函数，方便后续替换求解器。

function [x, fval, history] = frank_wolfe(f, grad_f, C_linprog, x0, opts) % F = objective handle, grad_f = gradient handle % C_linprog = struct with A, b, Aeq, beq, lb, ub for linprog % history = struct storing fval, iterate gap, gradient norm x = x0; history.fval = zeros(opts.maxiter, 1); history.gap = zeros(opts.maxiter, 1); for k = 1:opts.maxiter g = grad_f(x); s = solve_direction_subproblem(g, C_linprog); if norm(s - x, 2) < opts.tol % 见2.4，终止条件用迭代间隙 break end % 步长策略由opts.linesearch控制 gamma = select_stepsize(k, f, x, s, g, opts); x = x + gamma * (s - x); history.fval(k) = f(x); history.gap(k) = abs(norm(s - x, 2)); end end

代码的逻辑不复杂：grad_f(x)计算当前点梯度；solve_direction_subproblem求解线性方向子问题返回顶点$s$；select_stepsize返回步长；x = x + gamma * (s - x)完成更新。这个更新的几何含义是在当前点$x$和顶点$s$之间走一步，因此$x_k$永远留在可行域内部或边界上——这是Frank-Wolfe最大的优点：迭代点自然满足约束，不需要任何投影步骤。

2.2 方向子问题：为什么线性目标在凸多面体上永远有顶点解

方向子问题是整个算法里最重要的一步：

$$ s_k = \arg\min_{s \in C} \left\langle abla f(x_k), s \right\rangle $$

之所以这个子问题好求解，是因为线性目标函数在凸多面体上的最小值一定可以在某个顶点上取到。这个性质是最优性条件直接推论：若最优解是内部点，则梯度为零，否则最优点必然推至边界上某顶点。因此无论$C$是单纯形还是超立方体，只要你能枚举顶点或者用线性规划求解器找到最优顶点，方向子问题就解决了。

在MATLAB里最常见的求解方式是通过linprog。假设约束是多面体$Ax \le b$, $A_{eq}x = b_{eq}$，那么方向子问题的输入就是梯度向量：

function s = solve_direction_subproblem(g, C) % g: 当前梯度列向量 % C: 包含约束的结构体 f_lin = g; % linprog 求解 min g' * s optsLP = optimoptions('linprog', 'Display', 'off', 'Algorithm', 'dual-simplex'); [s, ~, exitflag] = linprog(f_lin, C.A, C.b, C.Aeq, C.beq, C.lb, C.ub, optsLP); if exitflag <= 0 error('Direction subproblem failed: exitflag=%d', exitflag); end end

参数说明：f_lin是线性目标系数，就是梯度向量本身；C.A和C.b对应不等式约束；C.Aeq和C.beq对应等式约束；lb和ub是变量边界。使用dual-simplex算法在线性规划问题中通常比内点法快，尤其在多面体约束矩阵是稀疏的大型问题上。每次迭代调用一次linprog是Frank-Wolfe的标准操作，在几百维到几千维的问题上，这个开销是完全可以接受的。

有个工程上的细节值得注意：linprog求出的解在数值上可能只是近似顶点，特别是当多面体有退化顶点时，返回点可能与理论顶点有微小偏差。这种误差在迭代中一般会被Frank-Wolfe的方向更新吸收掉，但如果你的约束条件有近似的线性依赖，建议在方向子问题里加入一个极小量的正则化，让数值更稳定。

2.3 步长规则：精确线搜索与衰减步长的取舍

Frank-Wolfe的步长选择有两个主流方式。

第一个是精确线搜索。沿方向$d_k = s_k - x_k$求最优步长：

$$ \gamma_k = \arg\min_{\gamma \in [0,1]} f(x_k + \gamma d_k) $$

这一步在$f$是二次函数时可以解析求解，在更复杂的目标函数上可以用MATLAB的fminbnd做一维搜索，每次迭代多花一点计算量，但能显著减少迭代次数。

第二个是通用衰减步长$\gamma_k = \frac{2}{k+2}$。这个步长不依赖目标函数的具体形式，理论收敛率$O(1/k)$在光滑强凸问题上是最优的。工程上的经验是：如果目标函数是良态的（比如最小二乘加线性约束），用固定衰减步长就足够；如果目标函数曲率变化剧烈，用衰减步长可能前期快后期慢，换成精确线搜索更好。

步长选择的MATLAB实现：

function gamma = select_stepsize(k, f, x, s, g, opts) switch opts.linesearch case 'decay' gamma = 2 / (k + 2); % 经典衰减步长，不需要函数值信息 case 'exact' d = s - x; fun = @(t) f(x + t * d); % 步长限定在[0,1]，因为x和s都在凸集C中，x + t*d是凸组合 [gamma, ~] = fminbnd(fun, 0, 1, optimset('TolX', 1e-10, 'Display', 'off')); otherwise error('Unknown linesearch strategy'); end end

注意fminbnd搜索区间要限制在$[0,1]$，因为当$t=1$时$x + d = s$，还在可行域里；$t$超过1意味着沿着$s$方向继续向前走，那可能跑出可行域，Frank-Wolfe的迭代点就不再保证可行了。这是实现时最容易犯的错误之一。

2.4 收敛性结论与终止条件

收敛性方面有一个工程上很实用的结论：如果$f$是光滑凸函数（梯度Lipschitz连续），步长采用$\gamma_k = \frac{2}{k+2}$，那么

$$ f(x_k) - f(x^*) \le \frac{2L D^2}{k+2} $$

其中$D$是可行域的直径，$L$是梯度的Lipschitz常数。这个界不依赖问题的维度，这一点和投影梯度法有本质区别，也是Frank-Wolfe在高维问题上受欢迎的原因之一。

终止条件最常用的是迭代间隙（duality gap的上界）。因为$f$是凸函数，有

$$ f(x) - f(x^*) \le \left\langle abla f(x), x - s \right\rangle $$

右端这个量恰好是每轮迭代可以免费算出来的。在代码里就是norm(s - x, 2) * norm(g, 2)的一个界。如果这个值小于容差，说明当前点已经足够接近最优了。实际实现种可以用：

gap = g' * (x - s); % 原始-对偶间隙的上界，非负 if gap < opts.tol break; end

这个量比观察目标函数值本身更可靠，因为目标函数值可能前期下降很快但后期卡住不动，gap能更真实地反映当前点距离最优解还有多远。工程上一般设opts.tol = 1e-6到1e-8即可；要求高精度时可能需要迭代几万次，此时衬度较大的问题应该考虑加加速度技术（见最后一章）或换用投影梯度类方法。

3. 在MATLAB中实现Frank-Wolfe主循环：从linprog到自写方向求解器

3.1 匿名函数传参与目标函数封装

MATLAB中处理优化问题的常用风格是让目标函数和梯度都以函数句柄传入。不要用全局变量传数据，在循环里多次调用时全局变量容易出错且难以调试。推荐把数据和参数封装在一个结构体里，或者直接用匿名函数捕获外部变量。

% 示例：目标函数为 f(x) = 0.5 * x'*Q*x + c'*x Q = [4, 1; 1, 2]; c = [-2; -1]; f = @(x) 0.5 * x' * Q * x + c' * x; grad_f = @(x) Q * x + c;

匿名函数的好处是记录在history里的函数值和梯度值都来自同一个封装，避免出现函数值计算和梯度计算不一致的笔误。梯度尽量解析求导，不要用数值差分；数值梯度一次误差放大就可能让方向子问题解出完全不同的顶点。

3.2 方向子问题的两种解法和取舍

如果可行域是标准单纯形$\Delta_n = {x \ge 0, \sum x_i = 1}$，方向子问题根本不需要linprog。线性函数在单纯形上的最小值点必然在某坐标轴方向上取到，直接把梯度向量的最小分量找出来，把质量集中在那个坐标上，其余为0即可。这是Frank-Wolfe被广泛应用在稀疏问题和概率分布估计上的原因之一。

function s = solve_direction_simplex(g) % g: 梯度向量 % 单纯形上最小化 <g, s>，解是单位向量，分量在g最小处为1 [~, idx] = min(g); s = zeros(size(g)); s(idx) = 1; end

如果约束是一个一般多面体，就用上一章给的linprog方案。取舍的标准是问题规模：单纯形或超立方体这类结构简单的约束，纯逻辑代码就能搞定，速度比linprog快几个数量级；一般多面体约束则建议直接用成熟线性规划求解器，因为自写单纯形法在退化情形下实现完整的两阶段法很费时间，而且数值稳定性上未必能超过linprog的实现。

一个折中方案是缓存linprog返回的活跃约束信息。如果多面体的约束矩阵不变，迭代过程中很多方向子问题是相似的，可以在结构体里保存上次的活跃集合，作为热启动传入linprog。MATLAB的linprog老版本不完全支持热启动，新版本可以通过initialpoint字段传入，具体看你的优化工具箱版本。

3.3 主循环代码：一个可以直接运行的完整版本

把上面所有部分拼起来，加上一个测试问题，就是一个20行左右可运行的Matlab脚本。

function test_fw_qp() % 测试 Frank-Wolfe: 最小化带框约束的凸二次函数 % min 0.5 * x'*Q*x - c'*x subject to 0 <= x_i <= ub_i rng(42); n = 50; Q = randn(n); Q = Q' * Q + eye(n); % 对称正定 c = randn(n, 1); ub = 2 * ones(n, 1); f = @(x) 0.5 * x' * Q * x - c' * x; grad_f = @(x) Q * x - c; % 用 linprog 求解方向子问题的约束结构 C.A = []; C.b = []; C.Aeq = []; C.beq = []; C.lb = zeros(n, 1); C.ub = ub; opts.maxiter = 10000; opts.tol = 1e-8; opts.linesearch = 'decay'; % 先用衰减步长试 x0 = zeros(n, 1); [x, fval, hist] = frank_wolfe(f, grad_f, C, x0, opts); fprintf('fval=%.8f, iter=%d\n', fval, length(find(hist.gap > 0))); end

主循环函数已经在2.1里给出，配合2.2、2.3两个子函数即可组成完整的程序。注意这里gap的第一轮会有值，后面如果步长选择合理，gap会持续下降；如果gap在某一轮突然变大，通常是数值误差累积或者linprog返回了非顶点解。

3.4 关键参数速查与调试技巧

下表给出常用参数的推荐范围和调整方向，能帮你把程序跑起来后快速收敛。

调试时最容易遇到的现象是gap下降到一个平台附近不再变化。此时先检查方向子问题是否真的返回了顶点，把s打印出来看是否在$C$内部；其次考虑步长序列是否衰减太快，改成exact线搜索试试；最后检查终止阈值是否太小，数值精度已经到极限了。

4. 验证与对比测试：单纯形约束、框约束和交通分配场景

4.1 测试一：投影到标准单纯形——FW充当一步近似投影

虽然标题没有提投影问题，但Frank-Wolfe在单纯形上的行为是检验实现正确与否最简洁的方法。把目标函数设为

$$ f(x) = \frac{1}{2}|x - y|_2^2 $$

可行域是标准单纯形，问题的最优解就是$y$在单纯形上的欧几里得投影。梯度是$x - y$，步长用精确线搜索。解析解可以通过专门算法（如project_simplex）算出，用来验证主循环是否正确。

function check_projection() y = [2; -1; 0.5; 1.2]; % 任意向量 n = length(y); Q = eye(n); c = -y; f = @(x) 0.5 * x' * Q * x - y' * x; grad_f = @(x) x - y; % 可行域：单纯形 sum(x)=1, x>=0 C.Aeq = ones(1,n); C.beq = 1; C.A = []; C.b = []; C.lb = zeros(n,1); C.ub = ones(n,1); % 冗余但无害 opts.maxiter = 1000; opts.tol = 1e-9; opts.linesearch = 'exact'; x0 = ones(n,1)/n; [x, fval] = frank_wolfe(f, grad_f, C, x0, opts); % 用投影法做交叉验证 x_proj = project_simplex(y); disp([x, x_proj]); % 两列应几乎相同 end function w = project_simplex(v) % 一种常见的单纯形投影算法 u = sort(v, 'descend'); sv = cumsum(u); rho = find(u .* (1:length(v))' > sv, 1, 'last'); theta = (sv(rho) - 1) / rho; w = max(v - theta, 0); end

如果两列结果一致，说明主循环和方向子问题实现正确。这个测试的价值在于：投影问题有一阶最优性条件可以直接验证——投影点到原点的连线与单纯形超平面法向一致。

4.2 测试二：目标函数带强线性项时，衰减步长为什么会慢

设$f(x) = \frac{1}{2}|x|^2 - 2x_1$，可行域是$[0,1]^2$。最优解应该在靠近$x_1=1$的边界远处。衰减步长$\gamma = 2/(k+2)$在迭代几百次后步长变得非常小，收敛速度明显变慢。此时改用精确线搜索可以显著减少迭代次数。

写个对比脚本：

function compare_linesearch() f = @(x) 0.5 * (x(1)^2 + x(2)^2) - 2 * x(1); grad_f = @(x) [x(1)-2; x(2)]; % 框约束 [0,1]x[0,1] C.A = []; C.b = []; C.Aeq = []; C.beq = []; C.lb = [0; 0]; C.ub = [1; 1]; x0 = [0; 0]; opts = struct('maxiter', 200, 'tol', 1e-8, 'linesearch', 'decay'); [~, ~, h1] = frank_wolfe(f, grad_f, C, x0, opts); opts.linesearch = 'exact'; [~, ~, h2] = frank_wolfe(f, grad_f, C, x0, opts); semilogy(1:200, h1.fval, 'r-', 1:200, h2.fval, 'b--'); legend('decay', 'exact'); end

线性项占主导时，衰减步长会让迭代点近似沿着一条路径缓慢爬行，而精确线搜索每一步几乎都走到底再折线接近最优解。工程经验是：可行域维度高、约束简单时用decay；维度低、可以接受一维搜索开销时用exact。

4.3 测试三：交通网络用户均衡分配（FW经典战场）

Frank-Wolfe最早的大规模应用场景是交通分配问题。用户均衡分配求解的是

$$ \min_{f} \sum_{a} \int_0^{f_a} t_a(w) dw $$

约束是路径流之和等于出发地到目的地的需求。目标函数可微且约束是多面体，非常适合FW。在MATLAB里做一个简单路网测试不需要特别复杂，核心是你有一个路段时间函数$t_a(f_a)$，一次路径流分配计算就是一次方向子问题。这里不再贴完整代码，但提醒一点：方向子问题的解在这个语境下是「全有全无分配」（将所有需求分配到最短路径），每次迭代都是在最短路径方向和当前路径流之间做一个凸组合。这个视角是Frank-Wolfe算法在交通领域如此流行的原因。

4.4 三个常见错误与避坑

第一，线搜索区间超过$[0,1]$导致迭代点出厂。第二，将梯度方向误用为目标函数值本身，方向子问题的输入必须是梯度向量。第三，linprog的约束矩阵忘记把变量边界统一处理，导致方向子问题无界或不可行。这些错误的共同信号是：gap不为正数或迭代点出现NaN。一条调试路径：先打印每轮s和gamma，检查s是否在约束内、gamma是否在$[0,1]$内，用断言把这两个条件写进去，代码安全很多。

5. 一个收尾技巧：用KKT残差检验你的Frank-Wolfe实现是否真的收敛

写最后一个实用工具。前面用的gap只能说明当前点距离最优点的目标函数值差已经被控制，但不能说明你找到了精确的约束最优解。更强的一个验证手段是检查KKT条件。对问题$\min f(x), s.t. Ax \le b$，KKT条件说存在乘子$\lambda \ge 0$使得梯度满足

$$ abla f(x) + A^T\lambda = 0, \quad \lambda_i (b_i - A_i x) = 0 $$

在MATLAB里可以算一遍：

function res = kkt_residual(x, grad, A, b) % 计算KKT残差 % 找出活跃约束（等号近似成立） active_iter = abs(A * x - b) < 1e-6; A_active = A(active_iter, :); b_active = b(active_iter); % 最小二乘求解乘子：最小化 ||grad + A_active'*lambda||^2, lambda>=0 options = optimoptions('lsqlin', 'Display', 'off', 'Algorithm', 'active-set'); lambda = lsqlin(A_active', -grad, [], [], [], [], zeros(sum(active_iter),1), [], [], options); res = norm(grad + A_active'*lambda, inf); end

这个残差越低，说明当前点越接近真正的KKT点。对Frank-Wolfe的终止可以同时看两个量：gap < tol和这个KKT残差。两者同时满足，再谈收敛就比较可信了。这个验证方法在你拿到别人的Frank-Wolfe算法 matlab程序.zip时特别有用，第一步跑通，第二步把约束和目标函数换成自己的问题，第三步用这个残差验收，比单纯看目标函数曲线可靠得多。配合把history.fval导出成CSV画对数坐标图，整个验证链路就齐了。

本文还有配套的精品资源，点击获取