C++代码实现MATLAB中的fminunc函数功能

// =====================================================================
// fminunc.cpp
// 纯 C++ 实现的 BFGS 拟牛顿优化器，对应 MATLAB fminunc 的核心算法
// 仅依赖 C++ 标准库，可独立编译运行
// =====================================================================
#include <iostream>
#include <vector>
#include <cmath>
#include <functional>
#include <iomanip>
#include <algorithm>
#include <string>
#include <stdexcept>

using Vec = std::vector<double>;
using ObjFunc = std::function<double(const Vec&, Vec&)>;

// —————- 基本向量运算 —————-
static double dot(const Vec& a, const Vec& b) {
double s = 0.0;
for (size_t i = 0; i < a.size(); ++i) s += a[i] * b[i];
return s;
}
static double norm2(const Vec& a) { return std::sqrt(dot(a, a)); }

// —————- 线搜索：强 Wolfe 条件 —————-
// 前向声明
static double zoom_search(const ObjFunc& func,
const Vec& x, double f0, const Vec& g0,
const Vec& p, double gp0,
double alpha_lo, double alpha_hi);

// 主搜索：先用扩展步长找区间，再交给 zoom
static double line_search(const ObjFunc& func,
const Vec& x, double f0, const Vec& g0,
const Vec& p, double gp0) {
const double c1 = 1e-4; // Armijo 常数
const double c2 = 0.9; // 曲率常数
const double alpha_max = 10.0;
double alpha_prev = 0.0;
double alpha = 1.0;
double f_prev = f0;

Vec g_tmp(x.size());

for (int i = 0; i < 50; ++i) {
Vec x_try(x.size());
for (size_t j = 0; j < x.size(); ++j)
x_try[j] = x[j] + alpha * p[j];
double f_try = func(x_try, g_tmp);
double gp_try = dot(g_tmp, p);

// Armijo 不满足，或函数值反而上升
if (f_try > f0 + c1 * alpha * gp0 || (i > 0 && f_try >= f_prev)) {
return zoom_search(func, x, f0, g0, p, gp0, alpha_prev, alpha);
}
// 强 Wolfe 曲率条件
if (std::fabs(gp_try) <= –c2 * gp0) {
return alpha;
}
// 方向导数变号
if (gp_try >= 0.0) {
return zoom_search(func, x, f0, g0, p, gp0, alpha, alpha_prev);
}
alpha_prev = alpha;
f_prev = f_try;
alpha = std::min(alpha * 2.0, alpha_max);
}
return alpha;
}

// 在 [alpha_lo, alpha_hi] 内用二分法缩小区间
static double zoom_search(const ObjFunc& func,
const Vec& x, double f0, const Vec& g0,
const Vec& p, double gp0,
double alpha_lo, double alpha_hi) {
const double c1 = 1e-4;
const double c2 = 0.9;
Vec g_tmp(x.size());

double f_lo;
{
Vec x_lo(x.size());
for (size_t j = 0; j < x.size(); ++j)
x_lo[j] = x[j] + alpha_lo * p[j];
f_lo = func(x_lo, g_tmp);
}

for (int i = 0; i < 60; ++i) {
double alpha = 0.5 * (alpha_lo + alpha_hi);
Vec x_try(x.size());
for (size_t j = 0; j < x.size(); ++j)
x_try[j] = x[j] + alpha * p[j];
double f_try = func(x_try, g_tmp);
double gp_try = dot(g_tmp, p);

if (f_try > f0 + c1 * alpha * gp0 || f_try >= f_lo) {
alpha_hi = alpha;
} else {
if (std::fabs(gp_try) <= –c2 * gp0) return alpha;
if (gp_try * (alpha_hi – alpha_lo) >= 0.0) alpha_hi = alpha_lo;
alpha_lo = alpha;
f_lo = f_try;
}
if (std::fabs(alpha_hi – alpha_lo) < 1e-12 * (1.0 + std::fabs(alpha_lo)))
break;
}
return alpha_lo;
}

// —————- BFGS 主优化器 —————-
struct OptResult {
Vec x; // 最优解
double fval; // 最优函数值
int iterations; // 迭代次数
int fevals; // 函数评估次数
bool converged; // 是否收敛
std::string message; // 结束信息
};

static OptResult fminunc_bfgs(const ObjFunc& func,
const Vec& x0,
double gtol = 1e-6,
double xtol = 1e-10,
int maxiter = 1000,
bool verbose = true) {
const int n = static_cast<int>(x0.size());
Vec x = x0;
Vec g(n, 0.0);
double f = func(x, g);
int fevals = 1;

// 逆 Hessian 近似 H，初始为单位矩阵（按行主序存储 n×n）
Vec H(n * n, 0.0);
for (int i = 0; i < n; ++i) H[i * n + i] = 1.0;

if (verbose) {
std::cout << std::scientific << std::setprecision(6);
std::cout << "Iter Fcount f(x) |g|\\n";
std::cout << std::string(56, '-') << "\\n";
std::cout << std::setw(4) << 0
<< std::setw(10) << fevals
<< std::setw(18) << f
<< std::setw(16) << norm2(g) << "\\n";
}

for (int iter = 0; iter < maxiter; ++iter) {
// 梯度收敛判据
double gnorm = norm2(g);
if (gnorm < gtol) {
return {x, f, iter, fevals, true, "梯度范数小于容差"};
}

// 搜索方向 p = -H * g
Vec p(n, 0.0);
for (int i = 0; i < n; ++i) {
double s = 0.0;
for (int j = 0; j < n; ++j) s += H[i * n + j] * g[j];
p[i] = –s;
}

// 若方向不下降，重置 H 为单位阵
double gp = dot(g, p);
if (gp >= 0.0) {
std::fill(H.begin(), H.end(), 0.0);
for (int i = 0; i < n; ++i) H[i * n + i] = 1.0;
for (int i = 0; i < n; ++i) p[i] = –g[i];
gp = dot(g, p);
}

// 线搜索
double alpha = line_search(func, x, f, g, p, gp);

// 试探新点
Vec x_new(n), g_new(n, 0.0);
for (int i = 0; i < n; ++i) x_new[i] = x[i] + alpha * p[i];
double f_new = func(x_new, g_new);
++fevals;

// 计算 s = x_new – x, y = g_new – g
Vec s(n), y(n);
for (int i = 0; i < n; ++i) {
s[i] = x_new[i] – x[i];
y[i] = g_new[i] – g[i];
}

// 曲率条件满足则做 BFGS 更新，否则重置
double sy = dot(s, y);
if (sy > 1e-10 * norm2(s) * norm2(y)) {
double rho = 1.0 / sy;
// Hy
Vec Hy(n, 0.0);
for (int i = 0; i < n; ++i) {
double acc = 0.0;
for (int j = 0; j < n; ++j) acc += H[i * n + j] * y[j];
Hy[i] = acc;
}
double yHy = dot(y, Hy);
double coef = rho * rho * yHy + rho;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
H[i * n + j] += coef * s[i] * s[j]
– rho * (Hy[i] * s[j] + s[i] * Hy[j]);
}
}
} else {
std::fill(H.begin(), H.end(), 0.0);
for (int i = 0; i < n; ++i) H[i * n + i] = 1.0;
}

double xdiff = norm2(s);
x = x_new;
g = g_new;
f = f_new;

if (verbose) {
std::cout << std::setw(4) << iter + 1
<< std::setw(10) << fevals
<< std::setw(18) << f
<< std::setw(16) << norm2(g) << "\\n";
}

if (xdiff < xtol * (1.0 + norm2(x))) {
return {x, f, iter + 1, fevals, true, "x 变化小于容差"};
}
}

return {x, f, maxiter, fevals, false, "达到最大迭代次数"};
}

// —————- 主程序：测试 —————-
int main() {
// ———- 测试 1: Rosenbrock 函数 ———-
std::cout << "===== 测试 1: Rosenbrock 函数 =====\\n";
auto rosen = [](const Vec& x, Vec& g) -> double {
double x0 = x[0], x1 = x[1];
double f = 100.0 * std::pow(x1 – x0 * x0, 2) + std::pow(1.0 – x0, 2);
if (!g.empty()) {
g[0] = –400.0 * x0 * (x1 – x0 * x0) – 2.0 * (1.0 – x0);
g[1] = 200.0 * (x1 – x0 * x0);
}
return f;
};

Vec x0 = {–1.2, 1.0};
OptResult res = fminunc_bfgs(rosen, x0, 1e-8, 1e-12, 1000, true);
std::cout << "\\n收敛: " << (res.converged ? "是" : "否")
<< " (" << res.message << ")\\n";
std::cout << "迭代次数: " << res.iterations << "\\n";
std::cout << "函数评估: " << res.fevals << "\\n";
std::cout << std::fixed << std::setprecision(10);
std::cout << "最优解: x = [" << res.x[0] << ", " << res.x[1] << "]\\n";
std::cout << "最优值: f = " << res.fval << "\\n";

// ———- 测试 2: 二次函数 ———-
std::cout << "\\n===== 测试 2: 二次函数 =====\\n";
auto quad = [](const Vec& x, Vec& g) -> double {
double x0 = x[0], x1 = x[1];
double f = (x0 – 2.0) * (x0 – 2.0) + (x1 – 3.0) * (x1 – 3.0) + 1.0;
if (!g.empty()) {
g[0] = 2.0 * (x0 – 2.0);
g[1] = 2.0 * (x1 – 3.0);
}
return f;
};
Vec x1 = {0.0, 0.0};
OptResult res2 = fminunc_bfgs(quad, x1, 1e-8, 1e-12, 200, true);
std::cout << "\\n收敛: " << (res2.converged ? "是" : "否")
<< " (" << res2.message << ")\\n";
std::cout << "迭代次数: " << res2.iterations << "\\n";
std::cout << "函数评估: " << res2.fevals << "\\n";
std::cout << "最优解: x = [" << res2.x[0] << ", " << res2.x[1] << "]\\n";
std::cout << "最优值: f = " << res2.fval << "\\n";

return 0;
}