YXS
/
nmWTAI-Platform
8
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

feat(nmNum): 新增时间窗口 Fisher 参数筛选

Browse Source 
- 优先处理误差能量最大的时间窗口，结合 Fisher 矩阵与残差推荐参数
- 比较最多 3 个推荐参数的组合，按全局预测改善量选择调整方案
- 沿用 LM 步长约束和全局误差验收，无效或改善不足时切换窗口
- 增加全局选参回退，避免部分窗口无效导致提前终止
feature/AutoFit-Optimize-20260914
lvjunjie 1 month ago
parent 8f4e35a3ed
commit 046453ff0a
2 changed files with 284 additions and 364 deletions
Show Stats Download Patch File Download Diff File

2
Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h
Unescape Escape View File

@ -26,7 +26,7 @@ struct AutoFitTimeWindowLM {
};
// 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据，
// 其余诊断量用于有限差分灵敏度分析和信赖域选参。
// 时间窗口和残差用于 Fisher 选参，其余诊断量用于解释曲线失配。
struct AutoFitObjectiveBreakdownLM {
bool valid;
double total;

646
Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp
Unescape Escape View File

@ -335,9 +335,8 @@ static int trustRegionDominantComponent(
return component;
}
// 求解选中参数对应的阻尼正规方程。上下和左右诊断量保留方向；形状没有
// 天然正负，因此使用 shapeLoss 对参数的局部导数。参数最多七维，使用带
// 部分主元的高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。
// 求解选中参数对应的阻尼正规方程。参数最多七维，使用带部分主元的
// 高斯消元处理该小矩阵，并在主元退化时明确返回失败。
static bool solveTrustRegionLinearSystem(
QVector<QVector<double> > matrix,
QVector<double> rightHandSide,
@ -405,35 +404,218 @@ static bool solveTrustRegionLinearSystem(
return true;
}
// 计算两个 Jacobian 列向量的绝对余弦相似度。接近 1 表示两个参数在当前
// 工作点对曲线的影响几乎相同，联合调整容易产生不可辨识方向。
static double trustRegionJacobianColumnCorrelation(
// Fisher 仅作为当前归一化坐标下的局部信息矩阵，不用于统计置信区间。
struct TrustRegionFisher
{
QVector<QVector<double> > matrix;
QVector<double> gradient;
explicit TrustRegionFisher(int dimensions = 0)
: matrix(dimensions, QVector<double>(dimensions, 0.0))
, gradient(dimensions, 0.0)
{}
};
static QVector<TrustRegionFisher> buildTrustRegionFisher(
const QVector<QVector<double> >& jacobian,
int leftColumn,
int rightColumn)
const QVector<double>& residual,
const QVector<bool>& columnValid)
{
double product = 0.0;
double leftNorm = 0.0;
double rightNorm = 0.0;
for(int row = 0; row < jacobian.size(); ++row) {
if(leftColumn >= jacobian[row].size() ||
rightColumn >= jacobian[row].size()) {
return 1.0;
// 残差和 J 已包含压力/导数及点数归一化，只再乘一次窗口权重。
// 前半段是压力，后半段是导数；同一时间点的两行使用相同权重。
const int dimensions = columnValid.size();
const int pointCount = residual.size() / 2;
QVector<TrustRegionFisher> information(
kAutoFitTimeWindowCount + 1, TrustRegionFisher(dimensions));
for(int k = 0; k < kAutoFitTimeWindowCount; ++k) {
TrustRegionFisher& local = information[k];
for(int row = 0; row < residual.size(); ++row) {
const double weight = autoFitTimeWindowWeight(
static_cast<double>(row % pointCount) / (pointCount - 1), k);
for(int p = 0; p < dimensions; ++p) {
if(!columnValid[p]) {
continue;
}
local.gradient[p] += weight * jacobian[row][p] * residual[row];
for(int q = 0; q <= p; ++q) {
if(columnValid[q]) {
local.matrix[p][q] +=
weight * jacobian[row][p] * jacobian[row][q];
}
}
}
}
// 互补窗口之和恢复全局 J^T*J 和 J^T*r，不能对各窗口再单独平均。
TrustRegionFisher& global = information[kAutoFitTimeWindowCount];
for(int p = 0; p < dimensions; ++p) {
global.gradient[p] += local.gradient[p];
for(int q = 0; q <= p; ++q) {
local.matrix[q][p] = local.matrix[p][q];
global.matrix[p][q] += local.matrix[p][q];
global.matrix[q][p] = global.matrix[p][q];
}
}
double left = jacobian[row][leftColumn];
double right = jacobian[row][rightColumn];
product += left * right;
leftNorm += left * left;
rightNorm += right * right;
}
return information;
}
if(leftNorm <= 1.0e-20 || rightNorm <= 1.0e-20) {
return 0.0;
static int nextTrustRegionWindow(
const QVector<AutoFitTimeWindowLM>& windows,
const QVector<bool>& attempted)
{
// 每次都按当前误差能量重新排序；本轮无效的窗口暂时跳过，-1 表示全局回退。
int selected = -1;
double largestEnergy = 0.0;
for(int k = 0; k < windows.size(); ++k) {
if(!attempted[k] && isFiniteNumber(windows[k].energy) &&
windows[k].energy > largestEnergy) {
largestEnergy = windows[k].energy;
selected = k;
}
}
return qAbs(product) / qSqrt(leftNorm * rightNorm);
return selected;
}
// 每次接受一个真实候选后，使用满足最新割线条件的秩一修正更新完整残差
static QVector<int> selectTrustRegionFisherColumns(
const TrustRegionFisher& local,
const TrustRegionFisher& global,
const QVector<bool>& columnValid,
const QVector<double>& coordinates,
double minimumStep)
{
// g^2/F 衡量修正当前残差的潜力，避免仅按影响强度选中无助于降误差的参数。
const int dimensions = columnValid.size();
double maximumDiagonal = 0.0;
for(int p = 0; p < dimensions; ++p) {
if(columnValid[p] && isFiniteNumber(local.matrix[p][p])) {
maximumDiagonal = qMax(maximumDiagonal, local.matrix[p][p]);
}
}
const double diagonalFloor = qMax(1.0e-30, maximumDiagonal * 1.0e-12);
QVector<double> scores(dimensions, 0.0);
double maximumScore = 0.0;
for(int p = 0; p < dimensions; ++p) {
if(!columnValid[p] || local.matrix[p][p] <= diagonalFloor ||
!isFiniteNumber(local.matrix[p][p]) ||
!isFiniteNumber(local.gradient[p]) ||
!isFiniteNumber(global.gradient[p])) {
continue;
}
// 全局下降方向指向边界外时暂不选该参数，避免它挤占可移动参数的位置。
if((coordinates[p] <= minimumStep && global.gradient[p] > 0.0) ||
(coordinates[p] >= 1.0 - minimumStep && global.gradient[p] < 0.0)) {
continue;
}
double scaledGradient = local.gradient[p] /
qSqrt(local.matrix[p][p] + diagonalFloor);
double score = scaledGradient * scaledGradient;
if(isFiniteNumber(score)) {
scores[p] = score;
maximumScore = qMax(maximumScore, score);
}
}
QVector<int> selected;
for(int selection = 0; selection < qMin(3, dimensions); ++selection) {
int best = -1;
double bestScore = maximumScore * 1.0e-12;
for(int p = 0; p < dimensions; ++p) {
if(selected.contains(p) || scores[p] <= bestScore) {
continue;
}
// 沿用现有全局 0.995 共线限制，不因单个窗口内响应相似而冻结参数。
bool correlated = false;
for(int i = 0; i < selected.size(); ++i) {
int q = selected[i];
double denominator = qSqrt(global.matrix[p][p]) *
qSqrt(global.matrix[q][q]);
if(denominator > 0.0 &&
qAbs(global.matrix[p][q]) / denominator > 0.995) {
correlated = true;
break;
}
}
if(!correlated) {
best = p;
bestScore = scores[p];
}
}
if(best < 0) {
break;
}
selected.append(best);
}
return selected;
}
static bool buildTrustRegionFisherStep(
const TrustRegionFisher& global,
const QVector<int>& selected,
const QVector<double>& coordinates,
double damping, double trustRadius, double minimumStep,
QVector<double>* step, double* predictedReduction)
{
// 复用原有全局 LM 方程、阻尼缩放、退化回退、信赖半径和边界投影。
const int dimensions = coordinates.size();
const int count = selected.size();
QVector<QVector<double> > normal(count, QVector<double>(count, 0.0));
QVector<double> rhs(count, 0.0);
for(int i = 0; i < count; ++i) {
const int p = selected[i];
rhs[i] = -global.gradient[p];
for(int j = 0; j < count; ++j) {
normal[i][j] = global.matrix[p][selected[j]];
}
normal[i][i] += damping * qMax(1.0e-10, normal[i][i]);
}
QVector<double> solution;
bool solved = solveTrustRegionLinearSystem(normal, rhs, &solution);
step->fill(0.0, dimensions);
if(solved) {
for(int i = 0; i < count; ++i) {
(*step)[selected[i]] = solution[i];
}
}
double norm = qSqrt(trustRegionSquaredNorm(*step));
if(!solved || !isFiniteNumber(norm) || norm < minimumStep) {
step->fill(0.0, dimensions);
for(int i = 0; i < count; ++i) {
const int p = selected[i];
double direction = -global.gradient[p];
if((coordinates[p] <= minimumStep && direction < 0.0) ||
(coordinates[p] >= 1.0 - minimumStep && direction > 0.0)) {
direction = 0.0;
}
(*step)[p] = direction;
}
norm = qSqrt(trustRegionSquaredNorm(*step));
if(norm > minimumStep) {
for(int p = 0; p < dimensions; ++p) {
(*step)[p] *= trustRadius / norm;
}
norm = trustRadius;
}
}
if(norm > trustRadius) {
for(int p = 0; p < dimensions; ++p) {
(*step)[p] *= trustRadius / norm;
}
}
for(int p = 0; p < dimensions; ++p) {
(*step)[p] = qBound(0.0, coordinates[p] + (*step)[p], 1.0) - coordinates[p];
}
// 只比较实际可执行步长的全局预测下降，阻尼项不属于真实拟合目标。
*predictedReduction = -trustRegionDotProduct(global.gradient, *step);
for(int p = 0; p < dimensions; ++p) {
for(int q = 0; q < dimensions; ++q) {
*predictedReduction -= 0.5 * (*step)[p] * global.matrix[p][q] * (*step)[q];
}
}
return isFiniteNumber(*predictedReduction) && *predictedReduction > 1.0e-14 &&
qSqrt(trustRegionSquaredNorm(*step)) >= minimumStep;
}
// 每次得到有效真实候选后，使用满足最新割线条件的秩一修正更新完整残差
// Jacobian。这样模型吸收了刚得到的真实变化，又不必立即逐参数重新试算。
static void updateTrustRegionJacobian(
QVector<QVector<double> >* jacobian,
@ -471,33 +653,6 @@ static void updateTrustRegionJacobian(
}
}
// 对上下偏差、左右偏差和形状损失的梯度执行同样的割线秩一修正，使诊断
// 选参模型与完整残差 Jacobian 保持在同一个已接受工作点。
static void updateTrustRegionScalarGradient(
QVector<double>* gradient,
double oldValue,
double newValue,
const QVector<double>& coordinateStep)
{
if(!gradient || gradient->size() != coordinateStep.size() ||
!isFiniteNumber(oldValue) || !isFiniteNumber(newValue)) {
return;
}
double denominator = trustRegionSquaredNorm(coordinateStep);
if(denominator <= 1.0e-12) {
return;
}
double predictedChange = trustRegionDotProduct(
*gradient, coordinateStep);
double correction =
(newValue - oldValue - predictedChange) / denominator;
for(int i = 0; i < gradient->size(); ++i) {
(*gradient)[i] += correction * coordinateStep[i];
}
}
static QStringList traceParameterNames()
{
QStringList names;
@ -1588,7 +1743,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
const double minimumCoordinateStep = 1.0e-5;
const double minimumTrustRadius = 2.0e-3;
const double maximumTrustRadius = 0.30;
const double columnCorrelationLimit = 0.995;
const double diagnosisThreshold = 1.0e-5;
// 误差下降至少达到绝对 1e-5 且相对当前有效基准 0.2% 才算有效改善。
// 更小的下降仍保留为最佳解，但不能反复清除停滞状态、延长拟合时间。
@ -1608,14 +1762,13 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
bool rebuildRequested = true;
bool modelRebuiltAtMinimumRadius = false;
bool stagnationConfirmationRequested = false;
QVector<bool> attemptedWindows(kAutoFitTimeWindowCount, false);
bool globalFallbackAttempted = false;
StopReasonLM stopReason = LM_MAX_ITERATIONS;
// jacobian 的行对应固定 160 维残差，列对应用户勾选的参数。
// 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。
// Fisher 直接复用残差 Jacobian，上下/左右/形状诊断仅保留用于结果说明。
QVector<QVector<double> > jacobian;
QVector<double> verticalGradient(dimensions, 0.0);
QVector<double> horizontalGradient(dimensions, 0.0);
QVector<double> shapeGradient(dimensions, 0.0);
QVector<bool> jacobianColumnValid(dimensions, false);
// 参数向量的顺序始终与 m_enabledParamIndices 一致，不能按完整参数索引
@ -1734,16 +1887,24 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
return false;
}
// 有效改善后重新开放全部窗口，并按新工作点的能量排序。
attemptedWindows.fill(false);
globalFallbackAttempted = false;
effectiveImprovementBaseline = fitness;
consecutiveIneffectiveSteps = 0;
stagnationConfirmationRequested = false;
return true;
};
// 连续三次没有有效改善时只请求一次灵敏度重建。重建完成后由主循环
// 直接检查累计改善，仍达不到门槛就判定局部收敛，不再继续微小试探。
// 无效尝试仍累计，但所有窗口及一次全局回退尝试完毕前，不触发停滞收敛。
// 完整一轮仍无改善时沿用原有重建确认，避免某个难处理窗口提前终止拟合。
auto recordIneffectiveStep = [&]() -> bool {
++consecutiveIneffectiveSteps;
if(!globalFallbackAttempted) {
return false;
}
attemptedWindows.fill(false);
globalFallbackAttempted = false;
if(consecutiveIneffectiveSteps < maximumIneffectiveSteps) {
return false;
}
@ -1752,13 +1913,13 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
return true;
}
consecutiveIneffectiveSteps = 0;
stagnationConfirmationRequested = true;
rebuildRequested = true;
emit logMessageGenerated(
tr("No effective improvement for %1 consecutive steps; "
"rebuilding sensitivity model for confirmation")
.arg(maximumIneffectiveSteps));
.arg(consecutiveIneffectiveSteps));
consecutiveIneffectiveSteps = 0;
return false;
};
@ -1784,9 +1945,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
jacobian = QVector<QVector<double> >(
residualCount, QVector<double>(dimensions, 0.0));
verticalGradient.fill(0.0, dimensions);
horizontalGradient.fill(0.0, dimensions);
shapeGradient.fill(0.0, dimensions);
jacobianColumnValid.fill(false, dimensions);
TrustRegionEvaluation bestProbe;
@ -1872,27 +2030,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
base.breakdown.residualVector[row]) / delta;
}
// 有符号诊断量只有在基点和试算点都可靠时才能计算方向梯度；
// shapeLoss 无方向可靠性标志，始终记录其局部变化率。
if(base.breakdown.verticalReliable &&
probe.breakdown.verticalReliable &&
!base.breakdown.registrationAmbiguous &&
!probe.breakdown.registrationAmbiguous) {
verticalGradient[column] =
(probe.breakdown.verticalCommonBias -
base.breakdown.verticalCommonBias) / delta;
}
if(base.breakdown.horizontalReliable &&
probe.breakdown.horizontalReliable &&
!base.breakdown.registrationAmbiguous &&
!probe.breakdown.registrationAmbiguous) {
horizontalGradient[column] =
(probe.breakdown.horizontalPhysicalShift -
base.breakdown.horizontalPhysicalShift) / delta;
}
shapeGradient[column] =
(probe.breakdown.shapeLoss -
base.breakdown.shapeLoss) / delta;
jacobianColumnValid[column] = true;
columnBuilt = true;
@ -1928,28 +2065,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
base.breakdown.residualVector,
bestProbe.breakdown.residualVector,
acceptedStep);
if(base.breakdown.verticalReliable &&
bestProbe.breakdown.verticalReliable) {
updateTrustRegionScalarGradient(
&verticalGradient,
base.breakdown.verticalCommonBias,
bestProbe.breakdown.verticalCommonBias,
acceptedStep);
}
if(base.breakdown.horizontalReliable &&
bestProbe.breakdown.horizontalReliable) {
updateTrustRegionScalarGradient(
&horizontalGradient,
base.breakdown.horizontalPhysicalShift,
bestProbe.breakdown.horizontalPhysicalShift,
acceptedStep);
}
updateTrustRegionScalarGradient(
&shapeGradient,
base.breakdown.shapeLoss,
bestProbe.breakdown.shapeLoss,
acceptedStep);
current = bestProbe;
publishAcceptedPoint(current);
restoreEvaluationState(current);
@ -2025,126 +2140,64 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
}
}
// 先确定当前最突出的可靠诊断误差，用其梯度回答“哪些参数最能改善
// 当前问题”；实际 LM 方向仍由完整残差梯度和 Jacobian 共同计算。
int dominantComponent = trustRegionDominantComponent(
current.breakdown, diagnosisThreshold);
const QVector<double>* componentGradient = nullptr;
if(dominantComponent == TRUST_REGION_VERTICAL_COMPONENT) {
componentGradient = &verticalGradient;
} else if(dominantComponent == TRUST_REGION_HORIZONTAL_COMPONENT) {
componentGradient = &horizontalGradient;
} else if(dominantComponent == TRUST_REGION_SHAPE_COMPONENT) {
componentGradient = &shapeGradient;
}
// 主目标采用 0.5*||r||^2，其对参数的梯度为 J^T*r。这里不再叠加
// vertical/horizontal/shape，保证诊断分量不会改变真实接受目标。
QVector<double> totalGradient(dimensions, 0.0);
for(int column = 0; column < dimensions; ++column) {
if(!jacobianColumnValid[column]) {
continue;
}
for(int row = 0; row < jacobian.size(); ++row) {
totalGradient[column] +=
jacobian[row][column] *
current.breakdown.residualVector[row];
}
}
// 每轮最多联合调整三个灵敏参数。按当前诊断梯度绝对值由大到小选取，
// 并剔除 Jacobian 响应过度共线的列，降低弱可辨识参数互相补偿的风险。
// 每轮从最新 J 和当前残差重算窗口 Fisher，包含有限差分与割线更新的变化。
const QVector<TrustRegionFisher> information = buildTrustRegionFisher(
jacobian, current.breakdown.residualVector, jacobianColumnValid);
const TrustRegionFisher& global = information[kAutoFitTimeWindowCount];
const int dominantComponent = trustRegionDominantComponent(
current.breakdown, diagnosisThreshold); // 仅用于现有诊断日志。
QVector<int> selectedColumns;
QVector<bool> alreadyConsidered(dimensions, false);
for(int selection = 0; selection < qMin(3, dimensions); ++selection) {
int bestColumn = -1;
double bestScore = 0.0;
for(int column = 0; column < dimensions; ++column) {
if(alreadyConsidered[column] ||
!jacobianColumnValid[column]) {
continue;
QVector<double> coordinateStep;
double predictedReduction = 0.0;
int selectedWindow = -1;
// 最大能量窗口优先。某窗口无可行全局下降步时立即换下一个，不调用 DLL。
// 全部窗口都处理不动后，再用全局 Fisher 作一次补充选参。
while(selectedColumns.isEmpty()) {
selectedWindow = nextTrustRegionWindow(
current.breakdown.timeWindows, attemptedWindows);
if(selectedWindow >= 0) {
attemptedWindows[selectedWindow] = true;
} else {
globalFallbackAttempted = true;
}
const TrustRegionFisher& local = selectedWindow >= 0
? information[selectedWindow] : global;
const QVector<int> proposed = selectTrustRegionFisherColumns(
local, global, jacobianColumnValid,
current.coordinates, minimumCoordinateStep);
// 最多三个推荐参数，比较其全部非空子集（最多七组），避免首参数必选。
// 这里只做小矩阵运算，真正的候选评价每轮仍至多一次。
for(int mask = 1; mask < (1 << proposed.size()); ++mask) {
QVector<int> columns;
for(int i = 0; i < proposed.size(); ++i) {
if(mask & (1 << i)) {
columns.append(proposed[i]);
}
}
double score = componentGradient
? qAbs((*componentGradient)[column])
: qAbs(totalGradient[column]);
if(!isFiniteNumber(score) || score <= bestScore) {
QVector<double> step;
double reduction = 0.0;
if(!buildTrustRegionFisherStep(
global, columns, current.coordinates,
damping, trustRadius, minimumCoordinateStep,
&step, &reduction)) {
continue;
}
bool excessivelyCorrelated = false;
for(int selectedIndex = 0;
selectedIndex < selectedColumns.size();
++selectedIndex) {
if(trustRegionJacobianColumnCorrelation(
jacobian,
column,
selectedColumns[selectedIndex]) >
columnCorrelationLimit) {
excessivelyCorrelated = true;
break;
}
}
if(!excessivelyCorrelated) {
bestColumn = column;
bestScore = score;
const double tolerance = 1.0e-12 * qMax(predictedReduction, reduction);
if(selectedColumns.isEmpty() || reduction > predictedReduction + tolerance ||
(qAbs(reduction - predictedReduction) <= tolerance &&
columns.size() < selectedColumns.size())) {
selectedColumns = columns;
coordinateStep = step;
predictedReduction = reduction;
}
}
if(bestColumn < 0 || bestScore <= 1.0e-12) {
if(selectedWindow < 0) {
break;
}
selectedColumns.append(bestColumn);
alreadyConsidered[bestColumn] = true;
}
// 诊断梯度接近零时，说明该分量在当前局部无法可靠选参，退回完整残差
// 梯度，但接受标准仍然只有真实 total，诊断值不会重复计入目标函数。
if(selectedColumns.isEmpty() && componentGradient) {
dominantComponent = TRUST_REGION_TOTAL_COMPONENT;
componentGradient = nullptr;
alreadyConsidered.fill(false, dimensions);
for(int selection = 0;
selection < qMin(3, dimensions);
++selection) {
int bestColumn = -1;
double bestScore = 0.0;
for(int column = 0; column < dimensions; ++column) {
if(alreadyConsidered[column] ||
!jacobianColumnValid[column]) {
continue;
}
double score = qAbs(totalGradient[column]);
if(score <= bestScore) {
continue;
}
bool excessivelyCorrelated = false;
for(int selectedIndex = 0;
selectedIndex < selectedColumns.size();
++selectedIndex) {
if(trustRegionJacobianColumnCorrelation(
jacobian,
column,
selectedColumns[selectedIndex]) >
columnCorrelationLimit) {
excessivelyCorrelated = true;
break;
}
}
if(!excessivelyCorrelated) {
bestColumn = column;
bestScore = score;
}
}
if(bestColumn < 0 || bestScore <= 1.0e-12) {
break;
}
selectedColumns.append(bestColumn);
alreadyConsidered[bestColumn] = true;
}
}
// 当前局部没有可用方向时先缩小半径并重建灵敏度；只有已经在最小
// 半径完整重建后仍无方向，才把它判定为局部最优。
// 只有窗口候选与全局回退均无方向，才收缩半径并进入原有重建/收敛处理。
if(selectedColumns.isEmpty()) {
if(trustRadius <= minimumTrustRadius * 1.01 &&
modelRebuiltAtMinimumRadius) {
@ -2161,126 +2214,19 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
continue;
}
// 在选中参数子空间构造 LM 正规方程：
// (J^T*J + damping*diag(J^T*J))*step = -J^T*r。
// 对角缩放使不同参数列的灵敏度量级差异不会直接改变阻尼强弱。
const int selectedCount = selectedColumns.size();
QVector<QVector<double> > normalMatrix(
selectedCount, QVector<double>(selectedCount, 0.0));
QVector<double> rightHandSide(selectedCount, 0.0);
for(int left = 0; left < selectedCount; ++left) {
int leftColumn = selectedColumns[left];
rightHandSide[left] = -totalGradient[leftColumn];
for(int right = 0; right < selectedCount; ++right) {
int rightColumn = selectedColumns[right];
for(int row = 0; row < jacobian.size(); ++row) {
normalMatrix[left][right] +=
jacobian[row][leftColumn] *
jacobian[row][rightColumn];
}
}
double diagonalScale = qMax(
1.0e-10, normalMatrix[left][left]);
normalMatrix[left][left] += damping * diagonalScale;
}
QVector<double> selectedStep;
bool solved = solveTrustRegionLinearSystem(
normalMatrix, rightHandSide, &selectedStep);
QVector<double> coordinateStep(dimensions, 0.0);
if(solved) {
for(int i = 0; i < selectedCount; ++i) {
coordinateStep[selectedColumns[i]] = selectedStep[i];
}
}
double stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
if(!solved || !isFiniteNumber(stepNorm) ||
stepNorm < minimumCoordinateStep) {
// 正规方程退化时使用投影最速下降方向，仍只移动本轮已选择的参数。
coordinateStep.fill(0.0, dimensions);
double gradientNormSquared = 0.0;
for(int i = 0; i < selectedCount; ++i) {
int column = selectedColumns[i];
double stepDirection = -totalGradient[column];
if((current.coordinates[column] <= minimumCoordinateStep &&
stepDirection < 0.0) ||
(current.coordinates[column] >=
1.0 - minimumCoordinateStep &&
stepDirection > 0.0)) {
stepDirection = 0.0;
}
coordinateStep[column] = stepDirection;
gradientNormSquared += stepDirection * stepDirection;
}
double gradientNorm = qSqrt(gradientNormSquared);
if(gradientNorm > minimumCoordinateStep) {
double scale = trustRadius / gradientNorm;
for(int i = 0; i < selectedCount; ++i) {
int column = selectedColumns[i];
coordinateStep[column] *= scale;
}
}
stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
}
// LM 解只给出局部模型建议方向；若超出当前信赖半径，保持方向不变并
// 等比例截短，避免一次试算离开 Jacobian 有效的局部区域。
if(stepNorm > trustRadius && stepNorm > 0.0) {
double scale = trustRadius / stepNorm;
for(int i = 0; i < coordinateStep.size(); ++i) {
coordinateStep[i] *= scale;
}
}
// 将 LM 步长投影到用户给定的参数范围，实际用于预测下降的也是投影后步长。
const double stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
QVector<double> candidateCoordinates = current.coordinates;
for(int i = 0; i < dimensions; ++i) {
candidateCoordinates[i] = qBound(
0.0,
current.coordinates[i] + coordinateStep[i],
1.0);
coordinateStep[i] = candidateCoordinates[i] -
current.coordinates[i];
}
stepNorm = qSqrt(trustRegionSquaredNorm(coordinateStep));
// 用线性模型 r_new ~= r_current + J*step 预测残差，再用平方能量
// 的下降量与真实候选下降量比较，作为调整阻尼和半径的依据。
QVector<double> predictedResidual =
current.breakdown.residualVector;
for(int row = 0; row < jacobian.size(); ++row) {
for(int column = 0; column < dimensions; ++column) {
predictedResidual[row] +=
jacobian[row][column] * coordinateStep[column];
}
for(int column = 0; column < dimensions; ++column) {
candidateCoordinates[column] += coordinateStep[column];
}
double predictedReduction = 0.5 *
(trustRegionSquaredNorm(current.breakdown.residualVector) -
trustRegionSquaredNorm(predictedResidual));
// 无实际移动或模型预测不下降时没有必要调用昂贵求解器。将它按一次
// 拒绝处理，并在连续发生后重建灵敏度，防止继续沿失效模型试算。
if(stepNorm < minimumCoordinateStep ||
!isFiniteNumber(predictedReduction) ||
predictedReduction <= 1.0e-14) {
trustRadius = qMax(minimumTrustRadius, trustRadius * 0.5);
damping = qMin(1.0e8, damping * 4.0);
++consecutiveRejectedSteps;
if(consecutiveRejectedSteps >= 2) {
if(trustRadius <= minimumTrustRadius * 1.01 &&
modelRebuiltAtMinimumRadius) {
stopReason = LM_LOCAL_OPTIMUM;
break;
}
rebuildRequested = true;
}
if(recordIneffectiveStep()) {
stopReason = LM_LOCAL_OPTIMUM;
break;
}
continue;
// 窗口和参数索引仅写入已有跟踪文件，不向拟合日志窗口增加分段信息。
QStringList selectedParameterIndices;
for(int i = 0; i < selectedColumns.size(); ++i) {
selectedParameterIndices << QString::number(m_enabledParamIndices[selectedColumns[i]]);
}
const QString selectionName = (selectedWindow >= 0
? QString("window_%1").arg(selectedWindow + 1) : QString("global")) +
"_params_" + selectedParameterIndices.join("_");
TrustRegionEvaluation candidate;
candidate.coordinates = candidateCoordinates;
@ -2306,7 +2252,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
candidate.fitness,
false,
candidate.elapsedMs,
"solver_invalid",
"solver_invalid_" + selectionName,
nullptr);
restoreEvaluationState(current);
if(consecutiveRejectedSteps >= 2) {
@ -2332,32 +2278,6 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
oldBreakdown.residualVector,
candidate.breakdown.residualVector,
coordinateStep);
if(oldBreakdown.verticalReliable &&
candidate.breakdown.verticalReliable &&
!oldBreakdown.registrationAmbiguous &&
!candidate.breakdown.registrationAmbiguous) {
updateTrustRegionScalarGradient(
&verticalGradient,
oldBreakdown.verticalCommonBias,
candidate.breakdown.verticalCommonBias,
coordinateStep);
}
if(oldBreakdown.horizontalReliable &&
candidate.breakdown.horizontalReliable &&
!oldBreakdown.registrationAmbiguous &&
!candidate.breakdown.registrationAmbiguous) {
updateTrustRegionScalarGradient(
&horizontalGradient,
oldBreakdown.horizontalPhysicalShift,
candidate.breakdown.horizontalPhysicalShift,
coordinateStep);
}
updateTrustRegionScalarGradient(
&shapeGradient,
oldBreakdown.shapeLoss,
candidate.breakdown.shapeLoss,
coordinateStep);
// reductionRatio 衡量局部线性模型的可信度：接近 1 表示预测准确；
// 值较小表示虽然可能下降，但模型低估了非线性，需要收紧下一步。
double actualReduction = 0.5 *
@ -2424,8 +2344,8 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
true,
candidate.elapsedMs,
accepted
? QString("accepted_%1").arg(componentName)
: QString("rejected_%1").arg(componentName),
? "accepted_" + selectionName
: "rejected_" + selectionName,
&candidate.breakdown);
QString componentDisplayName = componentName;
@ -2456,7 +2376,7 @@ StopReasonLM nmCalculationAutoFitLM::runTrustRegionFitting()
stopReason = LM_TARGET_ACHIEVED;
break;
}
if(trustRadius <= minimumTrustRadius * 1.01 &&
if(selectedWindow < 0 && trustRadius <= minimumTrustRadius * 1.01 &&
consecutiveRejectedSteps >= 2) {
if(modelRebuiltAtMinimumRadius) {
stopReason = LM_LOCAL_OPTIMUM;

Write Preview
Loading…
Cancel
Save