From 046453ff0a4ed7f1d528d9de59c188243498470d Mon Sep 17 00:00:00 2001 From: lvjunjie Date: Mon, 14 Sep 2026 16:07:44 +0800 Subject: [PATCH] =?UTF-8?q?feat(nmNum):=20=E6=96=B0=E5=A2=9E=E6=97=B6?= =?UTF-8?q?=E9=97=B4=E7=AA=97=E5=8F=A3=20Fisher=20=E5=8F=82=E6=95=B0?= =?UTF-8?q?=E7=AD=9B=E9=80=89?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 优先处理误差能量最大的时间窗口，结合 Fisher 矩阵与残差推荐参数 - 比较最多 3 个推荐参数的组合，按全局预测改善量选择调整方案 - 沿用 LM 步长约束和全局误差验收，无效或改善不足时切换窗口 - 增加全局选参回退，避免部分窗口无效导致提前终止 --- .../nmCalculation/nmCalculationAutoFitLM.h | 2 +- .../nmCalculation/nmCalculationAutoFitLM.cpp | 646 ++++++++---------- 2 files changed, 284 insertions(+), 364 deletions(-) diff --git a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h index deff5329..bea69df0 100644 --- a/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h +++ b/Include/nmNum/nmCalculation/nmCalculationAutoFitLM.h @@ -26,7 +26,7 @@ struct AutoFitTimeWindowLM { }; // 双对数曲线误差分解。total 是 LM 候选接受和排序的唯一依据， -// 其余诊断量用于有限差分灵敏度分析和信赖域选参。 +// 时间窗口和残差用于 Fisher 选参，其余诊断量用于解释曲线失配。 struct AutoFitObjectiveBreakdownLM { bool valid; double total; diff --git a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp index c892d6de..edea6655 100644 --- a/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp +++ b/Src/nmNum/nmCalculation/nmCalculationAutoFitLM.cpp @@ -335,9 +335,8 @@ static int trustRegionDominantComponent( return component; } -// 求解选中参数对应的阻尼正规方程。上下和左右诊断量保留方向；形状没有 -// 天然正负，因此使用 shapeLoss 对参数的局部导数。参数最多七维，使用带 -// 部分主元的高斯消元即可处理该小矩阵，并在主元退化时明确返回失败。 +// 求解选中参数对应的阻尼正规方程。参数最多七维，使用带部分主元的 +// 高斯消元处理该小矩阵，并在主元退化时明确返回失败。 static bool solveTrustRegionLinearSystem( QVector > matrix, QVector rightHandSide, @@ -405,35 +404,218 @@ static bool solveTrustRegionLinearSystem( return true; } -// 计算两个 Jacobian 列向量的绝对余弦相似度。接近 1 表示两个参数在当前 -// 工作点对曲线的影响几乎相同，联合调整容易产生不可辨识方向。 -static double trustRegionJacobianColumnCorrelation( +// Fisher 仅作为当前归一化坐标下的局部信息矩阵，不用于统计置信区间。 +struct TrustRegionFisher +{ + QVector > matrix; + QVector gradient; + + explicit TrustRegionFisher(int dimensions = 0) + : matrix(dimensions, QVector(dimensions, 0.0)) + , gradient(dimensions, 0.0) + {} +}; + +static QVector buildTrustRegionFisher( const QVector >& jacobian, - int leftColumn, - int rightColumn) + const QVector& residual, + const QVector& 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 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(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& windows, + const QVector& 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 selectTrustRegionFisherColumns( + const TrustRegionFisher& local, + const TrustRegionFisher& global, + const QVector& columnValid, + const QVector& 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 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 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& selected, + const QVector& coordinates, + double damping, double trustRadius, double minimumStep, + QVector* step, double* predictedReduction) +{ + // 复用原有全局 LM 方程、阻尼缩放、退化回退、信赖半径和边界投影。 + const int dimensions = coordinates.size(); + const int count = selected.size(); + QVector > normal(count, QVector(count, 0.0)); + QVector 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 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 >* jacobian, @@ -471,33 +653,6 @@ static void updateTrustRegionJacobian( } } -// 对上下偏差、左右偏差和形状损失的梯度执行同样的割线秩一修正，使诊断 -// 选参模型与完整残差 Jacobian 保持在同一个已接受工作点。 -static void updateTrustRegionScalarGradient( - QVector* gradient, - double oldValue, - double newValue, - const QVector& 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 attemptedWindows(kAutoFitTimeWindowCount, false); + bool globalFallbackAttempted = false; StopReasonLM stopReason = LM_MAX_ITERATIONS; // jacobian 的行对应固定 160 维残差，列对应用户勾选的参数。 - // 三个 gradient 单独描述诊断分量对参数的局部变化，只用于本轮选参。 + // Fisher 直接复用残差 Jacobian，上下/左右/形状诊断仅保留用于结果说明。 QVector > jacobian; - QVector verticalGradient(dimensions, 0.0); - QVector horizontalGradient(dimensions, 0.0); - QVector shapeGradient(dimensions, 0.0); QVector 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 >( residualCount, QVector(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* 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 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 information = buildTrustRegionFisher( + jacobian, current.breakdown.residualVector, jacobianColumnValid); + const TrustRegionFisher& global = information[kAutoFitTimeWindowCount]; + const int dominantComponent = trustRegionDominantComponent( + current.breakdown, diagnosisThreshold); // 仅用于现有诊断日志。 QVector selectedColumns; - QVector 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 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 proposed = selectTrustRegionFisherColumns( + local, global, jacobianColumnValid, + current.coordinates, minimumCoordinateStep); + // 最多三个推荐参数，比较其全部非空子集（最多七组），避免首参数必选。 + // 这里只做小矩阵运算，真正的候选评价每轮仍至多一次。 + for(int mask = 1; mask < (1 << proposed.size()); ++mask) { + QVector 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 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 > normalMatrix( - selectedCount, QVector(selectedCount, 0.0)); - QVector 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 selectedStep; - bool solved = solveTrustRegionLinearSystem( - normalMatrix, rightHandSide, &selectedStep); - QVector 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 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 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;