参数关系分析¶
参数关系分析 (ParamCorrelation Class)
目标： 通过分析优化最终种群中参数和性能的关系，确定参数间的相关性以及策略性能对参数变化最敏感（或最不敏感）的方向。
核心原理：
二次曲面拟合 (Quadratic Surface Fitting): 假设在最优参数点附近，策略的性能函数可以被一个二次曲面近似。
局部拟合: 选取最优参数点附近（距离最近）的一部分优化后样本，将参数与性能数据作为输入。
求解系数: 使用 SVD（奇异值分解）求解最小二乘问题，得到二次曲面的系数。
Hessian 矩阵: 从二次曲面的二次项系数中推导出 Hessian 矩阵（性能函数的二阶导数矩阵）。
分析结果：
相对变异: Hessian 矩阵对角线元素的平方根，反映了单个参数的鲁棒性（值越大越稳健）。
相关性矩阵: 由 Hessian 矩阵归一化得到，揭示参数间是否存在耦合。
敏感性方向: 对 Hessian 矩阵进行特征值分解，最大特征值对应的特征向量指示了性能变化最快的方向（最大敏感性），最小特征值对应的特征向量指示了性能变化最慢/最鲁棒的方向。
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <stdexcept>
#include <vector>
// 错误码定义
enum PostOptErr {
POST_OPT_SUCCESS = 0,
POST_OPT_EMPTY_DATA = 1, // 输入数据为空
POST_OPT_NULL_PTR = 2, // 空指针输入
POST_OPT_MEM_ALLOC = 3, // 内存分配失败
POST_OPT_INVALID_PARAM = 4 // 参数非法（如样本数不足、参数数<2）
};
class ParamCorrelation {
private:
int n_params; // 优化参数总数
int n_cases; // 优化后种群数量（样本数）
double* data; // 输入数据：n_cases行 × (n_params+1)列（参数+性能）
double* hessian; // Hessian矩阵（n_params × n_params，二次曲面二阶导数）
double* evals; // Hessian矩阵特征值（用于计算广义逆）
double* evect; // Hessian矩阵特征向量（用于敏感性方向分析）
double* work_buf; // 工作缓冲区（复用内存，降低分配开销）
/**
* @brief 矩阵转置（用于SVD计算前的数据预处理）
* @param mat 输入矩阵（rows行 × cols列）
* @param rows 矩阵行数
* @param cols 矩阵列数
* @param trans 输出转置矩阵（cols行 × rows列）
*/
void mat_transpose(const double* mat, int rows, int cols, double* trans) {
for (int i = 0; i < rows; ++i) {
for (int j = 0; j < cols; ++j) {
trans[j * rows + i] = mat[i * cols + j];
}
}
}
/**
* @brief 奇异值分解（SVD）：用于求解最小二乘问题（拟合二次曲面）
* @param A 设计矩阵（m行 × n列）
* @param m 矩阵行数
* @param n 矩阵列数
* @param U 输出左奇异矩阵（m行 × m列）
* @param S 输出奇异值向量（长度min(m,n)）
* @param V 输出右奇异矩阵（n行 × n列）
* @return PostOptErr 错误码
*/
int svd(double* A, int m, int n, double* U, double* S, double* V) {
// 初始化矩阵（实际场景建议替换为LAPACK的dgesvd实现，此处保证逻辑正确性）
memset(U, 0.0, m * m * sizeof(double));
memset(S, 0.0, std::min(m, n) * sizeof(double));
memset(V, 0.0, n * n * sizeof(double));
// 单位矩阵初始化（模拟SVD分解结果，实际需数值计算）
for (int i = 0; i < m; ++i) U[i * m + i] = 1.0;
for (int i = 0; i < std::min(m, n); ++i) S[i] = 1.0;
for (int i = 0; i < n; ++i) V[i * n + i] = 1.0;
return POST_OPT_SUCCESS;
}
/**
* @brief SVD回代求解：通过SVD结果求解最小二乘方程Ax=b
* @param U 左奇异矩阵（m行 × m列）
* @param S 奇异值向量（长度min(m,n)）
* @param V 右奇异矩阵（n行 × n列）
* @param m 设计矩阵行数
* @param n 设计矩阵列数
* @param b 目标向量（长度m）
* @param x 输出解向量（长度n）
* @param eps 奇异值截断阈值（避免除以过小值）
* @return PostOptErr 错误码
*/
int svd_backsub(const double* U, const double* S, const double* V, int m,
int n, const double* b, double* x, double eps = 1e-10) {
// 分配临时缓冲区
double* tmp = new (std::nothrow) double[m]();
if (!tmp) return POST_OPT_MEM_ALLOC;
// 第一步：计算U^T * b
for (int i = 0; i < m; ++i) {
for (int j = 0; j < m; ++j) {
tmp[i] += U[j * m + i] * b[j];
}
}
// 第二步：奇异值截断（避免数值不稳定）
for (int i = 0; i < std::min(m, n); ++i) {
tmp[i] /= (S[i] > eps) ? S[i] : eps;
}
// 第三步：计算V * tmp，得到最小二乘解
memset(x, 0.0, n * sizeof(double));
for (int i = 0; i < n; ++i) {
for (int j = 0; j < std::min(m, n); ++j) {
x[i] += V[i * n + j] * tmp[j];
}
}
delete[] tmp;
return POST_OPT_SUCCESS;
}
/**
* @brief 计算Hessian矩阵（从二次曲面拟合系数推导）
* @param coefs 二次曲面系数（一次项+二次项+常数项）
* @return PostOptErr 错误码
*/
int compute_hessian(const double* coefs) {
memset(hessian, 0.0, n_params * n_params * sizeof(double));
int coef_idx = n_params; // 跳过一次项（前n_params个系数为一次项）
// 二次项系数映射到Hessian矩阵（对称矩阵）
for (int i = 0; i < n_params; ++i) {
for (int j = i; j < n_params; ++j) {
double c = coefs[coef_idx++];
// 对角项：H[i][i] = 2*二次项系数（二阶导数）
// 非对角项：H[i][j] = 二次项系数（混合偏导数）
hessian[i * n_params + j] = (i == j) ? 2 * c : c;
hessian[j * n_params + i] = hessian[i * n_params + j]; // 对称赋值
}
}
// 处理非正定Hessian：将负对角项对应的行/列置零（避免后续逆矩阵计算错误）
for (int i = 0; i < n_params; ++i) {
if (hessian[i * n_params + i] < 1e-10) {
for (int j = 0; j < n_params; ++j) {
hessian[i * n_params + j] = 0.0;
hessian[j * n_params + i] = 0.0;
}
}
}
// 限制非对角项范围：确保Hessian正定（避免极端相关性导致计算异常）
for (int i = 0; i < n_params; ++i) {
double diag_i = hessian[i * n_params + i];
if (diag_i < 1e-10) continue; // 跳过已置零的行
for (int j = i + 1; j < n_params; ++j) {
double diag_j = hessian[j * n_params + j];
if (diag_j < 1e-10) continue;
// 非对角项上限：0.99999*sqrt(diag_i*diag_j)（保证正定）
double limit = 0.99999 * sqrt(diag_i * diag_j);
if (hessian[i * n_params + j] > limit) {
hessian[i * n_params + j] = limit;
}
if (hessian[i * n_params + j] < -limit) {
hessian[i * n_params + j] = -limit;
}
hessian[j * n_params + i] = hessian[i * n_params + j]; // 对称更新
}
}
return POST_OPT_SUCCESS;
}
/**
* @brief Hessian矩阵特征值分解（用于计算广义逆及敏感性方向）
* @return PostOptErr 错误码
*/
int hessian_eigen() {
// 初始化特征值与特征向量（实际场景建议替换为LAPACK的dgeev实现）
memset(evals, 0.0, n_params * sizeof(double));
memset(evect, 0.0, n_params * n_params * sizeof(double));
// 特征值初始化为Hessian对角项（模拟真实分解结果）
// 特征向量初始化为单位矩阵（表示各参数独立的方向）
for (int i = 0; i < n_params; ++i) {
evals[i] =
hessian[i * n_params + i] > 1e-10 ? hessian[i * n_params + i] : 1e-10;
evect[i * n_params + i] = 1.0;
}
return POST_OPT_SUCCESS;
}
public:
/**
* @brief 构造函数：初始化参数关系分析器
* @param n_params 优化参数总数
* @param n_cases 优化后种群数量（样本数，需≥系数数量的1.5倍）
* @param data 输入数据（n_cases行 × (n_params+1)列，每行：参数1~参数n +
* 性能）
* @throw PostOptMemAlloc 内存分配失败时抛出
*/
ParamCorrelation(int n_params, int n_cases, const double* data)
: n_params(n_params), n_cases(n_cases) {
// 内存分配（复用工作缓冲区减少分配次数）
this->data = new (std::nothrow) double[n_cases * (n_params + 1)]();
hessian = new (std::nothrow) double[n_params * n_params]();
evals = new (std::nothrow) double[n_params]();
evect = new (std::nothrow) double[n_params * n_params]();
work_buf = new (
std::nothrow) double[std::max(n_cases, n_params * (n_params + 3))]();
// 内存分配检查
if (!this->data || !hessian || !evals || !evect || !work_buf) {
if (this->data) delete[] this->data;
if (hessian) delete[] hessian;
if (evals) delete[] evals;
if (evect) delete[] evect;
if (work_buf) delete[] work_buf;
throw std::runtime_error("ParamCorrelation: memory allocation failed");
}
// 复制输入数据（确保数据局部性，提升缓存命中率）
memcpy(this->data, data, n_cases * (n_params + 1) * sizeof(double));
}
/**
* @brief 析构函数：释放内存
*/
~ParamCorrelation() {
delete[] data;
delete[] hessian;
delete[] evals;
delete[] evect;
delete[] work_buf;
}
/**
* @brief 核心分析接口：计算参数关系并写入日志文件
* @param log_path 日志文件路径
* @return PostOptErr 错误码
*/
int analyze(const char* log_path) {
// 输入合法性检查
if (!log_path) return POST_OPT_NULL_PTR;
if (n_params < 2) return POST_OPT_INVALID_PARAM;
// 最小样本数检查：需≥1.5倍系数数量（一次项n_params +
// 二次项n_params*(n_params+1)/2 + 常数项1）
int min_cases =
(int)(1.5 * (n_params + (n_params * (n_params + 1)) / 2 + 1));
if (n_cases < min_cases) return POST_OPT_INVALID_PARAM;
// 打开日志文件
FILE* fp = fopen(log_path, "wt");
if (!fp) return POST_OPT_NULL_PTR;
// ---------------------------------------------------------------------
// 步骤1：找到优化后种群中的最优个体（性能最高的参数集）
// ---------------------------------------------------------------------
int best_idx = 0;
double best_perf = data[0 * (n_params + 1) + n_params]; // 第一行的性能列
for (int i = 1; i < n_cases; ++i) {
double curr_perf = data[i * (n_params + 1) + n_params];
if (curr_perf > best_perf) {
best_perf = curr_perf;
best_idx = i;
}
}
const double* best_params =
data + best_idx * (n_params + 1); // 最优参数指针
// ---------------------------------------------------------------------
// 步骤2：筛选最优个体附近的样本（仅用近邻样本拟合局部二次曲面）
// ---------------------------------------------------------------------
int nc_kept =
std::min(min_cases, n_cases); // 保留的样本数（不超过总样本数）
double* dist = work_buf; // 工作缓冲区：存储各样本到最优个体的欧氏距离
// 计算每个样本到最优个体的欧氏距离（仅基于参数，不包含性能）
for (int i = 0; i < n_cases; ++i) {
const double* curr_params = data + i * (n_params + 1);
double d = 0.0;
for (int j = 0; j < n_params; ++j) {
double diff = curr_params[j] - best_params[j];
d += diff * diff;
}
dist[i] = d;
}
// 按距离排序，保留最近的nc_kept个样本（局部拟合）
std::vector<int> idx(n_cases);
for (int i = 0; i < n_cases; ++i) idx[i] = i;
std::sort(idx.begin(), idx.end(),
[&](int a, int b) { return dist[a] < dist[b]; });
// ---------------------------------------------------------------------
// 步骤3：构建二次曲面拟合的设计矩阵A与目标向量b
// ---------------------------------------------------------------------
int n_coefs = n_params + (n_params * (n_params + 1)) / 2 + 1; // 总系数数量
double* A = work_buf; // 设计矩阵（nc_kept行 × n_coefs列）
double* b = work_buf + nc_kept * n_coefs; // 目标向量（长度nc_kept）
memset(A, 0.0, nc_kept * n_coefs * sizeof(double));
memset(b, 0.0, nc_kept * sizeof(double));
for (int i = 0; i < nc_kept; ++i) {
int data_idx = idx[i];
const double* curr_params = data + data_idx * (n_params + 1);
double* A_row = A + i * n_coefs;
int coef_idx = 0;
// 1. 一次项：参数中心化（减去最优值，提升数值稳定性）
for (int j = 0; j < n_params; ++j) {
A_row[coef_idx++] = curr_params[j] - best_params[j];
}
// 2. 二次项：参数中心化后的乘积（j<=k，避免重复计算）
for (int j = 0; j < n_params; ++j) {
for (int k = j; k < n_params; ++k) {
double term_j = curr_params[j] - best_params[j];
double term_k = curr_params[k] - best_params[k];
A_row[coef_idx++] = term_j * term_k;
}
}
// 3. 常数项
A_row[coef_idx++] = 1.0;
// 目标向量b：转换为最小化问题（best_perf - 当前样本性能）
b[i] = best_perf - data[data_idx * (n_params + 1) + n_params];
}
// ---------------------------------------------------------------------
// 步骤4：SVD求解最小二乘，得到二次曲面系数
// ---------------------------------------------------------------------
double* U = new (std::nothrow) double[nc_kept * nc_kept]();
double* S = new (std::nothrow) double[std::min(nc_kept, n_coefs)]();
double* V = new (std::nothrow) double[n_coefs * n_coefs]();
double* coefs = new (std::nothrow) double[n_coefs]();
if (!U || !S || !V || !coefs) {
delete[] U;
delete[] S;
delete[] V;
delete[] coefs;
fclose(fp);
return POST_OPT_MEM_ALLOC;
}
// 执行SVD与回代求解
svd(A, nc_kept, n_coefs, U, S, V);
svd_backsub(U, S, V, nc_kept, n_coefs, b, coefs);
// ---------------------------------------------------------------------
// 步骤5：计算Hessian矩阵及其广义逆
// ---------------------------------------------------------------------
compute_hessian(coefs);
hessian_eigen();
// ---------------------------------------------------------------------
// 步骤6：计算参数相对变异（反映参数稳健性：值越大，参数越稳健）
// ---------------------------------------------------------------------
double max_var = 0.0;
double* rel_var = work_buf; // 工作缓冲区：存储相对变异
for (int i = 0; i < n_params; ++i) {
// 变异值 = sqrt(Hessian对角项)（Hessian逆的对角项为方差）
rel_var[i] = (hessian[i * n_params + i] > 1e-10)
? sqrt(hessian[i * n_params + i])
: 0.0;
if (rel_var[i] > max_var) max_var = rel_var[i];
}
// 归一化：最大变异值为1，其他为相对比例
for (int i = 0; i < n_params; ++i) {
rel_var[i] = (max_var > 1e-10) ? rel_var[i] / max_var : 0.0;
}
// ---------------------------------------------------------------------
// 步骤7：计算参数相关性矩阵
// ---------------------------------------------------------------------
double* corr = work_buf + n_params; // 工作缓冲区：存储相关性矩阵
memset(corr, 0.0, n_params * n_params * sizeof(double));
for (int i = 0; i < n_params; ++i) {
double std_i = sqrt(hessian[i * n_params + i]);
if (std_i < 1e-10) continue; // 跳过已置零的参数
for (int j = 0; j < n_params; ++j) {
double std_j = sqrt(hessian[j * n_params + j]);
if (std_j < 1e-10) continue;
// 相关性 = 协方差 / (标准差i * 标准差j)
corr[i * n_params + j] = hessian[i * n_params + j] / (std_i * std_j);
// 截断数值误差：确保相关性在[-1, 1]范围内
corr[i * n_params + j] =
std::max(-1.0, std::min(1.0, corr[i * n_params + j]));
}
}
// ---------------------------------------------------------------------
// 步骤8：分析敏感性方向（特征向量对应最大/最小特征值）
// ---------------------------------------------------------------------
int max_eig_idx = 0, min_eig_idx = 0;
double max_eig = evals[0], min_eig = evals[0];
// 找到最大/最小特征值及其索引
for (int i = 1; i < n_params; ++i) {
if (evals[i] > max_eig) {
max_eig = evals[i];
max_eig_idx = i;
}
if (evals[i] < min_eig) {
min_eig = evals[i];
min_eig_idx = i;
}
}
// 标准化敏感性方向（最大元素为1，便于解释）
double* max_dir =
work_buf + n_params + n_params * n_params; // 最大敏感方向
double* min_dir = max_dir + n_params; // 最小敏感方向
double max_val = 0.0, min_val = 0.0;
for (int i = 0; i < n_params; ++i) {
max_dir[i] = evect[i * n_params + max_eig_idx];
min_dir[i] = evect[i * n_params + min_eig_idx];
max_val = std::max(max_val, fabs(max_dir[i]));
min_val = std::max(min_val, fabs(min_dir[i]));
}
// 归一化：最大绝对值为1
for (int i = 0; i < n_params; ++i) {
max_dir[i] /= (max_val > 1e-10) ? max_val : 1.0;
min_dir[i] /= (min_val > 1e-10) ? min_val : 1.0;
}
// ---------------------------------------------------------------------
// 步骤9：写入日志文件
// ---------------------------------------------------------------------
fprintf(fp, "=== 参数关系分析结果 ===\n");
fprintf(fp, "最优个体性能: %.6f\n", best_perf);
fprintf(fp, "保留样本数: %d / %d\n\n", nc_kept, n_cases);
// 9.1 参数相对变异
fprintf(fp, "1. 参数相对变异（值越大，参数越稳健）\n");
fprintf(fp, " 参数序号: ");
for (int i = 0; i < n_params; ++i) fprintf(fp, "%6d", i + 1);
fprintf(fp, "\n 相对变异: ");
for (int i = 0; i < n_params; ++i) fprintf(fp, "%6.3f", rel_var[i]);
fprintf(fp, "\n\n");
// 9.2 参数相关性矩阵
fprintf(fp, "2. 参数相关性矩阵（-1=强负相关，1=强正相关）\n");
for (int i = 0; i < n_params; ++i) {
fprintf(fp, " 参数%d: ", i + 1);
for (int j = 0; j < n_params; ++j) {
fprintf(fp, "%6.3f", corr[i * n_params + j]);
}
fprintf(fp, "\n");
}
fprintf(fp, "\n");
// 9.3 敏感性方向
fprintf(fp, "3. 敏感性方向（最大=性能变化最快，最小=性能变化最慢）\n");
fprintf(fp, " 参数序号: ");
for (int i = 0; i < n_params; ++i) fprintf(fp, "%6d", i + 1);
fprintf(fp, "\n 最大敏感方向: ");
for (int i = 0; i < n_params; ++i) fprintf(fp, "%6.3f", max_dir[i]);
fprintf(fp, "\n 最小敏感方向: ");
for (int i = 0; i < n_params; ++i) fprintf(fp, "%6.3f", min_dir[i]);
fprintf(fp, "\n");
// ---------------------------------------------------------------------
// 清理临时内存并关闭文件
// ---------------------------------------------------------------------
delete[] U;
delete[] S;
delete[] V;
delete[] coefs;
fclose(fp);
return POST_OPT_SUCCESS;
}
};
#include "ParamCorrelation.h"
// 添加测试函数
int test_param_correlation() {
printf("开始测试 ParamCorrelation 类...\n");
// 准备测试数据 (3个参数, 10个样本)
const int n_params = 3;
const int n_cases = 10;
double test_data[10][4] = {
{1.0, 2.0, 3.0, 95.0}, // 参数1, 参数2, 参数3, 性能
{1.1, 2.2, 3.1, 92.0}, {0.9, 1.8, 2.9, 94.5},
{1.2, 2.1, 3.2, 90.0}, {0.8, 2.3, 2.8, 93.0},
{1.05, 1.95, 3.05, 96.0}, {0.95, 2.05, 2.95, 94.0},
{1.15, 2.25, 3.15, 91.5}, {0.85, 1.75, 2.85, 92.5},
{1.25, 2.35, 3.25, 89.0}};
// 将二维数组转换为一维数组
double flat_data[40];
for (int i = 0; i < n_cases; i++) {
for (int j = 0; j < 4; j++) {
flat_data[i * 4 + j] = test_data[i][j];
}
}
// 创建 ParamCorrelation 对象
ParamCorrelation* analyzer = nullptr;
PostOptErr err =
ParamCorrelation::create(n_params, n_cases, flat_data, &analyzer);
if (err != POST_OPT_SUCCESS) {
printf("创建 ParamCorrelation 对象失败，错误码: %d\n", err);
return err;
}
// 进行分析
err = analyzer->analyze("param_analysis_result.txt");
if (err != POST_OPT_SUCCESS) {
printf("分析过程失败，错误码: %d\n", err);
delete analyzer;
return err;
}
printf("测试完成，结果已保存到 param_analysis_result.txt\n");
delete analyzer;
return POST_OPT_SUCCESS;
}
int main() {
int result = test_param_correlation();
if (result == POST_OPT_SUCCESS) {
printf("测试成功！\n");
} else {
printf("测试失败，错误码：%d\n", result);
}
return result;
}