adaptive step + retry for num grad

Author: hexqnt

src/fit.rs

@@ -138,6 +138,8 @@ const MAX_POLYNOMIAL_PARAMS: usize = 10;
 const STEEPEST_DESCENT_GRAD_TOL: f64 = 1e-12;
 const HESSIAN_FD_REL_STEP: f64 = 1e-4;
 const HESSIAN_FD_MIN_STEP: f64 = 1e-6;
+// Пробуем базовый шаг, затем уменьшаем/увеличиваем его, чтобы переживать локальные NaN/Inf.
+const FD_STEP_RETRY_FACTORS: [f64; 5] = [1.0, 0.5, 2.0, 0.25, 4.0];
 pub(crate) const HESSIAN_DIAGONAL_JITTER: f64 = models::HESSIAN_DIAGONAL_JITTER;
 
 fn positive_x(value: f64) -> f64 {
@@ -153,6 +155,24 @@ fn gradient_l2_norm(values: &[f64]) -> f64 {
     values.iter().map(|value| value * value).sum::<f64>().sqrt()
 }
 
+#[inline]
+fn finite_central_difference(value_plus: f64, value_minus: f64, step: f64) -> Option<f64> {
+    if !value_plus.is_finite() || !value_minus.is_finite() {
+        return None;
+    }
+    let derivative = (value_plus - value_minus) / (2.0 * step);
+    if derivative.is_finite() {
+        Some(derivative)
+    } else {
+        None
+    }
+}
+
+#[inline]
+fn finite_array1(values: &Array1<f64>) -> bool {
+    values.iter().all(|value| value.is_finite())
+}
+
 fn vec_to_array1(values: &[f64]) -> Array1<f64> {
     Array1::from_vec(values.to_vec())
 }
@@ -194,20 +214,48 @@ where
     let dimension = param_slice.len();
     let mut hessian = Array2::zeros((dimension, dimension));
     let mut probe = param.clone();
+    let mut column_values = vec![0.0; dimension];
 
     for column in 0..dimension {
-        let step =
+        let base_step =
             ((param_slice[column].abs() + 1.0) * HESSIAN_FD_REL_STEP).max(HESSIAN_FD_MIN_STEP);
-        probe[column] = param[column] + step;
-        let grad_plus = problem.gradient(&probe)?;
-        probe[column] = param[column] - step;
-        let grad_minus = problem.gradient(&probe)?;
-        probe[column] = param[column];
+        let mut computed = false;
+        for factor in FD_STEP_RETRY_FACTORS {
+            let step = base_step * factor;
+            probe[column] = param[column] + step;
+            let grad_plus = problem.gradient(&probe)?;
+            probe[column] = param[column] - step;
+            let grad_minus = problem.gradient(&probe)?;
+            probe[column] = param[column];
+
+            if !finite_array1(&grad_plus) || !finite_array1(&grad_minus) {
+                continue;
+            }
 
-        let denom = 2.0 * step;
-        for row in 0..dimension {
-            let value = (grad_plus[row] - grad_minus[row]) / denom;
-            hessian[[row, column]] = if value.is_finite() { value } else { 0.0 };
+            let denom = 2.0 * step;
+            let mut column_is_finite = true;
+            for row in 0..dimension {
+                let value = (grad_plus[row] - grad_minus[row]) / denom;
+                if !value.is_finite() {
+                    column_is_finite = false;
+                    break;
+                }
+                column_values[row] = value;
+            }
+
+            if column_is_finite {
+                for row in 0..dimension {
+                    hessian[[row, column]] = column_values[row];
+                }
+                computed = true;
+                break;
+            }
+        }
+
+        if !computed {
+            for row in 0..dimension {
+                hessian[[row, column]] = 0.0;
+            }
         }
     }
 
@@ -2023,19 +2071,18 @@ impl Gradient for SplineProblem {
         let mut gradient = Array1::zeros(param_slice.len());
         for (index, gradient_value) in gradient.iter_mut().enumerate() {
             // Численный градиент по центральной схеме конечной разности.
-            let step =
+            let base_step =
                 ((param_slice[index].abs() + 1.0) * SPLINE_FD_REL_STEP).max(SPLINE_FD_MIN_STEP);
-            probe[index] = param_slice[index] + step;
-            let cost_plus = self.evaluate_objective(array1_as_slice(&probe));
-            probe[index] = param_slice[index] - step;
-            let cost_minus = self.evaluate_objective(array1_as_slice(&probe));
-            probe[index] = param_slice[index];
-            let derivative = (cost_plus - cost_minus) / (2.0 * step);
-            *gradient_value = if derivative.is_finite() {
-                derivative
-            } else {
-                LARGE_COST
-            };
+            let derivative = FD_STEP_RETRY_FACTORS.iter().copied().find_map(|factor| {
+                let step = base_step * factor;
+                probe[index] = param_slice[index] + step;
+                let cost_plus = self.evaluate_objective(array1_as_slice(&probe));
+                probe[index] = param_slice[index] - step;
+                let cost_minus = self.evaluate_objective(array1_as_slice(&probe));
+                probe[index] = param_slice[index];
+                finite_central_difference(cost_plus, cost_minus, step)
+            });
+            *gradient_value = derivative.unwrap_or(LARGE_COST);
         }
         Ok(gradient)
     }

src/fit/tests.rs

@@ -1,24 +1,26 @@
 use super::simd;
 use super::{
     CurveProblem, CurveProblemPredictionLoss, DEFAULT_SPLINE_KNOTS, FitError,
-    IncrementalSplineFitRunner, IncrementalSplineFitStep, MetricQuantization,
-    MetricQuantizationDecimalPlaces, OptimizationLossMetric, SplineConfig, SplineDuplicateXPolicy,
-    SplineExtrapolation, SplineFamilyKind, SplineFinalizeContext, SplineKnotStrategy,
-    approximate_spline_knots, build_spline_initial_curve_from_knot_y,
-    build_spline_result_from_knot_y, calculate_iteration_metrics,
-    calculate_iteration_metrics_with_quantization, calculate_metrics, evaluate_linear_spline,
-    expanded_spline_curve_x_bounds, fit_akima_spline, fit_akima_spline_with_config, fit_curve,
-    fit_curve_with_optimizer_config, fit_curve_with_progress,
-    fit_curve_with_progress_and_optimizer_config,
+    HESSIAN_DIAGONAL_JITTER, HESSIAN_FD_MIN_STEP, HESSIAN_FD_REL_STEP, IncrementalSplineFitRunner,
+    IncrementalSplineFitStep, MetricQuantization, MetricQuantizationDecimalPlaces,
+    OptimizationLossMetric, SplineConfig, SplineDuplicateXPolicy, SplineExtrapolation,
+    SplineFamilyKind, SplineFinalizeContext, SplineKnotStrategy, approximate_spline_knots,
+    build_spline_initial_curve_from_knot_y, build_spline_result_from_knot_y,
+    calculate_iteration_metrics, calculate_iteration_metrics_with_quantization, calculate_metrics,
+    evaluate_linear_spline, expanded_spline_curve_x_bounds, fit_akima_spline,
+    fit_akima_spline_with_config, fit_curve, fit_curve_with_optimizer_config,
+    fit_curve_with_progress, fit_curve_with_progress_and_optimizer_config,
     fit_curve_with_progress_and_optimizer_config_and_loss_metric, fit_linear_spline,
     fit_linear_spline_with_config, fit_monotone_cubic_spline, fit_natural_cubic_spline,
-    sorted_points_with_duplicate_policy,
+    numerical_hessian_from_gradient, sorted_points_with_duplicate_policy,
 };
 use crate::domain::{
     AdamConfig, CurveFamily, CurveParams, InputError, LbfgsConfig, NelderMeadConfig,
     NewtonCgConfig, OptimizerConfig, Point, Points, SgdConfig, SteepestDescentConfig,
 };
 use crate::models::{self, ObjectiveGrad, ObjectiveHessian, ObjectiveValue, PredictionLoss};
+use argmin::core::Gradient;
+use ndarray::Array1;
 
 #[derive(Clone, Copy)]
 struct MsePredictionLoss;
@@ -65,6 +67,115 @@ fn quantization(decimal_places: u8) -> MetricQuantization {
     )
 }
 
+struct RetryGradientProblem {
+    center: f64,
+    invalid_step: f64,
+}
+
+struct AlwaysInvalidGradientProblem;
+
+impl Gradient for RetryGradientProblem {
+    type Param = Array1<f64>;
+    type Gradient = Array1<f64>;
+
+    fn gradient(&self, param: &Self::Param) -> Result<Self::Gradient, argmin::core::Error> {
+        let delta = (param[0] - self.center).abs();
+        if (delta - self.invalid_step).abs() <= 1e-14 {
+            return Ok(Array1::from_vec(vec![f64::NAN]));
+        }
+        Ok(Array1::from_vec(vec![2.0 * param[0]]))
+    }
+}
+
+impl Gradient for AlwaysInvalidGradientProblem {
+    type Param = Array1<f64>;
+    type Gradient = Array1<f64>;
+
+    fn gradient(&self, param: &Self::Param) -> Result<Self::Gradient, argmin::core::Error> {
+        Ok(Array1::from_vec(vec![f64::NAN; param.len()]))
+    }
+}
+
+#[test]
+fn central_diff_gradient_retries_step_when_primary_step_is_invalid() {
+    let param = [1.0_f64];
+    let rel_step = 1e-3;
+    let min_step = 1e-3;
+    let base_step = ((param[0].abs() + 1.0) * rel_step).max(min_step);
+    let mut gradient = [0.0];
+    models::central_diff_gradient_from_value(
+        &param,
+        rel_step,
+        min_step,
+        |probe| {
+            let delta = (probe[0] - param[0]).abs();
+            if (delta - base_step).abs() <= 1e-14 {
+                f64::NAN
+            } else {
+                probe[0] * probe[0]
+            }
+        },
+        &mut gradient,
+    );
+
+    assert_near(gradient[0], 2.0, 1e-10);
+}
+
+#[test]
+fn central_diff_hessian_retries_step_when_primary_step_is_invalid() {
+    let param = [1.0_f64];
+    let rel_step = 1e-3;
+    let min_step = 1e-3;
+    let base_step = ((param[0].abs() + 1.0) * rel_step).max(min_step);
+    let hessian = models::central_diff_hessian_from_gradient(
+        &param,
+        rel_step,
+        min_step,
+        |probe, gradient_out| {
+            let delta = (probe[0] - param[0]).abs();
+            gradient_out[0] = if (delta - base_step).abs() <= 1e-14 {
+                f64::NAN
+            } else {
+                2.0 * probe[0]
+            };
+        },
+    );
+
+    assert_near(hessian[[0, 0]], 2.0, 1e-10);
+}
+
+#[test]
+fn fit_numerical_hessian_retries_step_when_primary_step_is_invalid() {
+    let param = Array1::from_vec(vec![1.0_f64]);
+    let base_step = ((param[0].abs() + 1.0) * HESSIAN_FD_REL_STEP).max(HESSIAN_FD_MIN_STEP);
+    let problem = RetryGradientProblem {
+        center: param[0],
+        invalid_step: base_step,
+    };
+    let hessian =
+        numerical_hessian_from_gradient(&problem, &param).expect("hessian must be computed");
+
+    assert_near(hessian[[0, 0]], 2.0 + HESSIAN_DIAGONAL_JITTER, 1e-10);
+}
+
+#[test]
+fn central_diff_gradient_falls_back_to_zero_when_all_retry_steps_invalid() {
+    let param = [1.0_f64];
+    let mut gradient = [123.0];
+    models::central_diff_gradient_from_value(&param, 1e-3, 1e-3, |_probe| f64::NAN, &mut gradient);
+
+    assert_eq!(gradient[0], 0.0);
+}
+
+#[test]
+fn fit_numerical_hessian_falls_back_to_diagonal_jitter_when_all_retry_steps_invalid() {
+    let param = Array1::from_vec(vec![1.0_f64]);
+    let hessian = numerical_hessian_from_gradient(&AlwaysInvalidGradientProblem, &param)
+        .expect("hessian must be computed even with invalid gradients");
+
+    assert_near(hessian[[0, 0]], HESSIAN_DIAGONAL_JITTER, 1e-15);
+}
+
 #[test]
 fn curve_objective_arrhenius_is_consistent_across_levels() {
     let x_values = [0.4, 0.8, 1.4, 2.5, 4.0];