GSL库简介与安装

GNU科学计算库（GSL，GNU Scientific Library）是一个为C和C++程序员提供的数值例程集合。它涵盖了从基本数学函数到复杂的统计和优化算法的广泛领域，是进行科学计算的强大工具。

在基于Ubuntu的系统中安装GSL库非常简单，只需执行以下命令即可：

sudo apt-get install libgsl-dev

GSL多变量函数最小化器：gsl_multimin_fdfminimizer

在多变量函数优化中，共轭梯度法（Conjugate Gradient Method）因其高效性而广受欢迎。GSL库提供了gsl_multimin_fdfminimizer模块，实现了多种共轭梯度算法，用于寻找多变量函数的局部最小值。使用此模块进行优化的基本步骤如下：

1. 定义目标函数： 创建一个C函数，计算给定输入向量 \\(\vec{x}\\) 的函数值。此函数通常包含一个指向额外参数的void *指针。
2. 定义梯度函数： 创建一个C函数，计算给定输入向量 \\(\vec{x}\\) 的梯度向量 \\(\nabla f(\\vec{x})\\)。同样，也包含一个void *参数指针。
3. 整合函数与梯度： 将目标函数和梯度函数封装到一个gsl_multimin_function_fdf结构体中。这个结构体将作为优化器的输入。
4. 初始化并迭代优化器： 实例化一个gsl_multimin_fdfminimizer对象，设置初始点、步长和容差，然后进行迭代，直到梯度向量的范数小于预设阈值。

示例一：抛物面函数的最小化

我们首先通过一个简单的二维抛物面函数来演示gsl_multimin_fdfminimizer的使用。目标函数定义为：

\\(f(x, y) = A(x - x_0)^2 + B(y - y_0)^2 + C\\)

在这个例子中，我们设定 \\(x_0=1, y_0=2, A=10, B=20, C=30\\)，因此函数的理论最小值为30，在 \\((1, 2)\\) 处取得。

C++代码实现

#include <iostream>
#include <gsl/gsl_multimin.h>

// 全局计数器，用于统计目标函数被调用的次数
static int obj_func_call_count = 0;

/* 定义目标函数：一个中心在(p[0], p[1])，缩放因子为(p[2], p[3])，最小值为p[4]的抛物面 */
double paraboloid_objective(const gsl_vector *input_vec, void *params_ptr) {
    obj_func_call_count++;
    double x_val = gsl_vector_get(input_vec, 0);
    double y_val = gsl_vector_get(input_vec, 1);
    double *p_coeffs = (double *)params_ptr; // 获取参数数组

    return p_coeffs[2] * (x_val - p_coeffs[0]) * (x_val - p_coeffs[0]) +
           p_coeffs[3] * (y_val - p_coeffs[1]) * (y_val - p_coeffs[1]) +
           p_coeffs[4];
}

/* 定义梯度函数：计算目标函数在(x, y)处的梯度向量 */
void paraboloid_gradient(const gsl_vector *input_vec, void *params_ptr, gsl_vector *grad_vec) {
    double x_val = gsl_vector_get(input_vec, 0);
    double y_val = gsl_vector_get(input_vec, 1);
    double *p_coeffs = (double *)params_ptr;

    // df/dx
    gsl_vector_set(grad_vec, 0, 2.0 * p_coeffs[2] * (x_val - p_coeffs[0]));
    // df/dy
    gsl_vector_set(grad_vec, 1, 2.0 * p_coeffs[3] * (y_val - p_coeffs[1]));
}

/* 同时计算目标函数值和梯度向量 */
void paraboloid_obj_and_grad(const gsl_vector *input_vec, void *params_ptr, double *func_val, gsl_vector *grad_vec) {
    *func_val = paraboloid_objective(input_vec, params_ptr);
    paraboloid_gradient(input_vec, params_ptr, grad_vec);
}

int main() {
    // 设置目标函数的参数：中心(1,2)，缩放因子(10,20)，最小值30
    double function_params[5] = {1.0, 2.0, 10.0, 20.0, 30.0};

    // 初始化 gsl_multimin_function_fdf 结构体
    gsl_multimin_function_fdf my_obj_fdf;
    my_obj_fdf.n = 2; // 函数的变量数量
    my_obj_fdf.f = ¶boloid_objective;
    my_obj_fdf.df = ¶boloid_gradient;
    my_obj_fdf.fdf = ¶boloid_obj_and_grad;
    my_obj_fdf.params = (void *)function_params; // 传递函数参数

    // 设置初始迭代点 (x, y)
    gsl_vector *initial_point = gsl_vector_alloc(2);
    gsl_vector_set(initial_point, 0, 5.0);
    gsl_vector_set(initial_point, 1, 7.0);

    // 选择共轭梯度法类型，这里使用 Fletcher-Reeves 算法
    const gsl_multimin_fdfminimizer_type *minimizer_type = gsl_multimin_fdfminimizer_conjugate_fr;

    // 分配优化器空间
    gsl_multimin_fdfminimizer *optimizer_state = gsl_multimin_fdfminimizer_alloc(minimizer_type, 2);
    // 设置优化器：函数、初始点、初始步长、一维搜索精度
    gsl_multimin_fdfminimizer_set(optimizer_state, &my_obj_fdf, initial_point, 0.01, 1E-4);

    int iteration = 0;
    int status_code;
    do {
        iteration++;
        status_code = gsl_multimin_fdfminimizer_iterate(optimizer_state); // 执行一次迭代

        if (status_code) { // 如果迭代出错
            break;
        }

        // 检查梯度范数是否小于容差，判断是否收敛
        status_code = gsl_multimin_test_gradient(optimizer_state->gradient, 1E-3);

        if (status_code == GSL_SUCCESS) {
            printf("最小点已找到：\n");
        }
        printf("%5d %.10f %.10f %10.5f\n", iteration,
               gsl_vector_get(optimizer_state->x, 0),
               gsl_vector_get(optimizer_state->x, 1),
               optimizer_state->f);

    } while (status_code == GSL_CONTINUE && iteration < 100); // 最多迭代100次

    std::cout << "目标函数调用次数: " << obj_func_call_count << std::endl;

    // 释放GSL分配的内存
    gsl_multimin_fdfminimizer_free(optimizer_state);
    gsl_vector_free(initial_point);

    return 0;
}

编译、运行与结果

使用GCC编译上述代码，并链接GSL和BLAS库：

g++ paraboloid_minimize.cpp -o paraboloid_minimize -lgsl -lblas
./paraboloid_minimize

运行结果示例：

    1 4.9962860932 6.9907152331 687.84780
    2 4.9888582797 6.9721456993 683.55456
    ...
最小点已找到：
   13 1.0000000000 2.0000000000   30.00000
目标函数调用次数: 28

从输出可以看出，优化器经过13次迭代找到了函数的最小值30，对应的位置是 \\((1.0, 2.0)\\)，与理论值完全一致。目标函数被调用了28次。

结果可视化

为了直观地展示优化过程，我们可以将迭代路径绘制在函数的等高线上。保存上述输出的迭代点数据，例如到optimization_path.txt文件中，然后使用Python脚本进行可视化：

import numpy as np
import matplotlib.pyplot as plt

# 从文件中加载迭代路径数据
path_data = np.loadtxt('optimization_path.txt', usecols=(1, 2)) # 提取x和y坐标

# 绘制迭代点
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(path_data[:, 0], path_data[:, 1], color='red', marker='o', label='迭代路径')
ax.plot(path_data[:, 0], path_data[:, 1], color='red', linestyle='--', linewidth=0.8)

# 生成等高线数据
step_size = 0.05
x_range = np.arange(-5, 10, step_size)
y_range = np.arange(-5, 10, step_size)
X, Y = np.meshgrid(x_range, y_range)

# 定义抛物面函数 (注意与C++代码中的参数保持一致)
# f(x, y) = 10*(x-1)^2 + 20*(y-2)^2 + 30
Z = 10 * (X - 1)**2 + 20 * (Y - 2)**2 + 30

# 绘制等高线
ax.contour(X, Y, Z, levels=np.logspace(0, 3, 20), cmap='viridis') # 使用对数间隔的等高线
ax.set_xlabel("X轴")
ax.set_ylabel("Y轴")
ax.set_title("抛物面函数最小化迭代路径")
ax.legend()
ax.grid(True)
plt.show()

![image]()

图中红点代表了优化器在每次迭代中探索到的位置，红线则连接了这些点，形成了收敛路径。可以看到，路径有效地沿着梯度下降方向向函数最小值收敛。

示例二：多变量维度下的性能表现

接下来，我们测试当函数变量数量增加时，共轭梯度法所需的迭代次数。我们将抛物面函数推广到 \\(N\\) 维：

\\(f(\\vec{x}) = \sum_{i=0}^{N-1} p_{N+i} (x_i - p_i)^2 + p_{2N}\\)

其中 \\(p_i\\) 是每个维度的中心值，\\(p_{N+i}\\) 是缩放因子，\\(p_{2N}\\) 是常数项。我们将 \\(N\\) 从50增加到1000，观察收敛所需的迭代次数。

C++代码实现

#include <iostream>
#include <iomanip>
#include <fstream>
#include <gsl/gsl_multimin.h>
#include <vector> // 使用vector代替C风格数组

static int current_dim_n = 0; // 全局变量，表示当前函数维度

/* N维抛物面目标函数 */
double n_dim_paraboloid_obj(const gsl_vector *v, void *params) {
    double *p_data = (double *)params;
    double x_coord, func_sum = 0.0;
    for (int i = 0; i < current_dim_n; i++) {
        x_coord = gsl_vector_get(v, i);
        func_sum += p_data[current_dim_n + i] * (x_coord - p_data[i]) * (x_coord - p_data[i]);
    }
    func_sum += p_data[2 * current_dim_n];
    return func_sum;
}

/* N维抛物面梯度函数 */
void n_dim_paraboloid_grad(const gsl_vector *v, void *params, gsl_vector *df_out) {
    double *p_data = (double *)params;
    double x_coord;
    for (int i = 0; i < current_dim_n; i++) {
        x_coord = gsl_vector_get(v, i);
        gsl_vector_set(df_out, i, 2.0 * p_data[current_dim_n + i] * (x_coord - p_data[i]));
    }
}

/* N维抛物面函数与梯度联合计算 */
void n_dim_paraboloid_obj_grad(const gsl_vector *x, void *params, double *f_out, gsl_vector *df_out) {
    *f_out = n_dim_paraboloid_obj(x, params);
    n_dim_paraboloid_grad(x, params, df_out);
}

// 封装最小化过程的函数
void execute_dim_test(std::ofstream &output_file) {
    gsl_multimin_function_fdf test_func_fdf;
    
    // 动态分配参数数组
    std::vector<double> p_coeffs(2 * current_dim_n + 1);
    for (int i = 0; i < current_dim_n; i++) {
        p_coeffs[i] = i + 1.0;         // 中心点坐标
        p_coeffs[current_dim_n + i] = 10.0 * i + 10.0; // 缩放因子
    }
    p_coeffs[2 * current_dim_n] = 100.0; // 常数项

    test_func_fdf.n = current_dim_n;
    test_func_fdf.f = &n_dim_paraboloid_obj;
    test_func_fdf.df = &n_dim_paraboloid_grad;
    test_func_fdf.fdf = &n_dim_paraboloid_obj_grad;
    test_func_fdf.params = (void *)p_coeffs.data(); // 使用vector的底层数据指针

    gsl_vector *initial_guess_vec = gsl_vector_alloc(current_dim_n);
    for (int i = 0; i < current_dim_n; i++) {
        gsl_vector_set(initial_guess_vec, i, 10.0); // 初始点设置为所有维度都为10.0
    }

    const gsl_multimin_fdfminimizer_type *type = gsl_multimin_fdfminimizer_conjugate_fr;
    gsl_multimin_fdfminimizer *minimizer_state = gsl_multimin_fdfminimizer_alloc(type, current_dim_n);
    gsl_multimin_fdfminimizer_set(minimizer_state, &test_func_fdf, initial_guess_vec, 0.01, 1E-4);

    int iter_count = 0;
    int status_val;
    do {
        iter_count++;
        status_val = gsl_multimin_fdfminimizer_iterate(minimizer_state);
        if (status_val) break; // 迭代出错则退出

        status_val = gsl_multimin_test_gradient(minimizer_state->gradient, 1E-3); // 检查梯度收敛

        if (status_val == GSL_SUCCESS) {
            output_file << current_dim_n << "\t" << iter_count << std::endl;
            printf("维度 n=%d, 经过 %d 次迭代找到最小值\n", current_dim_n, iter_count);
        }
    } while (status_val == GSL_CONTINUE && iter_count < 10000); // 最大迭代10000次

    gsl_multimin_fdfminimizer_free(minimizer_state);
    gsl_vector_free(initial_guess_vec);
}

int main() {
    std::ofstream data_out("iteration_scaling.txt");

    for (current_dim_n = 50; current_dim_n <= 1000; current_dim_n += 50) {
        execute_dim_test(data_out);
    }

    data_out.close();
    return 0;
}

结果可视化

编译并运行上述代码，会将每个维度数对应的迭代次数输出到iteration_scaling.txt文件中。然后使用Python脚本绘制散点图：

import numpy as np
import matplotlib.pyplot as plt

# 从文件中加载数据
data = np.loadtxt('iteration_scaling.txt')
dimensions = data[:, 0]
iterations = data[:, 1]

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(dimensions, iterations, color='blue', alpha=0.7)
ax.set_xlim(0, 1050)
ax.set_ylim(0, 750)
ax.set_xlabel("多元函数的自由参数数量 (N)")
ax.set_ylabel("梯度范数收敛所需迭代次数 (|g| < 1E-3)")
ax.set_title("共轭梯度法迭代次数随维度N的变化")
ax.grid(True)
plt.show()

![image]()

从图中可以看出，对于这个特定的抛物面函数，共轭梯度法所需的迭代次数随着变量维度的增加而近似呈线性增长。这表明该算法在处理高维凸函数时仍然保持了较好的性能。

示例三：带软约束的优化问题（通过惩罚函数）

在实际应用中，我们经常需要对优化变量施加约束。虽然GSL的gsl_multimin_fdfminimizer主要用于无约束优化，但我们可以通过引入惩罚函数（Penalty Function）将约束问题转化为无约束问题。这里我们尝试优化函数：

\\(f(\\vec{r}) = \\frac{10(x^2 + 2y^2)}{x^2 + y^2}\\)

该函数在 \\(y=0\\) 时取得最小值10，但对 \\(x\\) 的大小没有限制。现在我们希望找到一个最小值点，同时要求向量的模长 \\(|\\vec{r}|\\) 接近1。为此，我们引入一个惩罚项：

\\(g(\\vec{r}) = 100 \cdot (||\\vec{r}|| - 1)^2\\)，当 \\(||\\vec{r}|| > 1\\) 或 \\(||\\vec{r}|| < 0.5\\) 时；否则为0。

最终优化的函数是 \\(f'(\\vec{r}) = f(\\vec{r}) + g(\\vec{r})\\)。

C++代码实现

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath> // for sqrt
#include <gsl/gsl_multimin.h>

/* 定义带惩罚项的目标函数 */
double penalized_objective(const gsl_vector *v, void *params) {
    double x = gsl_vector_get(v, 0);
    double y = gsl_vector_get(v, 1);
    double r_squared = x * x + y * y; // 向量模长的平方
    double r_norm = sqrt(r_squared); // 向量模长

    double base_val = 10.0 * (x * x + 2.0 * y * y) / r_squared; // 原始目标函数

    // 添加惩罚项
    if (r_norm > 1.0 || r_norm < 0.5) {
        return base_val + 100.0 * (r_norm - 1.0) * (r_norm - 1.0);
    } else {
        return base_val;
    }
}

/* 定义带惩罚项的梯度函数 */
void penalized_gradient(const gsl_vector *v, void *params, gsl_vector *df_out) {
    double x = gsl_vector_get(v, 0);
    double y = gsl_vector_get(v, 1);
    double r_squared = x * x + y * y;
    double r_norm = sqrt(r_squared);

    // 计算原始函数的梯度分量
    // df/dx = (20xy^2) / (x^2+y^2)^2 - (20x) / (x^2+y^2)
    // 简化后: df/dx = -20 * x * y^2 / (x^2+y^2)^2
    // df/dy = (20x^2y) / (x^2+y^2)^2
    double df_dx_base = -20.0 * x * y * y / (r_squared * r_squared);
    double df_dy_base = 20.0 * x * x * y / (r_squared * r_squared);

    // 添加惩罚项的梯度分量
    double penalty_grad_x = 0.0;
    double penalty_grad_y = 0.0;
    if (r_norm > 1.0 || r_norm < 0.5) {
        // g(r) = 100 * (r_norm - 1)^2
        // dg/dx = dg/dr_norm * dr_norm/dx
        // dg/dr_norm = 200 * (r_norm - 1)
        // dr_norm/dx = x / r_norm
        // dr_norm/dy = y / r_norm
        penalty_grad_x = 200.0 * (r_norm - 1.0) * (x / r_norm);
        penalty_grad_y = 200.0 * (r_norm - 1.0) * (y / r_norm);
    }

    gsl_vector_set(df_out, 0, df_dx_base + penalty_grad_x);
    gsl_vector_set(df_out, 1, df_dy_base + penalty_grad_y);
}

/* 同时计算带惩罚项的目标函数值和梯度 */
void penalized_obj_grad_combined(const gsl_vector *x, void *params, double *f_out, gsl_vector *df_out) {
    *f_out = penalized_objective(x, params);
    penalized_gradient(x, params, df_out);
}

/* 通用的GSL FDF最小化器封装 */
void generic_gsl_minimizer(int num_vars, double *initial_coords, double init_step_size, double line_search_tol, double grad_eps_abs, int max_iterations,
                           double (*obj_func)(const gsl_vector *, void *),
                           void (*grad_func)(const gsl_vector *, void *, gsl_vector *),
                           void (*obj_grad_func)(const gsl_vector *, void *, double *, gsl_vector *)) {

    gsl_multimin_function_fdf target_func_fdf;
    target_func_fdf.n = num_vars;
    target_func_fdf.f = obj_func;
    target_func_fdf.df = grad_func;
    target_func_fdf.fdf = obj_grad_func;
    target_func_fdf.params = nullptr; // 本例中不需要额外参数

    gsl_vector *current_x = gsl_vector_alloc(num_vars);
    for (int i = 0; i < num_vars; i++) {
        gsl_vector_set(current_x, i, initial_coords[i]);
    }

    const gsl_multimin_fdfminimizer_type *type = gsl_multimin_fdfminimizer_conjugate_fr;
    gsl_multimin_fdfminimizer *minimizer_inst = gsl_multimin_fdfminimizer_alloc(type, num_vars);
    gsl_multimin_fdfminimizer_set(minimizer_inst, &target_func_fdf, current_x, init_step_size, line_search_tol);

    int iter_count = 0;
    int status_code;
    do {
        iter_count++;
        status_code = gsl_multimin_fdfminimizer_iterate(minimizer_inst);
        if (status_code) break;

        status_code = gsl_multimin_test_gradient(minimizer_inst->gradient, grad_eps_abs);

        if (status_code == GSL_SUCCESS) {
            printf("在 %d 次迭代后找到最小值！\n", iter_count);
        }
        std::cout << "\t " << iter_count;
        for (int j = 0; j < num_vars; j++) {
            std::cout << "\t " << gsl_vector_get(minimizer_inst->x, j);
        }
        std::cout << "\t " << minimizer_inst->f << std::endl;
    } while (status_code == GSL_CONTINUE && iter_count < max_iterations);

    if (iter_count == max_iterations) {
        std::cout << "未能收敛，已达到最大迭代次数 " << max_iterations << std::endl;
    }

    gsl_multimin_fdfminimizer_free(minimizer_inst);
    gsl_vector_free(current_x);
}

int main() {
    int variables_count = 2; // 二维变量
    double initial_values[] = {5.0, 7.0}; // 初始位置
    double initial_step = 0.1;
    double line_search_precision = 0.001; // 一维搜索精度
    double gradient_tolerance = 1E-3;     // 梯度范数收敛阈值
    int maximum_iterations = 10000;       // 最大迭代次数

    generic_gsl_minimizer(variables_count, initial_values, initial_step, line_search_precision, gradient_tolerance, maximum_iterations,
                          penalized_objective, penalized_gradient, penalized_obj_grad_combined);

    return 0;
}

运行结果

编译并运行上述代码，输出结果示例如下：

	 1	 4.94188	 6.91863	 508236
	 2	 4.82563	 6.75588	 461446
	 ...
	 11	 0.694989	 -0.131158	 10.3439
	 12	 0.702905	 -0.0892135	 10.1585
	 13	 0.718737	 -0.00532433	 10.0005
在 14 次迭代后找到最小值！
	 14	 0.719741	 -3.544e-11	 10

从结果可以看出，优化器成功收敛。最终找到的最小点近似为 \\((0.719741, 0.0)\\)，函数值为10。该点位于 \\(y=0\\) 轴上，满足原始函数的最小值条件。同时，向量的模长约为 \\(\\sqrt{0.719741^2 + 0^2} \\approx 0.719741\\)，这个值介于0.5和1之间，说明惩罚函数有效地引导优化器找到了满足约束的区域内的最小值点。