基于c++ eigen的Nelder-Mead算法（仿照scipy）
·
Nelder-Mead算法简介
Nelder-Mead算法（也称单纯形法）是一种常用的无导数优化方法，特别适合于：
- 函数导数难以计算或不存在的情况
- 寻找非线性函数的局部最小值
- 维度适中的问题
这种算法被广泛应用于机器学习、计算机视觉、信号处理等领域，特别是当目标函数计算成本高昂时尤为有用。
SciPy版本的特点：
- 边界约束支持：SciPy版本添加了对变量边界的支持，通过bounds参数实现。它使用np.clip确保所有单纯形顶点都在指定边界内，这是标准Nelder-Mead算法的扩展。
- 自适应参数策略：通过adaptive参数启用，基于Gao和Han
2012年的论文。这种策略根据问题维度自动调整反射、扩展和收缩参数，使算法在高维问题上表现更好。 - 双重收敛条件：同时使用函数值容差(fatol)和参数值容差(xatol)作为终止条件，这比只用单一条件更稳健。
- 完善的终止处理：支持通过最大迭代次数(maxiter)或最大函数评估次数(maxfev)限制算法运行时间，并提供详细的终止状态。
更详细的分析见本文结尾。
代码
纯头文件库，包含opti.h,opti.inl。还有一个测试程序main.cpp
opti.h:
#ifndef MY_OPTIMIZE_H |
#define MY_OPTIMIZE_H |
#include <Eigen/Core> |
#include <algorithm> |
#include <cmath> |
#include <functional> |
#include <iostream> |
#include <limits> |
#include <vector> |
// NelderMead算法选项 |
struct NelderMeadOptions { |
double xatol = 1e-4; // x收敛容差 |
double fatol = 1e-4; // 函数值收敛容差 |
int maxiter = -1; // 最大迭代次数，-1表示使用默认值 |
int maxfev = -1; // 最大函数评估次数，-1表示使用默认值 |
bool disp = false; // 是否显示收敛消息 |
bool return_all = false; // 是否返回所有迭代的结果 |
bool adaptive = false; // 是否使用自适应参数 |
// 边界，每个变量的(min, max)对，std::nullopt表示无边界 |
std::vector<std::pair<double, double>> bounds; |
}; |
// NelderMead算法结果 |
struct NelderMeadResult { |
Eigen::VectorXd x; // 最优解 |
double fun; // 最优解对应的函数值 |
int nit; // 迭代次数 |
int nfev; // 函数评估次数 |
int status; // 状态码：0=成功，1=达到最大函数评估次数，2=达到最大迭代次数 |
bool success; // 是否成功 |
std::string message; // 状态信息 |
std::vector<Eigen::VectorXd> allvecs; // 所有迭代的结果（如果return_all=true） |
}; |
/** |
* Nelder-Mead单纯形算法 |
* |
* @param func 目标函数 |
* @param x0 初始猜测点 |
* @param initial_simplex 初始单纯形，如果提供则覆盖x0 |
* @param options 算法参数选项 |
* @param callback 回调函数，每次迭代后调用 |
* @return 优化结果 |
*/ |
template <typename Scalar = double> |
NelderMeadResult |
minimize_neldermead(std::function<Scalar(const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>&)> func, |
const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& x0, |
const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>* initial_simplex = nullptr, |
const NelderMeadOptions& options = NelderMeadOptions(), |
std::function<bool(const NelderMeadResult&)> callback = nullptr); |
/** |
* 简化版Nelder-Mead优化函数，类似于scipy的fmin |
* |
* @param func 目标函数 |
* @param x0 初始猜测点 |
* @param initial_simplex 初始单纯形，如果提供则覆盖x0 |
* @param xtol x收敛容差 |
* @param ftol 函数值收敛容差 |
* @param maxiter 最大迭代次数 |
* @param maxfun 最大函数评估次数 |
* @param full_output 是否返回完整输出 |
* @param disp 是否显示收敛消息 |
* @param retall 是否返回所有迭代的结果 |
* @param callback 回调函数 |
* @return 优化结果点或完整结果（取决于full_output） |
*/ |
template <typename Scalar = double> |
Eigen::Matrix<Scalar, Eigen::Dynamic, 1> |
fmin(std::function<Scalar(const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>&)> func, |
const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& x0, |
const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>* initial_simplex = nullptr, |
double xtol = 1e-4, |
double ftol = 1e-4, |
int maxiter = -1, |
int maxfun = -1, |
bool full_output = false, |
bool disp = true, |
bool retall = false, |
std::function<bool(const NelderMeadResult&)> callback = nullptr); |
#include "opti.inl" |
#endif //MY_OPTIMIZE_H |
opti.inl:
#ifndef MY_OPTIMIZE_INL |
#define MY_OPTIMIZE_INL |
#include <algorithm> |
#include <cmath> |
#include <iostream> |
template <typename Scalar> |
NelderMeadResult minimize_neldermead(std::function<Scalar(const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>&)> func, |
const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& x0, |
const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>* initial_simplex, |
const NelderMeadOptions& options, |
std::function<bool(const NelderMeadResult&)> callback) { |
using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>; |
using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>; |
// 提取选项参数 |
double xatol = options.xatol; |
double fatol = options.fatol; |
int maxiter = options.maxiter; |
int maxfev = options.maxfev; |
bool disp = options.disp; |
bool return_all = options.return_all; |
bool adaptive = options.adaptive; |
const auto& bounds = options.bounds; |
// 设置算法参数 |
double rho, chi, psi, sigma; |
if (adaptive) { |
double dim = static_cast<double>(x0.size()); |
rho = 1.0; |
chi = 1.0 + 2.0 / dim; |
psi = 0.75 - 1.0 / (2.0 * dim); |
sigma = 1.0 - 1.0 / dim; |
} else { |
rho = 1.0; |
chi = 2.0; |
psi = 0.5; |
sigma = 0.5; |
} |
double nonzdelt = 0.05; |
double zdelt = 0.00025; |
// 检查边界 |
bool has_bounds = !bounds.empty(); |
Vector lower_bound, upper_bound; |
if (has_bounds) { |
if (bounds.size() != static_cast<size_t>(x0.size())) { |
throw std::invalid_argument("Bounds size does not match x0 size"); |
} |
lower_bound(x0.size()); |
upper_bound(x0.size()); |
for (int i = 0; i < x0.size(); ++i) { |
lower_bound(i) = bounds[i].first; |
upper_bound(i) = bounds[i].second; |
} |
if ((lower_bound.array() > upper_bound.array()).any()) { |
throw std::invalid_argument("Nelder Mead - one of the lower bounds is greater than an upper bound."); |
} |
if (((x0.array() < lower_bound.array()) || (x0.array() > upper_bound.array())).any() && disp) { |
std::cerr << "Warning: Initial guess is not within the specified bounds" << std::endl; |
} |
} |
// 将点裁剪到边界内 |
auto clip = [&](const Vector& x) -> Vector { |
if (!has_bounds) |
return x; |
return x.cwiseMax(lower_bound).cwiseMin(upper_bound); |
}; |
// 按函数值排序单纯形顶点 |
auto sort_simplex = [](Vector& fsim, Matrix& sim) { |
const int N = sim.cols(); |
std::vector<std::pair<Scalar, int>> pairs(N + 1); |
for (int i = 0; i < N + 1; ++i) { |
pairs[i] = {fsim(i), i}; |
} |
std::sort(pairs.begin(), pairs.end()); |
// 重排fsim和sim |
Vector fsim_sorted(N + 1); |
Matrix sim_sorted(N + 1, N); |
for (int i = 0; i < N + 1; ++i) { |
fsim_sorted(i) = pairs[i].first; |
sim_sorted.row(i) = sim.row(pairs[i].second); |
} |
fsim = fsim_sorted; |
sim = sim_sorted; |
}; |
// 创建初始单纯形 |
Vector x0_copy = x0; |
if (has_bounds) { |
x0_copy = clip(x0_copy); |
} |
int N = x0_copy.size(); |
Matrix sim; |
if (initial_simplex == nullptr) { |
sim = Matrix(N + 1, N); |
sim.row(0) = x0_copy; |
for (int k = 0; k < N; ++k) { |
Vector y = x0_copy; |
if (y(k) != 0) { |
y(k) = (1.0 + nonzdelt) * y(k); |
} else { |
y(k) = zdelt; |
} |
sim.row(k + 1) = y; |
} |
} else { |
if (initial_simplex->rows() != N + 1 || initial_simplex->cols() != N) { |
throw std::invalid_argument("`initial_simplex` should be an array of shape (N+1,N)"); |
} |
sim = *initial_simplex; |
} |
// 如果边界存在，确保所有单纯形顶点都在边界内 |
if (has_bounds) { |
for (int i = 0; i < sim.rows(); ++i) { |
Vector vertex = sim.row(i); |
sim.row(i) = clip(vertex); |
} |
} |
// 如果既没有设置maxiter也没有设置maxfev，设置它们为默认值 |
if (maxiter < 0 && maxfev < 0) { |
maxiter = N * 200; |
maxfev = N * 200; |
} else if (maxiter < 0) { |
if (maxfev == std::numeric_limits<int>::max()) { |
maxiter = N * 200; |
} else { |
maxiter = std::numeric_limits<int>::max(); |
} |
} else if (maxfev < 0) { |
if (maxiter == std::numeric_limits<int>::max()) { |
maxfev = N * 200; |
} else { |
maxfev = std::numeric_limits<int>::max(); |
} |
} |
// 计算初始单纯形中每个点的函数值 |
Vector fsim = Vector::Constant(N + 1, std::numeric_limits<Scalar>::max()); |
int fcalls = 0; |
for (int k = 0; k < N + 1; ++k) { |
fsim(k) = func(sim.row(k).transpose()); |
fcalls++; |
if (fcalls >= maxfev) { |
break; |
} |
} |
// 按函数值排序单纯形顶点 |
sort_simplex(fsim, sim); |
// 保存所有迭代结果 |
std::vector<Vector> allvecs; |
if (return_all) { |
allvecs.push_back(sim.row(0).transpose()); |
} |
// 开始迭代 |
int iterations = 1; |
while (fcalls < maxfev && iterations < maxiter) { |
// 检查收敛 |
double max_dist = 0.0; |
for (int i = 1; i < N + 1; ++i) { |
max_dist = std::max(max_dist, (sim.row(i) - sim.row(0)).norm()); |
} |
double max_diff = 0.0; |
for (int i = 1; i < N + 1; ++i) { |
max_diff = std::max(max_diff, std::abs(fsim(i) - fsim(0))); |
} |
if (max_dist <= xatol && max_diff <= fatol) { |
break; |
} |
// 计算质心（除了最差点） |
// 1. sim.topRows(N) - 选择前N行，等价于Python的sim[:-1] |
// (因为单纯形有N+1个顶点，选择前N行就是除了最后一行外的所有行) |
// 2. .colwise().sum() - 对每一列分别求和，结果是一个行向量 |
// (等价于Python的np.add.reduce(..., 0)) |
// 3. .transpose() - 将行向量转置为列向量 |
// (因为在Eigen中Vector通常表示列向量) |
// 4. / static_cast<Scalar>(N) - 除以N得到平均值(质心) |
Vector xbar = sim.topRows(N).colwise().sum().transpose() / static_cast<Scalar>(N); |
// 执行反射 |
Vector xr = (1 + rho) * xbar - rho * sim.row(N).transpose(); |
if (has_bounds) { |
xr = clip(xr); |
} |
Scalar fxr = func(xr); |
fcalls++; |
bool doshrink = false; |
if (fxr < fsim(0)) { |
// 反射点比最好点还好，尝试扩展 |
Vector xe = (1 + rho * chi) * xbar - rho * chi * sim.row(N).transpose(); |
if (has_bounds) { |
xe = clip(xe); |
} |
Scalar fxe = func(xe); |
fcalls++; |
if (fxe < fxr) { |
// 扩展点更好，接受扩展 |
sim.row(N) = xe.transpose(); |
fsim(N) = fxe; |
} else { |
// 反射点更好，接受反射 |
更多推荐
6
0- 0
已为社区贡献4条内容
所有评论(0)