Anti-Degeneracy Defense System¶

This section provides comprehensive documentation for the Anti-Degeneracy Defense System introduced in Homodyne and enhanced with quantile-based per-angle scaling. The system addresses fundamental optimization challenges in laminar flow XPCS analysis where structural degeneracy and gradient cancellation prevent accurate parameter estimation.

Introduction¶

When fitting XPCS data under laminar flow conditions with many azimuthal angles, the optimizer often fails to correctly estimate shear parameters (particularly gamma_dot_t0). This “parameter collapse” occurs because the optimization landscape contains degenerate solutions where per-angle scaling parameters absorb physical signals.

Example Problem Scenario:

Dataset: 23 phi angles, 1M+ points per angle
Expected: gamma_dot_t0 ≈ 10^-3 (true physical value)
Observed: gamma_dot_t0 → 10^-6 (collapsed to lower bound)
Root cause: Per-angle contrast/offset parameters absorb angle-dependent shear signal

The Anti-Degeneracy Defense System provides a five-layer solution:

Defense Layers Summary¶
Layer	Solution	Root Cause Addressed	Status	Effectiveness
1	Constant/Individual Mode Selection	Structural Degeneracy	ENHANCED	High
2	Hierarchical Optimization	Gradient Cancellation	Local NLSQ only	High
3	Adaptive CV Regularization	Ineffective λ	Local NLSQ only	Medium
4	Gradient Collapse Detection	Runtime Monitoring	Local NLSQ only	Medium
5	Shear-Sensitivity Weighting	Gradient Cancellation	Local NLSQ only	High
—	Data Shuffling	Sequential Bias	EXISTING	Medium

Note

Cross-Path Uniformity: Layer 1 (Constant/Individual Mode) now applies uniformly across all optimization paths: local NLSQ, multi-start, and CMA-ES. Layers 2-5 are specific to local NLSQ where gradient-based optimizations benefit from hierarchical alternation, regularization, and gradient monitoring.

Theoretical Background¶

Understanding the Root Causes¶

To understand why the defense system is necessary, we must examine the three fundamental root causes of parameter collapse in laminar flow fitting.

Root Cause 1: Gradient Cancellation¶

The shear term in the XPCS model introduces angle-dependent dynamics:

\[c_1(\phi, t_1, t_2) = \exp\left(-q^2 J(t_1, t_2)\right) \times \text{sinc}\left[\frac{qh}{2\pi} \Gamma(t_1, t_2) \cos(\phi - \phi_0)\right]\]

where \(\Gamma(t_1, t_2)\) is the cumulative shear integral proportional to gamma_dot_t0.

The gradient of the loss with respect to gamma_dot_t0 is:

\[\frac{\partial L}{\partial \dot{\gamma}_0} \propto \sum_{\phi} \cos(\phi_0 - \phi) \times (\text{data terms})\]

The Problem: The term \(\cos(\phi_0 - \phi)\) changes sign across different angles. When summed globally (as in standard gradient descent), contributions from different angles partially cancel, producing a weak net gradient.

Example with 8 equally-spaced angles (0°, 45°, 90°, ...):

φ = 0°:   cos(0°)   = +1.0    ─┐
φ = 45°:  cos(45°)  = +0.71   ─┼─ These partially cancel when summed
φ = 90°:  cos(90°)  = 0.0     ─┤
φ = 135°: cos(135°) = -0.71   ─┤
φ = 180°: cos(180°) = -1.0    ─┘

Net effect: Weak gradient signal for gamma_dot_t0

Root Cause 2: Structural Degeneracy¶

With n_phi azimuthal angles, the per-angle scaling introduces 2 × n_phi free parameters (contrast and offset for each angle). This creates an over-determined fitting problem:

\[c_2^{\text{fitted}}(\phi_i) = \text{offset}_i + \text{contrast}_i \times [c_1(\phi_i)]^2\]

The Problem: Per-angle parameters can absorb nearly any angle-dependent signal. When gamma_dot_t0 produces angle-dependent \(c_1\) values, the optimizer finds it easier to adjust 2 × 23 = 46 per-angle parameters than to correctly estimate 7 physical parameters.

Gradient magnitude comparison (typical 23-angle fit):

|∇_{contrast}| ≈ 10⁻²  ─┐
|∇_{offset}|   ≈ 10⁻²  ─┼─ ~92% of gradient magnitude
|∇_{physical}| ≈ 10⁻⁴  ─┘   ~8% of gradient magnitude

Result: Optimizer preferentially adjusts per-angle params

Root Cause 3: Ineffective Regularization¶

Previous versions used group variance regularization with λ = 0.01:

\[L_{\text{reg}} = \lambda \times \left[\text{Var}(\text{contrast}) + \text{Var}(\text{offset})\right] \times N\]

With typical MSE values and data sizes, this contributed only ~0.05% to the total loss—providing no effective constraint on per-angle parameter variation.

The Math:

Typical values:
- MSE ≈ 10⁻⁴
- N_points = 23M
- Var(contrast) ≈ 0.01 (10% std on mean of 0.5)

L_data = MSE × N = 10⁻⁴ × 23×10⁶ = 2.3×10³
L_reg  = 0.01 × 0.01 × 23×10⁶ = 2.3×10³

Contribution = L_reg / L_total ≈ 0.0001 (0.01%)

The regularization is effectively invisible to the optimizer!

Layer 1: Fourier Reparameterization¶

Mathematical Foundation¶

Instead of treating per-angle contrast and offset as independent parameters, Layer 1 expresses them as truncated Fourier series:

\[\text{contrast}(\phi) = c_0 + \sum_{k=1}^{K} \left[c_k \cos(k\phi) + s_k \sin(k\phi)\right]\]

\[\text{offset}(\phi) = o_0 + \sum_{k=1}^{K} \left[o_k \cos(k\phi) + t_k \sin(k\phi)\right]\]

For Fourier order K=2:

Contrast: 5 coefficients [c₀, c₁, s₁, c₂, s₂]
Offset: 5 coefficients [o₀, o₁, t₁, o₂, t₂]
Total: 10 Fourier coefficients vs 2×n_phi independent parameters

Per-Angle Mode Selection Guide¶

The per_angle_mode setting controls how per-angle parameters (contrast/offset) are handled. This is the most important configuration choice for laminar flow fitting with many angles.

per_angle_mode Options¶
Mode	Description
`auto`	DEFAULT For n_phi >= 3: quantile-estimated per-angle scaling is averaged and optimized as 9 parameters (7 physical + 2 averaged scaling). For n_phi < 3: falls back to `individual`. Most robust against degeneracy.
`constant`	Per-angle scaling computed from quantile estimation and fixed (not optimized). Result: 7 physical parameters only. Fastest convergence but no scaling optimization.
`individual`	Each angle has independent contrast and offset (2 × n_phi + 7 optimized params). Full flexibility but high degeneracy risk for n_phi > 6.
`fourier`	Contrast/offset expressed as truncated Fourier series (2 × (2K+1) + 7 optimized params). Enforces smooth angular variation, reduces degeneracy risk.

Note

The default per_angle_mode is "auto", which for datasets with 3+ phi angles uses averaged scaling (auto_averaged internally): quantile-estimated contrast/offset are averaged and optimized alongside 7 physical parameters (9 total). Explicit "constant" mode (fixed_constant internally) instead fixes scaling from quantiles without optimizing it (7 params only).

Mode Comparison¶

per_angle_mode Comparison¶
Aspect	auto	constant	individual	fourier	Recommendation	Notes
Optimized Params	9 (N≥3) or 2N+7 (N<3)	7 (scaling fixed)	2N + 7	10-17 + 7	auto/constant	Lower = faster, more robust
Degeneracy Risk	Low	Low	High (N > 6)	Medium	auto/constant	Low = cannot easily absorb shear signal
Flexibility	Adaptive	Low	Maximum	Medium	depends	High flex = can fit noise
Speed	Fast	Fastest	Slowest	Medium	constant	Fewer params = faster convergence
Angular Variation	Uniform (averaged)	Fixed per-angle (not optimized)	Fully flexible	Smooth only	depends	Physical basis matters

Decision Flowchart¶

Use this flowchart to select the appropriate per_angle_mode:

                    ┌─────────────────────────┐
                    │ Start: laminar_flow     │
                    │ mode with N phi angles  │
                    └───────────┬─────────────┘
                                │
                    ┌───────────▼───────────┐
                    │ Use "auto" (DEFAULT)  │
                    │ N >= 3 → averaged     │
                    │ N < 3  → individual   │
                    └───────────┬───────────┘
                                │
          ┌─────────────────────┴─────────────────────┐
          │ Need explicit control?                    │
          ▼                                           ▼
┌─────────────────────┐               ┌───────────────────────────┐
│ NO: Use "auto"      │               │ YES: Choose explicit mode │
│ (RECOMMENDED)       │               │                           │
│                     │               └─────────────┬─────────────┘
│ • Adaptive to N     │                             │
│ • 9 params for N≥3  │               ┌─────────────┴─────────────┐
│ • Most robust       │               │                           │
└─────────────────────┘               ▼                           ▼
                          ┌─────────────────────┐   ┌─────────────────────┐
                          │ Smooth angular      │   │ Full flexibility    │
                          │ variation expected? │   │ needed? (N ≤ 3)     │
                          │ (sample asymmetry)  │   │                     │
                          └─────────┬───────────┘   └─────────┬───────────┘
                                    │                         │
                          ┌─────────▼─────────────┐ ┌─────────▼───────────┐
                          │ Use "fourier"         │ │ Use "individual"    │
                          │                       │ │                     │
                          │ • Smooth variation    │ │ • Full flexibility  │
                          │ • N > 6 recommended   │ │ • N ≤ 3 only        │
                          │ • 17 params (K=2)     │ │ • High degeneracy   │
                          └───────────────────────┘ └─────────────────────┘

Parameter Count Reduction¶

The different modes provide varying degrees of parameter reduction:

Parameter Count Comparison (laminar_flow, 7 physical params)¶
n_phi	Individual	Fourier (K=2)	Constant (fixed)	Auto Selection	Auto Params
2	4 + 7 = 11	4 + 7 = 11*	0 + 7 = 7	individual (N<3)	11
3	6 + 7 = 13	6 + 7 = 13*	0 + 7 = 7	averaged (N>=3)	9
6	12 + 7 = 19	10 + 7 = 17	0 + 7 = 7	averaged	9
10	20 + 7 = 27	10 + 7 = 17	0 + 7 = 7	averaged	9
23	46 + 7 = 53	10 + 7 = 17	0 + 7 = 7	averaged	9
100	200 + 7 = 207	10 + 7 = 17	0 + 7 = 7	averaged	9

* For n_phi <= 2×(order+1), Fourier mode provides no reduction over individual
Constant (fixed): Scaling fixed from quantile estimation, NOT optimized (7 params)
Auto averaged: 2 averaged scaling params (contrast + offset) ARE optimized (9 params)

Quantile-Based Scaling Estimation¶

Both auto and constant modes use quantile-based estimation to compute per-angle contrast/offset from data. The key difference is how these estimates are used:

Auto mode (auto_averaged internally): Averages per-angle estimates to single values and optimizes them alongside 7 physical parameters (9 total).
Constant mode (fixed_constant internally): Uses per-angle estimates as fixed values (not optimized), fitting only 7 physical parameters.

Physics Foundation:

The Siegert relation connects intensity autocorrelation (C2) to field autocorrelation (g1):

\[C_2 = \text{contrast} \times g_1^2 + \text{offset}\]

The key insight is that g1 decays with time lag:

At small time lags (Δt → 0): g1² ≈ 1, so C2 ≈ contrast + offset (the “ceiling”)
At large time lags (Δt → ∞): g1² → 0, so C2 → offset (the “floor”)

Quantile-Based Estimation Algorithm:

For each phi angle independently:

Step 1: Partition data by time lag
┌──────────────────────────────────────────────────────────────┐
│ Time Lag Distribution for Angle φᵢ                          │
│                                                              │
│    ▼ Small lags (bottom 20%)    ▼ Large lags (top 20%)      │
│ ├─────────────────────────────────────────────────────────┤ │
│ └──── Used for ceiling ────┘       └──── Used for floor ───┘│
└──────────────────────────────────────────────────────────────┘

Step 2: Estimate OFFSET (from large-lag region where g1² → 0)
offset = 10th percentile of C2 values at large lags

Step 3: Estimate CONTRAST (from small-lag region where g1² ≈ 1)
ceiling = 90th percentile of C2 values at small lags
contrast = ceiling - offset

Step 4: Apply bounds clipping
offset = clip(offset, offset_bounds)
contrast = clip(contrast, contrast_bounds)

Step 5: AVERAGE across all phi angles
avg_contrast = mean(contrast_per_angle)
avg_offset = mean(offset_per_angle)

Why Quantiles Instead of Min/Max:

Outlier Robustness: Quantiles ignore extreme values from noise spikes
Noise Tolerance: Less sensitive to measurement noise
Systematic Error Handling: Avoids contamination from partially decayed points

Default Quantile Parameters:

Parameter	Default	Description
`lag_floor_quantile`	0.80	Top 20% of lags for floor estimation
`lag_ceiling_quantile`	0.20	Bottom 20% of lags for ceiling estimation
`value_quantile_low`	0.10	10th percentile for robust floor
`value_quantile_high`	0.90	90th percentile for robust ceiling

9-Parameter Optimization:

AUTO (AVERAGED) MODE 9-PARAMETER OPTIMIZATION
=========================================================

Before optimization:
    1. Compute per-angle scaling from data quantiles
       contrast_per_angle[n_phi] ← quantile estimation
       offset_per_angle[n_phi]   ← quantile estimation

    2. AVERAGE to get initial estimates
       avg_contrast = mean(contrast_per_angle)
       avg_offset = mean(offset_per_angle)

During optimization:
    Parameter vector: [avg_contrast, avg_offset, D₀, α, D_offset, γ̇₀, β, γ̇_offset, φ₀]
                      (9 params: 2 averaged scaling + 7 physical)

    Model evaluation:
       contrast = broadcast(avg_contrast, n_phi)  # Same for all angles
       offset = broadcast(avg_offset, n_phi)      # Same for all angles

After optimization:
    Expand for output: [contrast..., offset..., physical_opt]
                       (2*n_phi + 7 values for backward compatibility)

CONSTANT (FIXED) MODE 7-PARAMETER OPTIMIZATION
==========================================================

Before optimization:
    1. Compute per-angle scaling from data quantiles
       contrast_per_angle[n_phi] ← quantile estimation
       offset_per_angle[n_phi]   ← quantile estimation

    2. FIX per-angle values (not optimized)

During optimization:
    Parameter vector: [D₀, α, D_offset, γ̇₀, β, γ̇_offset, φ₀]
                      (7 physical params only)

    Model evaluation:
       contrast = fixed_contrast_per_angle[n_phi]  # Fixed per angle
       offset = fixed_offset_per_angle[n_phi]      # Fixed per angle

After optimization:
    Expand for output: [fixed_contrast..., fixed_offset..., physical_opt]
                       (2*n_phi + 7 values for backward compatibility)

When to Use Auto Mode (Default):

Most XPCS datasets where detector response is reasonably uniform
Large datasets (>1M points) where parameter reduction improves convergence
Multi-start optimization with many angles (tractable parameter count)
When per-angle variation is primarily noise, not physics

When to Use Constant Mode:

When you want zero per-angle optimization (fastest, 7 physical params only)
When quantile estimates are trusted and no further tuning is needed

When to Switch to Other Modes:

Use fourier if smooth angular variation is physically expected (sample asymmetry)
Use individual for very small n_phi (≤ 3) or when full flexibility is needed

Configuration:

optimization:
  nlsq:
    anti_degeneracy:
      enable: true
      per_angle_mode: "auto"          # Default: averaged scaling, 9 params (N≥3)
      # per_angle_mode: "constant"    # Fixed scaling from quantiles, 7 params
      # per_angle_mode: "fourier"     # Truncated Fourier series
      # per_angle_mode: "individual"  # Independent per-angle params

Impact on Multi-Start and Global Optimization:

Auto/constant modes make multi-start and CMA-ES optimization tractable for many-angle datasets:

Multi-Start/CMA-ES Tractability (23-angle laminar_flow)¶
Mode	n_params	Min n_starts	Recommended n_starts
individual	53	53	100-150 (expensive)
fourier	17	17	20-40 (moderate)
constant (via auto)	9	9	10-15 (tractable)

Note

Uniform Anti-Degeneracy: Both local NLSQ and global optimization (CMA-ES, multi-start) now use the same 9-parameter constant mode. This ensures consistent behavior across all optimization paths.

Why Constant Mode Works Best¶

Significant Parameter Reduction: Only 9 parameters (7 physical + 2 averaged scaling) are optimized instead of 53 (23×2 per-angle + 7 physical) for a 23-angle dataset.
Low Degeneracy Risk: With only 2 averaged scaling parameters (shared across all angles), they cannot absorb angle-dependent shear signals that vary with φ.
Data-Driven Initialization: Quantile-based estimation uses actual experimental data, not arbitrary defaults, for contrast/offset values.
Robust to Noise: Quantile-based estimation is more robust to outliers than least-squares or min/max approaches.
Uniform Across Paths: Both local NLSQ, multi-start, and CMA-ES global optimization use the same 9-parameter constant mode.

Why Fourier Mode is Still Useful¶

Reduced Degrees of Freedom: 10 Fourier coefficients cannot absorb arbitrary angle-dependent signals that would require 46 free parameters.
Smooth Variation: Fourier series produces physically reasonable smooth variation in contrast/offset across angles, matching experimental reality.
Preserved Information: The dominant modes (DC + first two harmonics) capture physically meaningful variations (e.g., sample asymmetry, beam position offsets).
Breaks Degeneracy: With only 10 coefficients for per-angle scaling, the optimizer must rely on physical parameters to explain angle-dependent shear dynamics.

Implementation¶

The FourierReparameterizer class handles all conversions:

from homodyne.optimization.nlsq.fourier_reparam import (
    FourierReparameterizer,
    FourierReparamConfig,
)

# Configure Fourier reparameterization
config = FourierReparamConfig(
    mode="auto",         # "independent", "fourier", or "auto"
    fourier_order=2,     # Number of Fourier harmonics
    auto_threshold=6,    # Use Fourier when n_phi > threshold
)

# Create reparameterizer for your phi angles
phi_angles = np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5])  # radians
reparameterizer = FourierReparameterizer(phi_angles, config)

# Convert Fourier coefficients to per-angle values
fourier_coeffs = np.array([0.5, 0.01, -0.02, 0.005, 0.003,  # contrast
                           1.0, 0.02, 0.01, -0.01, 0.005])  # offset
contrast, offset = reparameterizer.fourier_to_per_angle(fourier_coeffs)

# Convert back (useful for initialization)
recovered_coeffs = reparameterizer.per_angle_to_fourier(contrast, offset)

Basis Matrix Computation¶

The Fourier basis matrix B satisfies values = B @ coefficients:

\[\begin{split}B = \begin{bmatrix} 1 & \cos(\phi_1) & \sin(\phi_1) & \cos(2\phi_1) & \sin(2\phi_1) \\ 1 & \cos(\phi_2) & \sin(\phi_2) & \cos(2\phi_2) & \sin(2\phi_2) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \cos(\phi_n) & \sin(\phi_n) & \cos(2\phi_n) & \sin(2\phi_n) \end{bmatrix}\end{split}\]

Covariance Transformation¶

After optimization, uncertainties must be transformed from Fourier space to per-angle space. The Jacobian of the transformation is used:

\[\Sigma_{\text{per-angle}} = J \times \Sigma_{\text{Fourier}} \times J^T\]

where \(J = \frac{\partial (\text{per-angle})}{\partial (\text{Fourier})}\).

ParameterIndexMapper Reference¶

The ParameterIndexMapper class provides consistent parameter index calculations across all modes. This is the single source of truth for parameter group boundaries.

from homodyne.optimization.nlsq.parameter_index_mapper import ParameterIndexMapper

# Constant mode (23 phi angles, 7 physical params)
mapper = ParameterIndexMapper(n_phi=23, n_physical=7, use_constant=True)
print(mapper.mode_name)           # "constant"
print(mapper.n_per_angle_total)   # 2 (single contrast + offset, shared)
print(mapper.total_params)        # 9 (2 + 7)
print(mapper.get_group_indices()) # [(0, 1), (1, 2)]

# Individual mode (23 phi angles, 7 physical params)
mapper = ParameterIndexMapper(n_phi=23, n_physical=7, use_constant=False)
print(mapper.mode_name)           # "individual"
print(mapper.n_per_angle_total)   # 46 (23 contrasts + 23 offsets)
print(mapper.total_params)        # 53 (46 + 7)
print(mapper.get_group_indices()) # [(0, 23), (23, 46)]

# Fourier mode (23 phi angles, order=2, 7 physical params)
mapper = ParameterIndexMapper(n_phi=23, n_physical=7, fourier=fourier_obj)
print(mapper.mode_name)           # "fourier"
print(mapper.n_per_angle_total)   # 10 (5 contrast + 5 offset coeffs)
print(mapper.total_params)        # 17 (10 + 7)
print(mapper.get_group_indices()) # [(0, 5), (5, 10)]

ParameterIndexMapper by Mode (n_phi=23, n_physical=7)¶
Mode	n_per_group	n_per_angle_total	get_group_indices()	total_params
`constant`	1	2	[(0, 1), (1, 2)]	9
`individual`	23	46	[(0, 23), (23, 46)]	53
`fourier` (K=2)	5	10	[(0, 5), (5, 10)]	17

Layer 2: Hierarchical Optimization¶

Algorithmic Approach¶

Hierarchical optimization breaks the gradient cancellation by alternating between two optimization stages:

Stage 1: Physical Parameters Only

Freeze per-angle parameters at current values
Optimize only physical parameters (D₀, α, γ̇₀, β, etc.)
Physical gradients are not diluted by per-angle updates

Stage 2: Per-Angle Parameters Only

Freeze physical parameters at current values
Optimize only per-angle scaling (contrast, offset)
Per-angle params adjust to match the fixed physics model

This alternation continues until convergence or maximum iterations.

HIERARCHICAL OPTIMIZATION ALGORITHM
====================================

Initialize: params = [per_angle_params, physical_params]

for outer_iter in range(max_outer_iterations):

    ┌──────────────────────────────────────────────────┐
    │ STAGE 1: Fit PHYSICAL params (per-angle frozen) │
    │                                                  │
    │ • Per-angle gradients = 0 (frozen)              │
    │ • Physical gradients get 100% of update         │
    │ • Uses L-BFGS with physical-only bounds         │
    └──────────────────────────────────────────────────┘
                           │
                           ▼
    ┌──────────────────────────────────────────────────┐
    │ STAGE 2: Fit PER-ANGLE params (physical frozen) │
    │                                                  │
    │ • Physical gradients = 0 (frozen)               │
    │ • Per-angle params adjust to fixed physics      │
    │ • Uses L-BFGS with per-angle-only bounds        │
    └──────────────────────────────────────────────────┘
                           │
                           ▼
    ┌──────────────────────────────────────────────────┐
    │ CHECK CONVERGENCE                               │
    │                                                  │
    │ if ||physical_new - physical_old|| < tolerance: │
    │     break                                        │
    └──────────────────────────────────────────────────┘

Return: optimized [per_angle_params, physical_params]

Why This Works¶

No Gradient Competition: In Stage 1, physical parameters receive 100% of the gradient update without competition from per-angle parameters.
Proper Attribution: In Stage 2, per-angle parameters must adjust to match a fixed physics model, correctly attributing angle-dependent variations.
Iterative Refinement: Multiple alternations allow both parameter groups to converge to a consistent solution.
Natural Regularization: The alternating freeze prevents either parameter group from dominating the optimization.

Implementation¶

from homodyne.optimization.nlsq.hierarchical import (
    HierarchicalOptimizer,
    HierarchicalConfig,
    HierarchicalResult,
)

# Configure hierarchical optimization
config = HierarchicalConfig(
    enable=True,
    max_outer_iterations=5,
    outer_tolerance=1e-6,
    physical_max_iterations=100,
    physical_ftol=1e-8,
    per_angle_max_iterations=50,
    per_angle_ftol=1e-6,
    log_stage_transitions=True,
)

# Create optimizer
optimizer = HierarchicalOptimizer(
    config=config,
    n_phi=23,
    n_physical=7,
    fourier_reparameterizer=fourier,  # Optional, for Layer 1 integration
)

# Run optimization
result = optimizer.fit(
    loss_fn=loss_function,
    grad_fn=gradient_function,
    p0=initial_params,
    bounds=(lower_bounds, upper_bounds),
)

# Access results
print(f"Converged in {result.n_outer_iterations} outer iterations")
print(f"Final loss: {result.fun:.6g}")

# Examine optimization history
for entry in result.history:
    print(f"Iter {entry['outer_iter']}: Stage1 loss={entry['stage1_loss']:.4g}, "
          f"Stage2 loss={entry['stage2_loss']:.4g}")

Layer 3: Adaptive CV Regularization¶

From Absolute to Relative Regularization¶

Previous regularization used absolute variance:

\[L_{\text{reg,old}} = \lambda \times \text{Var}(\text{params}) \times N\]

The new approach uses coefficient of variation (CV) with relative scaling:

\[CV = \frac{\text{std}(\text{params})}{|\text{mean}(\text{params})|}\]

\[L_{\text{reg,new}} = \lambda \times CV^2 \times \text{MSE} \times N\]

Auto-Tuned Lambda¶

The regularization strength is automatically computed to achieve target contribution to the loss:

\[\lambda = \frac{\text{target\_contribution}}{\text{target\_CV}^2}\]

Example:

Target: Allow 10% variation (CV = 0.1)
Target contribution: 10% of MSE
Auto-computed: λ = 0.1 / 0.01 = 10

This is 100× stronger than the previous λ = 0.01, providing effective constraint on per-angle parameter variation.

Why CV-Based Regularization¶

Scale-Invariant: CV measures relative variation, not absolute. A 5% variation is penalized the same whether the mean is 0.5 or 5.0.
Physically Meaningful: CV = 0.1 means 10% standard deviation relative to mean—easy to interpret and tune.
Properly Scaled: By scaling with MSE, the regularization maintains consistent importance regardless of data noise level.
Adaptive: Auto-tuning ensures regularization contributes meaningfully to the optimization objective.

Implementation¶

from homodyne.optimization.nlsq.adaptive_regularization import (
    AdaptiveRegularizer,
    AdaptiveRegularizationConfig,
)

# Configure adaptive regularization
config = AdaptiveRegularizationConfig(
    enable=True,
    mode="relative",            # "absolute", "relative", or "auto"
    lambda_base=1.0,            # Base lambda (used if auto_tune=False)
    target_cv=0.10,             # Target coefficient of variation
    target_contribution=0.10,   # Target fraction of MSE
    auto_tune_lambda=True,      # Compute lambda from targets
    max_cv=0.20,                # Maximum allowed CV
)

# Create regularizer
regularizer = AdaptiveRegularizer(config, n_phi=23)

# Compute regularization in loss function
def regularized_loss(params, data):
    mse = compute_mse(params, data)
    n_points = len(data)
    reg_term = regularizer.compute_regularization(params, mse, n_points)
    return mse * n_points + reg_term

# Check for constraint violations
violations = regularizer.check_constraint_violation(params)
if violations:
    print(f"Warning: CV exceeds max_cv threshold: {violations}")

Layer 4: Gradient Collapse Detection¶

Runtime Monitoring¶

Layer 4 monitors the optimization in real-time to detect gradient collapse— when physical parameter gradients become negligible compared to per-angle gradients.

Detection Criterion:

\[\text{ratio} = \frac{|\nabla_{\text{physical}}|}{|\nabla_{\text{per-angle}}|}\]

When ratio < threshold (default: 0.01) for N consecutive iterations (default: 5), gradient collapse is declared.

Response Actions¶

Upon detecting gradient collapse, the system can:

warn: Log warning, continue optimization
hierarchical: Switch to hierarchical optimization (Layer 2)
reset: Reset per-angle parameters to their mean values
abort: Terminate optimization with warning

Implementation¶

from homodyne.optimization.nlsq.gradient_monitor import (
    GradientCollapseMonitor,
    GradientMonitorConfig,
)

# Configure gradient monitoring
config = GradientMonitorConfig(
    enable=True,
    ratio_threshold=0.01,       # Collapse threshold
    consecutive_triggers=5,     # Must trigger N times
    response_mode="hierarchical",  # Action on collapse
    reset_per_angle_to_mean=True,
    lambda_multiplier_on_collapse=10.0,
)

# Create monitor
monitor = GradientCollapseMonitor(
    config=config,
    physical_indices=[10, 11, 12, 13, 14, 15, 16],  # Physical param indices
    per_angle_indices=list(range(10)),  # Per-angle param indices
)

# Check gradients during optimization
for iteration in range(max_iterations):
    gradients = compute_gradients(params)

    status = monitor.check(gradients, iteration, params)

    if status == "COLLAPSE_DETECTED":
        response = monitor.get_response()
        if response["mode"] == "hierarchical":
            # Switch to hierarchical optimization
            result = run_hierarchical_optimization(params)
            break

    # Normal optimization step
    params = update_params(params, gradients)

# Get final diagnostics
diagnostics = monitor.get_diagnostics()
print(f"Gradient monitoring: {diagnostics}")

Layer 5: Shear-Sensitivity Weighting¶

Mathematical Foundation¶

Layer 5 addresses gradient cancellation by weighting residuals according to their sensitivity to the shear parameter gamma_dot_t0. The shear term in the XPCS model produces gradients that scale with \(|\cos(\phi_0 - \phi)|\):

Angles where \(\phi \approx \phi_0\) or \(\phi \approx \phi_0 + \pi\) are highly sensitive to shear (parallel/anti-parallel flow)
Angles where \(\phi \approx \phi_0 \pm \pi/2\) are insensitive to shear (perpendicular flow)

Weight Formula:

\[w(\phi) = w_{\min} + (1 - w_{\min}) \times |\cos(\phi_0 - \phi)|^\alpha\]

where:

\(w_{\min}\) is the minimum weight (default: 0.3) ensuring perpendicular angles still contribute
\(\alpha\) is the sensitivity exponent (default: 1.0 for linear weighting)

The weights are normalized so their mean equals 1.0, preserving the overall loss scale.

Why This Works¶

Standard unweighted least squares gives equal importance to all angles. When angles span 360°, the positive and negative contributions to \(\partial L / \partial \dot{\gamma}_0\) partially cancel (94.6% cancellation for 23 uniformly-spaced angles).

Shear-sensitivity weighting addresses this by:

Amplifying informative residuals: Errors at shear-sensitive angles (parallel flow) contribute more to the loss gradient
Attenuating uninformative residuals: Errors at insensitive angles (perpendicular flow) contribute less
Breaking symmetry: The asymmetric weighting prevents gradient cancellation

Weight distribution for 8 angles with φ₀ = 0°:

φ = 0°:   |cos(0°)|   = 1.0  → w = 1.0   (HIGH weight)
φ = 45°:  |cos(45°)|  = 0.71 → w = 0.80
φ = 90°:  |cos(90°)|  = 0.0  → w = 0.30  (LOW weight)
φ = 135°: |cos(135°)| = 0.71 → w = 0.80
φ = 180°: |cos(180°)| = 1.0  → w = 1.0   (HIGH weight)
...

Net effect: Gradient from parallel angles dominates,
            preventing cancellation

Data Flow¶

Shear-sensitivity weighting is computed in Homodyne and passed to NLSQ as generic residual weights, maintaining separation of concerns:

┌─────────────────────────────────────────────────────────────────────────┐
│ Homodyne wrapper.py                                                     │
│  ┌──────────────────────┐    ┌─────────────────────────────────────┐   │
│  │ ShearSensitivity     │───▶│ get_weights() → numpy array (n_phi) │   │
│  │ Weighting            │    └─────────────────────────────────────┘   │
│  └──────────────────────┘                      │                        │
│                                                ▼                        │
│                                    .tolist() → Python list              │
│                                                │                        │
│  ┌─────────────────────────────────────────────▼────────────────────┐  │
│  │ HybridStreamingConfig(                                           │  │
│  │   enable_residual_weighting=True,                                │  │
│  │   residual_weights=[w0, w1, ..., w_{n_phi-1}]                    │  │
│  │ )                                                                │  │
│  └──────────────────────────────────────────────────────────────────┘  │
└─────────────────────────────────────────────────────────────────────────┘
                                    │
                                    ▼
┌─────────────────────────────────────────────────────────────────────────┐
│ NLSQ adaptive_hybrid.py                                                 │
│  ┌─────────────────────┐                                               │
│  │ _setup_residual_    │ Converts list → jnp.array(float64)            │
│  │ weights()           │                                               │
│  └─────────────────────┘                                               │
│           │                                                             │
│           ▼                                                             │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ loss_fn(params, x_batch, y_batch):                              │   │
│  │   group_idx = x_batch[:, 0].astype(int32)  # φ index per point  │   │
│  │   weights = residual_weights[group_idx]    # lookup per-point   │   │
│  │   wmse = sum(w * r²) / sum(w)              # weighted MSE       │   │
│  └─────────────────────────────────────────────────────────────────┘   │
└─────────────────────────────────────────────────────────────────────────┘

Dynamic Weight Updates¶

As the optimizer converges, \(\phi_0\) may change from its initial guess. The ShearSensitivityWeighter class supports dynamic weight updates:

from homodyne.optimization.nlsq.shear_weighting import (
    ShearSensitivityWeighter,
    ShearWeightingConfig,
)

# Create weighter with initial phi0 estimate
config = ShearWeightingConfig(
    enable=True,
    min_weight=0.3,
    alpha=1.0,
    update_frequency=1,  # Update every outer iteration
    normalize=True,
)
weighter = ShearSensitivityWeighter(phi_angles, config, phi0_initial=0.0)

# Get initial weights
weights = weighter.get_weights()
print(f"Initial weights: {weights}")

# During optimization, update phi0 as it converges
for outer_iter in range(max_iterations):
    # ... optimization step ...

    if outer_iter % config.update_frequency == 0:
        weighter.update_phi0(current_phi0_estimate)
        weights = weighter.get_weights()  # Weights recomputed

# Get diagnostics
diagnostics = weighter.get_diagnostics()
print(f"Weight range: [{diagnostics['weight_min']:.3f}, {diagnostics['weight_max']:.3f}]")

Implementation¶

from homodyne.optimization.nlsq.shear_weighting import (
    ShearSensitivityWeighter,
    ShearWeightingConfig,
    create_shear_weighting,
)

# High-level factory function (recommended)
phi_angles = np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5])  # radians
weighter = create_shear_weighting(
    phi_angles=phi_angles,
    phi0=0.0,
    config={
        'enable': True,
        'min_weight': 0.3,
        'alpha': 1.0,
    }
)

if weighter is not None:
    weights = weighter.get_weights()
    print(f"n_weights: {len(weights)}")
    print(f"Weight range: [{min(weights):.3f}, {max(weights):.3f}]")

Cross-Path Uniformity¶

The Anti-Degeneracy System applies across three optimization paths with different layer coverage:

Layer Coverage by Optimization Path¶
Optimization Path	Layer 1	Layer 2	Layer 3	Layer 4	Layer 5
Local NLSQ (NLSQWrapper)	✅	✅	✅	✅	✅
Multi-Start (fit_nlsq_multistart)	✅	✅	✅	✅	✅
CMA-ES (fit_nlsq_cmaes)	✅	—	—	—	—

Why Layers 2-5 Are Local NLSQ Only:

CMA-ES is a population-based evolutionary strategy that operates fundamentally differently from gradient-based local optimization:

Layer 2 (Hierarchical): Designed for alternating physical/per-angle optimization stages. CMA-ES already evolves all parameters simultaneously via population.
Layer 3 (Adaptive Regularization): Adds variance penalty to gradient-based loss. CMA-ES uses fitness ranking, not gradient-based updates.
Layer 4 (Gradient Monitor): Monitors gradient magnitude ratios. CMA-ES doesn’t compute gradients—it uses evolutionary selection.
Layer 5 (Shear Weighting): Weights residuals by abs(cos(φ₀-φ)). Could be added to CMA-ES fitness evaluation but not currently implemented.

Layer 1 Is Sufficient for CMA-ES:

Layer 1 (Constant Mode) addresses the primary degeneracy issue by:

Reducing 46 per-angle parameters to 2 averaged parameters
Eliminating the ability for per-angle params to absorb physical signals
Making the 9-parameter space tractable for population-based search

Multi-Start Inherits All Layers:

Multi-start optimization calls fit_nlsq_jax() for each starting point, which routes through NLSQWrapper. Therefore, multi-start inherits all 5 layers from the local NLSQ path.

Configuration Reference¶

YAML Configuration¶

Complete YAML configuration for the Anti-Degeneracy Defense System:

optimization:
  nlsq:
    # ... existing settings ...

    # Anti-Degeneracy Defense System
    anti_degeneracy:

      # Layer 1: Per-Angle Mode Selection (defaults)
      per_angle_mode: "auto"            # DEFAULT: selects constant (N>=3) or individual (N<3)
      # per_angle_mode: "constant"      # Explicit: 9-param optimization (2 averaged + 7 physical)
      # per_angle_mode: "fourier"       # Explicit: Fourier reparameterization
      # per_angle_mode: "individual"    # Explicit: Independent per-angle params
      constant_scaling_threshold: 3    # Threshold for auto mode: N >= 3 → constant
      fourier_order: 2                  # Number of Fourier harmonics (if fourier mode)
      fourier_auto_threshold: 6         # Legacy: only used if auto selected fourier

      # Layer 2: Hierarchical Optimization
      hierarchical:
        enable: true
        max_outer_iterations: 5
        outer_tolerance: 1.0e-6
        physical_max_iterations: 100
        per_angle_max_iterations: 50

      # Layer 3: Adaptive Relative Regularization
      regularization:
        mode: "relative"            # "absolute", "relative", "auto"
        lambda: 1.0                 # Base lambda (100x higher than previous)
        target_cv: 0.10             # 10% variation target
        target_contribution: 0.10   # 10% of MSE contribution
        max_cv: 0.20                # Hard limit: 20% max variation
        auto_tune_lambda: true      # Compute lambda automatically

      # Layer 4: Gradient Collapse Detection
      gradient_monitoring:
        enable: true
        ratio_threshold: 0.01       # |∇_physical|/|∇_per_angle| threshold
        consecutive_triggers: 5     # Must trigger N times
        response: "hierarchical"    # "warn", "hierarchical", "reset", "abort"

      # Layer 5: Shear-Sensitivity Weighting
      shear_weighting:
        enable: true              # Enable angle-dependent loss weighting
        min_weight: 0.3           # Minimum weight for perpendicular angles
        alpha: 1.0                # Shear sensitivity exponent (1 = linear)
        update_frequency: 1       # Update weights every N outer iterations
        normalize: true           # Normalize weights so mean = 1.0

Configuration Options Summary¶

Layer 1 Configuration¶
Option	Default	Description
`per_angle_mode`	“auto”	DEFAULT: “auto” selects constant (N>=3) or individual (N<3). Options: “auto”, “constant”, “fourier”, “individual”
`constant_scaling_threshold`	3	Threshold for auto mode: N >= threshold selects constant mode
`fourier_order`	2	Number of Fourier harmonics (K) for fourier mode
`fourier_auto_threshold`	6	Legacy: only used if auto mode selected fourier

Layer 2 Configuration¶
Option	Default	Description
`hierarchical.enable`	true	Enable hierarchical optimization
`hierarchical.max_outer_iterations`	5	Maximum alternating iterations
`hierarchical.outer_tolerance`	1e-6	Convergence tolerance for physical params

Layer 3 Configuration¶
Option	Default	Description
`regularization.mode`	“relative”	Regularization type: “absolute”, “relative”, “auto”
`regularization.lambda`	1.0	Base regularization strength (100× previous default)
`regularization.target_cv`	0.10	Target coefficient of variation (10%)
`regularization.target_contribution`	0.10	Target contribution to loss (10%)
`regularization.max_cv`	0.20	Maximum allowed CV before violation

Layer 4 Configuration¶
Option	Default	Description
`gradient_monitoring.enable`	true	Enable gradient collapse detection
`gradient_monitoring.ratio_threshold`	0.01	Collapse detection threshold
`gradient_monitoring.consecutive_triggers`	5	Required consecutive triggers
`gradient_monitoring.response`	“hierarchical”	Response action on collapse

Layer 5 Configuration¶
Option	Default	Description
`shear_weighting.enable`	true	Enable shear-sensitivity weighting
`shear_weighting.min_weight`	0.3	Minimum weight for perpendicular angles
`shear_weighting.alpha`	1.0	Sensitivity exponent (1.0 = linear)
`shear_weighting.update_frequency`	1	Update weights every N outer iterations
`shear_weighting.normalize`	true	Normalize weights so mean = 1.0

Usage Tutorial¶

Basic Usage¶

For most users, the default configuration provides robust fitting:

from homodyne.optimization import fit_nlsq_jax

# Load your data
data = load_xpcs_data("experiment.hdf5")

# Fit with default anti-degeneracy settings
result = fit_nlsq_jax(
    t1=data.t1,
    t2=data.t2,
    c2=data.c2,
    phi_rad=data.phi_angles,
    mode="laminar_flow",
    # Anti-degeneracy defense is enabled by default for n_phi > 6
)

# gamma_dot_t0 should now be correctly estimated
# result.parameters is an ndarray; use parameter_names for labels
from homodyne.config.parameter_names import get_physical_param_names
names = get_physical_param_names("laminar_flow")
for name, val in zip(names, result.parameters):
    print(f"{name} = {val:.4e}")

Advanced Usage: Fine-Tuning¶

For challenging datasets, you may need to adjust settings:

from homodyne.config import ConfigManager

# Load configuration from YAML
config = ConfigManager("my_config.yaml")

# Or modify programmatically
config.optimization.nlsq.anti_degeneracy.hierarchical.max_outer_iterations = 10
config.optimization.nlsq.anti_degeneracy.regularization.target_cv = 0.05  # Tighter

# Run fit with custom config
result = fit_nlsq_jax(
    data,
    config=config,
    mode="laminar_flow",
)

Diagnosing Problems¶

If fitting still produces collapsed parameters:

Step 1: Check diagnostics

# Access optimization diagnostics
diagnostics = result.diagnostics

# Check if gradient collapse was detected
if diagnostics.get("gradient_monitor", {}).get("collapse_detected"):
    print("Gradient collapse was detected!")
    print(f"Response taken: {diagnostics['gradient_monitor']['response_actions']}")

# Check per-angle contrast from diagnostics
contrast_vals = diagnostics.get("per_angle_contrast", [])
contrast_cv = np.nanstd(contrast_vals) / np.nanmean(contrast_vals) if len(contrast_vals) else 0.0
if contrast_cv > 0.20:
    print(f"Warning: High contrast variation (CV={contrast_cv:.2%})")

Step 2: Increase regularization

optimization:
  nlsq:
    anti_degeneracy:
      regularization:
        lambda: 10.0          # Increase from 1.0
        target_cv: 0.05       # Tighter constraint

Step 3: Force Fourier mode

optimization:
  nlsq:
    anti_degeneracy:
      per_angle_mode: "fourier"  # Force Fourier, even for small n_phi
      fourier_order: 1           # Reduce degrees of freedom further

Step 4: Increase hierarchical iterations

optimization:
  nlsq:
    anti_degeneracy:
      hierarchical:
        max_outer_iterations: 20  # More alternations
        outer_tolerance: 1.0e-8   # Tighter convergence

When to Disable¶

The anti-degeneracy system can be disabled for:

Small n_phi (≤ 3): Structural degeneracy is less problematic
Static mode: No shear parameters to collapse
Quick exploratory fits: When speed matters more than accuracy

optimization:
  nlsq:
    anti_degeneracy:
      per_angle_mode: "independent"
      hierarchical:
        enable: false
      regularization:
        enable: false
      gradient_monitoring:
        enable: false

API Reference¶

Parameter Utilities¶

Parameter Utilities for NLSQ Optimization.

Provides utility functions for parameter handling, labeling, status classification, and per-angle initialization in NLSQ optimization.

Key Functions: - build_parameter_labels: Create parameter labels with per-angle support - classify_parameter_status: Identify parameters at bounds - sample_xdata: Subsample x-data for diagnostic computations - compute_consistent_per_angle_init: Initialize per-angle params consistently - compute_jacobian_stats: Compute Jacobian-based statistics

homodyne.optimization.nlsq.parameter_utils.build_parameter_labels(per_angle_scaling, n_phi, physical_param_names)[source]

Build parameter labels including per-angle scaling parameters.

Parameters:

per_angle_scaling (bool) – Whether per-angle contrast/offset are used
n_phi (int) – Number of phi angles
physical_param_names (list[str]) – Names of physical parameters

Returns:

Full list of parameter labels

Return type:

list[str]

homodyne.optimization.nlsq.parameter_utils.classify_parameter_status(values, lower, upper, atol=1e-9)[source]

Classify parameters as active or at bounds.

Parameters:

values (ndarray) – Current parameter values
lower (ndarray | None) – Lower bounds
upper (ndarray | None) – Upper bounds
atol (float) – Absolute tolerance for bound comparison

Returns:

Status for each parameter: ‘active’, ‘at_lower_bound’, or ‘at_upper_bound’

Return type:

list[str]

homodyne.optimization.nlsq.parameter_utils.sample_xdata(xdata, max_points)[source]

Subsample x-data for diagnostic computations.

Parameters:

xdata (ndarray) – Input data
max_points (int) – Maximum number of points to return

Returns:

Subsampled data

Return type:

ndarray

homodyne.optimization.nlsq.parameter_utils.compute_jacobian_stats(residual_fn, x_subset, params, scaling_factor)[source]

Compute Jacobian statistics for diagnostics.

Parameters:

residual_fn (Callable[..., Any]) – Residual function
x_subset (ndarray) – Subset of x data for computation
params (ndarray) – Current parameters
scaling_factor (float) – Scaling factor for statistics

Returns:

(J^T J matrix, column norms) or (None, None) on failure

Return type:

tuple[ndarray | None, ndarray | None]

homodyne.optimization.nlsq.parameter_utils.compute_consistent_per_angle_init(stratified_data, physical_params, physical_param_names, default_contrast=0.5, default_offset=1.0, logger=None)[source]

Compute per-angle contrast/offset consistent with initial physical parameters.

This function solves a critical initialization problem in laminar_flow mode: when physical shear parameters (gamma_dot_t0) are nonzero, the model predicts DIFFERENT g2 values at different angles. If per-angle contrast/offset are initialized uniformly, large initial residuals can cause the optimizer to incorrectly reduce gamma_dot_t0 to zero.

Instead, we compute per-angle contrast/offset by fitting:: g2_data[angle] ≈ offset[angle] + contrast[angle] × g1_model²[angle]

where g1_model is computed using the initial physical parameters.

Parameters:

stratified_data (Any) – Data containing per-angle g2, phi, t1, t2 arrays
physical_params (ndarray) – Initial physical parameters [D0, alpha, D_offset, (gamma_dot_t0, beta, gamma_dot_t_offset, phi0)]
physical_param_names (list[str]) – Names of physical parameters to determine analysis mode
default_contrast (float) – Default contrast value if fitting fails
default_offset (float) – Default offset value if fitting fails
logger (Any) – Logger for diagnostic messages

Return type:

tuple[ndarray, ndarray]

Returns:

contrast_per_angle (np.ndarray) – Per-angle contrast values consistent with physical params
offset_per_angle (np.ndarray) – Per-angle offset values consistent with physical params

homodyne.optimization.nlsq.parameter_utils.compute_quantile_per_angle_scaling(stratified_data, contrast_bounds=(0.0, 1.0), offset_bounds=(0.5, 1.5), lag_floor_quantile=0.80, lag_ceiling_quantile=0.20, value_quantile_low=0.10, value_quantile_high=0.90, logger=None)[source]

Estimate per-angle contrast/offset from quantiles of c2_experimental values.

This function uses physics-informed quantile analysis to estimate contrast and offset for each phi angle independently. Unlike least-squares fitting, this approach does not require a model and directly extracts scaling from the data.

Physics basis:

C2 = contrast × g1² + offset

At large time lags, g1² → 0, so C2 → offset (the “floor”)
At small time lags, g1² ≈ 1, so C2 ≈ contrast + offset (the “ceiling”)

Parameters:

stratified_data (Any) – Data containing per-angle g2_flat, phi_flat, t1_flat, t2_flat arrays.
contrast_bounds (tuple[float, float]) – Valid bounds for contrast parameter.
offset_bounds (tuple[float, float]) – Valid bounds for offset parameter.
lag_floor_quantile (float) – Quantile threshold for “large lag” region (default: 0.80 = top 20% of lags).
lag_ceiling_quantile (float) – Quantile threshold for “small lag” region (default: 0.20 = bottom 20% of lags).
value_quantile_low (float) – Quantile for robust floor estimation (default: 0.10).
value_quantile_high (float) – Quantile for robust ceiling estimation (default: 0.90).
logger (Any) – Logger for diagnostic messages.

Return type:

tuple[ndarray, ndarray]

Returns:

contrast_per_angle (np.ndarray) – Per-angle contrast values from quantile estimation.
offset_per_angle (np.ndarray) – Per-angle offset values from quantile estimation.

Notes

The estimation is robust to outliers by using quantiles instead of min/max. The lag-based segmentation ensures we sample from appropriate regions of the correlation decay curve.

This function is designed for the “constant” mode in anti-degeneracy defense, where per-angle contrast/offset are estimated once and treated as fixed parameters during optimization.

Fourier Reparameterization¶

Fourier Reparameterization for Anti-Degeneracy Defense.

This module replaces n_phi independent per-angle contrast/offset values with truncated Fourier series, dramatically reducing structural degeneracy.

Part of Anti-Degeneracy Defense System v2.9.0. See: docs/specs/anti-degeneracy-defense-v2.9.0.md

Mathematical Formulation¶

contrast(φ) = c₀ + Σₖ[cₖ×cos(kφ) + sₖ×sin(kφ)] for k=1..order offset(φ) = o₀ + Σₖ[oₖ×cos(kφ) + tₖ×sin(kφ)] for k=1..order

For order=2: - Contrast: 5 coefficients [c₀, c₁, s₁, c₂, s₂] - Offset: 5 coefficients [o₀, o₁, t₁, o₂, t₂] - Total: 10 Fourier coefficients vs 2×n_phi independent params

Parameter Count Comparison:

n_phi | Independent | Fourier (order=2) | Reduction
------|-------------|-------------------|----------
 |     4       |        4          |    0%
 |     6       |        6          |    0%
 |    20       |       10          |   50%
 |    46       |       10          |   78%
 |   200       |       10          |   95%

Note: For n_phi <= 2*(order+1), independent mode is used.

class homodyne.optimization.nlsq.fourier_reparam.FourierReparamConfig[source]

Bases: object

Configuration for Fourier reparameterization.

mode

Per-angle parameter mode: - “independent”: Use n_phi independent contrast/offset values - “fourier”: Use truncated Fourier series - “auto”: Use Fourier when n_phi > auto_threshold

Type:: str

fourier_order

Number of Fourier harmonics. Default 2. order=2 gives 5 coefficients per parameter (c0, c1, s1, c2, s2).

Type:: int

auto_threshold

Use Fourier when n_phi > this threshold in auto mode. Default 6.

Type:: int

c0_bounds

Bounds for mean contrast coefficient. Default (0.1, 0.8).

Type:: tuple

ck_bounds

Bounds for harmonic contrast amplitudes. Default (-0.2, 0.2).

Type:: tuple

o0_bounds

Bounds for mean offset coefficient. Default (0.5, 1.5).

Type:: tuple

ok_bounds

Bounds for harmonic offset amplitudes. Default (-0.3, 0.3).

Type:: tuple

mode: Literal['independent', 'fourier', 'auto'] = 'auto'

fourier_order: int = 2

auto_threshold: int = 6

c0_bounds: tuple[float, float] = (0.1, 0.8)

ck_bounds: tuple[float, float] = (-0.2, 0.2)

o0_bounds: tuple[float, float] = (0.5, 1.5)

ok_bounds: tuple[float, float] = (-0.3, 0.3)

classmethod from_dict(config_dict)[source]

Create config from dictionary.

Return type:: FourierReparamConfig

__init__(mode='auto', fourier_order=2, auto_threshold=6, c0_bounds=(0.1, 0.8), ck_bounds=(-0.2, 0.2), o0_bounds=(0.5, 1.5), ok_bounds=(-0.3, 0.3))

class homodyne.optimization.nlsq.fourier_reparam.FourierReparameterizer[source]

Bases: object

Handles conversion between Fourier coefficients and per-angle values.

This class provides the core functionality for Fourier reparameterization: 1. Convert per-angle values to Fourier coefficients (initialization) 2. Convert Fourier coefficients to per-angle values (model evaluation) 3. Compute Jacobian for covariance transformation

The Fourier basis ensures smooth variation of contrast/offset with angle, preventing the optimizer from using per-angle parameters to absorb angle-dependent physical signals (like the shear term cos(φ₀-φ)).

Parameters:

phi_angles (ndarray) – Unique phi angles in radians, shape (n_phi,).
config (FourierReparamConfig) – Fourier configuration.

n_phi

Number of unique phi angles.

Type:: int

n_coeffs

Total number of Fourier coefficients (contrast + offset).

Type:: int

n_coeffs_per_param

Coefficients per parameter type (contrast or offset).

Type:: int

use_fourier

Whether Fourier mode is active.

Type:: bool

Examples

>>> phi_angles = np.linspace(-np.pi, np.pi, 23)
>>> config = FourierReparamConfig(mode="fourier", fourier_order=2)
>>> fourier = FourierReparameterizer(phi_angles, config)
>>> # Convert initial per-angle values to Fourier
>>> contrast = np.full(23, 0.3)
>>> offset = np.full(23, 1.0)
>>> fourier_coeffs = fourier.per_angle_to_fourier(contrast, offset)
>>> # Convert back during model evaluation
>>> contrast_out, offset_out = fourier.fourier_to_per_angle(fourier_coeffs)

__init__(phi_angles, config)[source]

Initialize Fourier reparameterizer.

Parameters:

phi_angles (ndarray) – Unique phi angles in radians, shape (n_phi,).
config (FourierReparamConfig) – Fourier configuration.

get_basis_matrix()[source]

Get the Fourier basis matrix for covariance transformation.

Returns:: Basis matrix of shape (n_phi, n_coeffs_per_param) if in Fourier mode, None if in independent mode. The basis matrix B satisfies: per_angle_values = B @ fourier_coeffs
Return type:: ndarray | None

Notes

Used for transforming covariance from Fourier space to per-angle space: pcov_per_angle = B @ pcov_fourier @ B.T

property order: int

Get the Fourier order (number of harmonics).

Returns:: Fourier order from config.
Return type:: int

fourier_to_per_angle(fourier_coeffs)[source]

Convert Fourier coefficients to per-angle contrast/offset.

Parameters:

fourier_coeffs (ndarray) – Shape (n_coeffs,) = [c0,c1,s1,c2,s2,…,o0,o1,t1,o2,t2,…].

Return type:

tuple[ndarray, ndarray]

Returns:

contrast (np.ndarray) – Per-angle contrast values, shape (n_phi,).
offset (np.ndarray) – Per-angle offset values, shape (n_phi,).

Raises:

ValueError – If fourier_coeffs has wrong shape.

per_angle_to_fourier(contrast, offset)[source]

Convert per-angle values to Fourier coefficients.

Uses least squares fitting when n_phi > n_coeffs_per_param.

Parameters:

contrast (ndarray) – Per-angle contrast values, shape (n_phi,).
offset (ndarray) – Per-angle offset values, shape (n_phi,).

Returns:

Fourier coefficients, shape (n_coeffs,).

Return type:

ndarray

Raises:

ValueError – If contrast or offset has wrong shape.

get_jacobian_transform()[source]

Get Jacobian of transformation: d(per_angle)/d(fourier).

Used for covariance transformation back to per-angle space:: Cov_per_angle = J @ Cov_fourier @ J.T

Returns:: Jacobian matrix of shape (2*n_phi, n_coeffs).
Return type:: ndarray

get_bounds()[source]

Get bounds for Fourier coefficients.

Return type:

tuple[ndarray, ndarray]

Returns:

lower (np.ndarray) – Lower bounds, shape (n_coeffs,).
upper (np.ndarray) – Upper bounds, shape (n_coeffs,).

get_initial_coefficients(contrast_init, offset_init)[source]

Get initial Fourier coefficients from initial values.

Parameters:

contrast_init (float | ndarray) – Initial contrast (scalar for uniform, array for per-angle).
offset_init (float | ndarray) – Initial offset (scalar for uniform, array for per-angle).

Returns:

Initial Fourier coefficients.

Return type:

ndarray

get_coefficient_labels()[source]

Get parameter labels for Fourier coefficients.

Returns:: Parameter labels.
Return type:: list[str]

to_fourier(per_angle_values)[source]

Convert a single per-angle array to Fourier coefficients.

Convenience method for transforming one group (contrast or offset) at a time, rather than both together.

Parameters:: per_angle_values (ndarray) – Per-angle values, shape (n_phi,).
Returns:: Fourier coefficients, shape (n_coeffs_per_param,).
Return type:: ndarray
Raises:: ValueError – If per_angle_values has wrong shape.

from_fourier(fourier_coeffs)[source]

Convert Fourier coefficients to per-angle values for a single group.

Convenience method for transforming one group (contrast or offset) at a time, rather than both together.

Parameters:: fourier_coeffs (ndarray) – Fourier coefficients, shape (n_coeffs_per_param,).
Returns:: Per-angle values, shape (n_phi,).
Return type:: ndarray
Raises:: ValueError – If fourier_coeffs has wrong shape.

get_diagnostics()[source]

Get Fourier reparameterization diagnostics.

Returns:: Diagnostic information.
Return type:: dict

homodyne.optimization.nlsq.fourier_reparam.create_fourier_model_wrapper(model_fn, fourier, n_physical)[source]

Create a model function wrapper that handles Fourier conversion.

The wrapper converts Fourier coefficients to per-angle values before calling the underlying model function.

Parameters:

model_fn (Callable[[ndarray, ndarray], ndarray]) – Original model function that expects per-angle parameters: f(params, x) where params = [contrast_per_angle, offset_per_angle, physical]
fourier (FourierReparameterizer) – Fourier reparameterizer instance.
n_physical (int) – Number of physical parameters.

Returns:

Wrapped model function that accepts Fourier parameters: f(params, x) where params = [fourier_coeffs, physical]

Return type:

Callable[[ndarray, ndarray], ndarray]

Hierarchical Optimization¶

Hierarchical Two-Stage Optimization for Anti-Degeneracy Defense.

This module implements alternating optimization between physical and per-angle parameters, breaking the gradient cancellation cycle that causes structural degeneracy in streaming optimization.

Part of Anti-Degeneracy Defense System v2.9.0. See: docs/specs/anti-degeneracy-defense-v2.9.0.md

Algorithm:

Initialize: params = [per_angle_params, physical_params]

for outer_iter in range(max_outer_iterations):

    # Stage 1: Fit PHYSICAL params only
    freeze(per_angle_params)
    result1 = L-BFGS(
        loss_fn(physical_params | frozen_per_angle),
        physical_params
    )
    physical_params = result1.x

    # Stage 2: Fit PER-ANGLE params only
    freeze(physical_params)
    result2 = L-BFGS(
        loss_fn(per_angle_params | frozen_physical),
        per_angle_params
    )
    per_angle_params = result2.x

    # Check convergence
    if converged(physical_params, previous_physical_params):
        break

return [per_angle_params, physical_params]

Why It Works¶

In Stage 1, there are NO per-angle DoF to compete with physical params
gamma_dot_t0 gradient CANNOT cancel (no per-angle params to absorb signal)
Physical params converge to true values
Stage 2 only cleans up residuals with physical interpretation fixed

class homodyne.optimization.nlsq.hierarchical.HierarchicalConfig[source]

Bases: object

Configuration for hierarchical optimization.

enable

Whether to enable hierarchical optimization. Default True.

Type:: bool

max_outer_iterations

Maximum outer iterations. Default 5.

Type:: int

outer_tolerance

Convergence tolerance for physical parameters. Default 1e-6.

Type:: float

physical_max_iterations

Max iterations for Stage 1 (physical params). Default 100.

Type:: int

physical_ftol

Function tolerance for Stage 1. Default 1e-8.

Type:: float

per_angle_max_iterations

Max iterations for Stage 2 (per-angle params). Default 50.

Type:: int

per_angle_ftol

Function tolerance for Stage 2. Default 1e-6.

Type:: float

log_stage_transitions

Whether to log stage transitions. Default True.

Type:: bool

save_intermediate_results

Whether to save intermediate results. Default False.

Type:: bool

enable: bool = True

max_outer_iterations: int = 5

outer_tolerance: float = 1e-06

physical_max_iterations: int = 100

physical_ftol: float = 1e-08

per_angle_max_iterations: int = 50

per_angle_ftol: float = 1e-06

log_stage_transitions: bool = True

save_intermediate_results: bool = False

classmethod from_dict(config_dict)[source]

Create config from dictionary with safe type conversion.

Return type:: HierarchicalConfig

__init__(enable=True, max_outer_iterations=5, outer_tolerance=1e-06, physical_max_iterations=100, physical_ftol=1e-08, per_angle_max_iterations=50, per_angle_ftol=1e-06, log_stage_transitions=True, save_intermediate_results=False)

class homodyne.optimization.nlsq.hierarchical.HierarchicalResult[source]

Bases: object

Result from hierarchical optimization.

x

Optimized parameters.

Type:: np.ndarray

fun

Final loss value.

Type:: float

success

Whether optimization succeeded.

Type:: bool

n_outer_iterations

Number of outer iterations performed.

Type:: int

history

History of each outer iteration.

Type:: list

total_time

Total optimization time in seconds.

Type:: float

message

Status message.

Type:: str

x: ndarray

fun: float

success: bool

n_outer_iterations: int

history: list[dict]

total_time: float = 0.0

message: str = ''

__init__(x, fun, success, n_outer_iterations, history=<factory>, total_time=0.0, message='')

class homodyne.optimization.nlsq.hierarchical.HierarchicalOptimizer[source]

Bases: object

Two-stage hierarchical optimizer for decoupled fitting.

This optimizer breaks the gradient cancellation problem by alternating between physical and per-angle parameter optimization:

Stage 1: Physical parameters only

Per-angle parameters are frozen
gamma_dot_t0 gradient cannot be cancelled by per-angle absorption
Physical params converge to true values

Stage 2: Per-angle parameters only

Physical parameters are frozen
Per-angle params absorb only experimental noise
Cannot change the physical interpretation

Parameters:

config (HierarchicalConfig) – Hierarchical optimization configuration.
n_phi (int) – Number of unique phi angles.
n_physical (int) – Number of physical parameters.
fourier_reparameterizer (FourierReparameterizer | None) – Fourier reparameterizer if using Fourier mode.

Examples

>>> config = HierarchicalConfig(max_outer_iterations=5)
>>> optimizer = HierarchicalOptimizer(config, n_phi=23, n_physical=7)
>>> result = optimizer.fit(loss_fn, grad_fn, p0, bounds)

__init__(config, n_phi, n_physical, fourier_reparameterizer=None)[source]

Initialize hierarchical optimizer.

Parameters:

config (HierarchicalConfig) – Configuration.
n_phi (int) – Number of unique phi angles.
n_physical (int) – Number of physical parameters.
fourier_reparameterizer (FourierReparameterizer | None) – Fourier reparameterizer for Fourier mode.

per_angle_indices: ndarray

physical_indices: ndarray

fit(loss_fn, grad_fn, p0, bounds, outer_iteration_callback=None)[source]

Run hierarchical optimization.

Parameters:

loss_fn (Callable[[ndarray], float]) – Loss function f(params) -> scalar.
grad_fn (Callable[[ndarray], ndarray] | None) – Gradient function g(params) -> gradient array. If None, uses finite differences.
p0 (ndarray) – Initial parameters.
bounds (tuple[ndarray, ndarray]) – (lower_bounds, upper_bounds).
outer_iteration_callback (Callable[[ndarray, int], None] | None) – Optional callback called at the start of each outer iteration. Signature: callback(current_params, outer_iter). Used for updating shear-sensitivity weights based on current phi0 estimate.

Returns:

Optimization result with diagnostics.

Return type:

HierarchicalResult

get_diagnostics()[source]

Get optimizer diagnostics.

Returns:: Diagnostic information.
Return type:: dict

Adaptive Regularization¶

Adaptive Relative Regularization for Anti-Degeneracy Defense.

This module implements CV-based (Coefficient of Variation) regularization that scales properly with data, replacing the ineffective absolute variance regularization.

Part of Anti-Degeneracy Defense System v2.9.0. See: docs/specs/anti-degeneracy-defense-v2.9.0.md

Mathematical Formulation:

Current (ineffective):
    L_reg = lambda * Var(params) * n_points

Proposed (CV-based):
    CV = std(params) / abs(mean(params))
    L_reg = lambda * CV^2 * MSE * n_points

Auto-tuned lambda:
    lambda = target_contribution / target_cv^2

    Example: Allow 10% variation (CV=0.1), contribute 10% to loss
    lambda = 0.1 / 0.01 = 10

class homodyne.optimization.nlsq.adaptive_regularization.AdaptiveRegularizationConfig[source]

Bases: object

Configuration for adaptive relative regularization.

enable

Whether to enable regularization. Default True.

Type:: bool

mode

Regularization mode: “absolute”, “relative”, or “auto”. - “absolute”: Original variance-based (L_reg = λ × Var × n) - “relative”: CV-based (L_reg = λ × CV² × MSE × n) - “auto”: Use relative for n_phi > 5, absolute otherwise

Type:: str

lambda_base

Base regularization strength. Default 1.0 (100× stronger than v2.8).

Type:: float

target_cv

Target coefficient of variation. Default 0.10 (10% variation allowed).

Type:: float

target_contribution

Target fraction of MSE to contribute. Default 0.10 (10% of loss).

Type:: float

auto_tune_lambda

Whether to auto-compute λ from target_cv and target_contribution.

Type:: bool

max_cv

Maximum allowed CV before hard constraint warning. Default 0.20.

Type:: float

group_indices

Parameter group indices [(start, end), …]. Auto-computed if None.

Type:: list of tuple, optional

enable: bool = True

mode: Literal['absolute', 'relative', 'auto'] = 'relative'

lambda_base: float = 1.0

target_cv: float = 0.1

target_contribution: float = 0.1

auto_tune_lambda: bool = True

max_cv: float = 0.2

group_indices: list[tuple[int, int]] | None = None

classmethod from_dict(config_dict)[source]

Create config from dictionary with safe type conversion.

Return type:: AdaptiveRegularizationConfig

__init__(enable=True, mode='relative', lambda_base=1.0, target_cv=0.1, target_contribution=0.1, auto_tune_lambda=True, max_cv=0.2, group_indices=None)

class homodyne.optimization.nlsq.adaptive_regularization.AdaptiveRegularizer[source]

Bases: object

CV-based adaptive regularization for per-angle parameters.

This regularizer addresses the fundamental problem where absolute variance regularization (λ=0.01) contributed only ~0.05% to total loss, providing no effective constraint on per-angle parameter variation.

The CV-based approach ensures regularization scales properly: - CV is dimensionless (ratio of std to mean) - Auto-tuned λ makes regularization ~10% of MSE - Prevents per-angle parameters from absorbing physical signals

Parameters:

config (AdaptiveRegularizationConfig) – Regularization configuration.
n_phi (int) – Number of unique phi angles.

lambda_value

Effective regularization strength (auto-tuned or from config).

Type:: float

group_indices

Parameter groups to regularize.

Type:: list of tuple

Examples

>>> config = AdaptiveRegularizationConfig(target_cv=0.10, target_contribution=0.10)
>>> regularizer = AdaptiveRegularizer(config, n_phi=23)
>>> reg_term = regularizer.compute_regularization(
...     params, mse=0.04, n_points=23_000_000
... )

__init__(config, n_phi, n_params=None)[source]

Initialize adaptive regularizer.

Parameters:

config (AdaptiveRegularizationConfig) – Regularization configuration.
n_phi (int) – Number of unique phi angles.
n_params (int | None) – Actual parameter vector length. When provided and less than 2 * n_phi + n_physical, auto_averaged mode is assumed (2 scaling params instead of 2 * n_phi).

compute_regularization(params, mse, n_points)[source]

Compute regularization term to add to loss.

Parameters:

params (ndarray) – Full parameter vector.
mse (float) – Current mean squared error.
n_points (int) – Number of data points.

Returns:

Regularization term to add to loss (SSE scale).

Return type:

float

compute_regularization_jax(params, mse, n_points)[source]

Compute regularization term using JAX for autodiff compatibility.

This method uses JAX operations (jnp) instead of NumPy, making it compatible with JAX’s JIT compilation and autodiff (jax.grad).

Use this method when the regularization needs to be part of a differentiable loss function.

Parameters:

params (Array) – Full parameter vector (JAX array, possibly traced).
mse (Array) – Current mean squared error (JAX scalar, possibly traced).
n_points (int) – Number of data points.

Returns:

Regularization term to add to loss (SSE scale, JAX scalar).

Return type:

Array

compute_regularization_gradient(params, mse, n_points)[source]

Compute gradient of regularization term.

Parameters:

params (ndarray) – Full parameter vector.
mse (float) – Current mean squared error.
n_points (int) – Number of data points.

Returns:

Gradient w.r.t. all parameters (zeros for non-regularized params).

Return type:

ndarray

check_constraint_violation(params)[source]

Check if CV exceeds max_cv threshold.

Parameters:: params (ndarray) – Full parameter vector.
Returns:: Dictionary of violations, empty if none.
Return type:: dict[str, dict]

get_diagnostics()[source]

Get regularization diagnostics for logging.

Returns:: Diagnostic information including CV values and contribution.
Return type:: dict

log_summary(params, mse, n_points)[source]

Log regularization summary.

Parameters:

params (ndarray) – Full parameter vector.
mse (float) – Current mean squared error.
n_points (int) – Number of data points.

Return type:

None

Gradient Monitoring¶

Gradient Collapse Monitor for Anti-Degeneracy Defense.

This module provides runtime detection of gradient collapse (physical params losing gradient signal) with automatic response actions.

Part of Anti-Degeneracy Defense System v2.9.0. See: docs/specs/anti-degeneracy-defense-v2.9.0.md

Detection Mechanism:

Monitor the ratio:
    ratio = norm(grad_physical) / norm(grad_per_angle)

If ratio < threshold for N consecutive iterations:
    - Gradient collapse detected
    - Physical params are losing signal to per-angle params

Response Actions¶

“warn”: Log warning only
“hierarchical”: Switch to hierarchical optimization mode
“reset”: Reset per-angle params to mean values
“abort”: Abort optimization and return best params so far

class homodyne.optimization.nlsq.gradient_monitor.GradientMonitorConfig[source]

Bases: object

Configuration for gradient collapse detection.

enable

Whether to enable gradient monitoring. Default True.

Type:: bool

ratio_threshold

Ratio of norm(grad_physical) / norm(grad_per_angle) below this triggers detection. Default 0.01 (physical gradient is 1% of per-angle gradient).

Type:: float

consecutive_triggers

Must trigger N consecutive times to confirm collapse. Default 5.

Type:: int

response_mode

Response action on collapse detection: - “warn”: Log warning only - “hierarchical”: Switch to hierarchical optimization - “reset”: Reset per-angle params to mean - “abort”: Abort and return best params

Type:: str

reset_per_angle_to_mean

When resetting, reset per-angle to mean values. Default True.

Type:: bool

lambda_multiplier_on_collapse

Multiply regularization λ by this on collapse. Default 10.0.

Type:: float

check_interval

Check every N iterations. Default 1 (every iteration).

Type:: int

enable: bool = True

ratio_threshold: float = 0.01

consecutive_triggers: int = 5

response_mode: Literal['warn', 'hierarchical', 'reset', 'abort'] = 'hierarchical'

reset_per_angle_to_mean: bool = True

lambda_multiplier_on_collapse: float = 10.0

check_interval: int = 1

watch_parameters: list[int] | None = None

watch_threshold: float = 1e-08

watch_consecutive_triggers: int = 3

watch_min_iteration: int = 5

classmethod from_dict(config_dict)[source]

Create config from dictionary with safe type conversion.

Return type:: GradientMonitorConfig

__init__(enable=True, ratio_threshold=0.01, consecutive_triggers=5, response_mode='hierarchical', reset_per_angle_to_mean=True, lambda_multiplier_on_collapse=10.0, check_interval=1, watch_parameters=None, watch_threshold=1e-08, watch_consecutive_triggers=3, watch_min_iteration=5)

class homodyne.optimization.nlsq.gradient_monitor.CollapseEvent[source]

Bases: object

Record of a gradient collapse event.

iteration

Iteration when collapse was detected.

Type:: int

ratio

Gradient ratio at detection.

Type:: float

physical_grad_norm

Physical parameter gradient norm.

Type:: float

per_angle_grad_norm

Per-angle parameter gradient norm.

Type:: float

response_mode

Response action taken.

Type:: str

iteration: int

ratio: float

physical_grad_norm: float

per_angle_grad_norm: float

response_mode: str

__init__(iteration, ratio, physical_grad_norm, per_angle_grad_norm, response_mode)

class homodyne.optimization.nlsq.gradient_monitor.GradientCollapseMonitor[source]

Bases: object

Monitor for detecting and responding to gradient collapse.

This monitor tracks the ratio of physical to per-angle gradient norms during optimization. When the ratio drops below a threshold for consecutive iterations, it indicates that physical parameters are losing gradient signal (being absorbed by per-angle parameters).

Parameters:

config (GradientMonitorConfig) – Monitor configuration.
physical_indices (Sequence[int] | ndarray) – Indices of physical parameters in the full parameter vector.
per_angle_indices (Sequence[int] | ndarray) – Indices of per-angle parameters in the full parameter vector.

collapse_detected

Whether gradient collapse has been detected.

Type:: bool

consecutive_count

Current count of consecutive low-ratio iterations.

Type:: int

Notes

History is capped at MAX_HISTORY_SIZE to prevent memory leaks during long-running optimizations. Older entries are discarded when the limit is reached.

Examples

>>> config = GradientMonitorConfig(ratio_threshold=0.01, consecutive_triggers=5)
>>> monitor = GradientCollapseMonitor(config, physical_indices=[6,7,8,9,10,11,12],
...                                    per_angle_indices=list(range(6)))
>>> for iter in range(100):
...     gradients = compute_gradients(params)
...     status = monitor.check(gradients, iter)
...     if status == "COLLAPSE_DETECTED":
...         response = monitor.get_response()
...         # Take action based on response

MAX_HISTORY_SIZE: int = 1000

__init__(config, physical_indices, per_angle_indices)[source]

Initialize gradient collapse monitor.

Parameters:

config (GradientMonitorConfig) – Monitor configuration.
physical_indices (Sequence[int] | ndarray) – Indices of physical parameters. Converted to numpy array internally to support both NumPy and JAX array indexing.
per_angle_indices (Sequence[int] | ndarray) – Indices of per-angle parameters (or Fourier coefficients when Fourier reparameterization is active). Converted to numpy array internally.

Notes

When Fourier reparameterization is active, per_angle_indices should correspond to Fourier coefficient indices (typically 10 for order=2), not independent per-angle indices (2 * n_phi).

physical_indices: ndarray

per_angle_indices: ndarray

history: deque[dict]

consecutive_count: int

collapse_detected: bool

collapse_events: list[CollapseEvent]

best_params: ndarray | None

best_loss: float

check(gradients, iteration, params=None, loss=None)[source]

Check for gradient collapse.

Parameters:

gradients (ndarray) – Full gradient vector.
iteration (int) – Current iteration number.
params (ndarray | None) – Current parameters (for response actions and tracking).
loss (float | None) – Current loss value (for tracking best params).

Returns:

Status: “OK”, “WARNING”, “COLLAPSE_DETECTED”

Return type:

str

get_response()[source]

Get response action after collapse detection.

Returns:: Response action dictionary, or None if no collapse.
Return type:: dict | None

compute_reset_params(params, n_phi)[source]

Compute parameters with per-angle values reset to mean.

Parameters:

params (ndarray) – Current parameter vector.
n_phi (int) – Number of phi angles.

Returns:

Parameters with per-angle values reset.

Return type:

ndarray

reset()[source]

Reset monitor state for new optimization run.

Return type:: None

get_diagnostics()[source]

Get monitoring diagnostics for logging.

Returns:: Diagnostic information.
Return type:: dict

log_summary()[source]

Log monitoring summary.

Return type:: None

homodyne.optimization.nlsq.gradient_monitor.create_gradient_function_with_monitoring(grad_fn, monitor)[source]

Wrap gradient function to include monitoring.

Parameters:

grad_fn (Callable[[ndarray], ndarray]) – Original gradient function.
monitor (GradientCollapseMonitor) – Monitor instance.

Returns:

Wrapped gradient function that records to monitor.

Return type:

Callable[[ndarray], ndarray]

Shear-Sensitivity Weighting¶

Shear-Sensitivity Weighting for Anti-Degeneracy Defense.

This module implements angle-dependent loss weighting to prevent gradient cancellation in the shear term during optimization.

Part of Anti-Degeneracy Defense System v2.9.1.

The Problem¶

The shear term gradient is:: d(g1_shear)/d(gamma_dot_t0) ~ cos(phi0 - phi)

When summed uniformly over all angles: - Angles near phi0: cos(phi0 - phi) ~ +1 (positive contribution) - Angles near phi0 +/- 90deg: cos ~ 0 (negligible) - Angles near phi0 +/- 180deg: cos ~ -1 (negative contribution)

With uniformly distributed angles, positive and negative contributions CANCEL, leading to near-zero net gradient for gamma_dot_t0. This causes the shear parameter to collapse to its lower bound.

The Solution¶

Use angle-dependent loss weighting:

L = sum_phi w(phi) * sum_tau (g2_model - g2_exp)^2

where w(phi) emphasizes shear-sensitive angles:

w(phi) = w_min + (1 - w_min) * abs(cos(phi0_current - phi))^alpha

This converts gradient cancellation into a weighted sum where shear-sensitive angles (parallel/antiparallel to flow) contribute more than perpendicular angles. All angles still contribute to prevent information loss.

Configuration¶

shear_weighting:: enable: true # Enable shear-sensitivity weighting min_weight: 0.3 # Minimum weight (0-1) alpha: 1.0 # Shear sensitivity exponent (1 = linear) update_frequency: 1 # Update weights every N outer iterations initial_phi0: null # Initial phi0 guess (null = use config)

class homodyne.optimization.nlsq.shear_weighting.ShearWeightingConfig[source]

Bases: object

Configuration for shear-sensitivity weighting.

enable

Enable shear-sensitivity weighting. Default True.

Type:: bool

min_weight

Minimum weight for perpendicular angles. Range [0, 1]. Default 0.3.

Type:: float

alpha

Shear sensitivity exponent. Higher = more aggressive weighting. Default 1.0 (linear).

Type:: float

update_frequency

Update weights every N outer iterations. Default 1.

Type:: int

initial_phi0

Initial phi0 guess in degrees. None = use config or 0.0.

Type:: float or None

normalize

Normalize weights so sum = n_phi. Default True.

Type:: bool

enable: bool = True

min_weight: float = 0.3

alpha: float = 1.0

update_frequency: int = 1

initial_phi0: float | None = None

normalize: bool = True

classmethod from_config(config)[source]

Create from configuration dictionary.

Parameters:: config (Mapping) – Configuration dictionary.
Returns:: Configuration object.
Return type:: ShearWeightingConfig

__init__(enable=True, min_weight=0.3, alpha=1.0, update_frequency=1, initial_phi0=None, normalize=True)

class homodyne.optimization.nlsq.shear_weighting.ShearSensitivityWeighting[source]

Bases: object

Shear-sensitivity weighted loss for anti-degeneracy defense.

This class manages angle-dependent weights that emphasize shear-sensitive angles during optimization, preventing gradient cancellation.

Parameters:

phi_angles (ndarray) – Array of phi angles in degrees.
n_physical (int) – Number of physical parameters.
phi0_index (int) – Index of phi0 in physical parameters (typically 6 for laminar_flow).
config (ShearWeightingConfig | None) – Weighting configuration.

Examples

>>> phi_angles = np.array([-30, 0, 30, 60, 90, 120])
>>> weighter = ShearSensitivityWeighting(phi_angles, n_physical=7, phi0_index=6)
>>> weights = weighter.get_weights(phi0_current=-5.0)
>>> # Angles near -5 deg and 175 deg get higher weight

__init__(phi_angles, n_physical, phi0_index, config=None)[source]

update_phi0(params, iteration=0)[source]

Update phi0 estimate from current parameters.

Parameters:

params (ndarray) – Current parameter vector. Physical parameters should be at the end.
iteration (int) – Current iteration number.

Return type:

None

get_weights(phi0_current=None)[source]

Get current angle weights.

Parameters:: phi0_current (float | None) – Override phi0 for weight computation. If None, uses stored value.
Returns:: Weight array of shape (n_phi,).
Return type:: ndarray

get_weights_jax()[source]

Get current angle weights as JAX array.

Returns:: Weight array of shape (n_phi,).
Return type:: Array

apply_weights_to_loss(residuals, phi_indices)[source]

Apply angle weights to residuals for loss computation.

Computes weighted mean squared error:: L = sum_i w[phi_idx[i]] * residuals[i]^2 / sum_i w[phi_idx[i]]

Parameters:

residuals (Array) – Residuals array of shape (n_data,).
phi_indices (Array) – Phi index for each data point, shape (n_data,).

Returns:

Weighted loss (scalar).

Return type:

Array

compute_weighted_mse(residuals, phi_indices)[source]

Compute weighted MSE (for gradient computation).

Parameters:

residuals (Array) – Residuals array of shape (n_data,).
phi_indices (Array) – Phi index for each data point, shape (n_data,).

Returns:

Weighted MSE (scalar).

Return type:

Array

get_diagnostics()[source]

Get weighting diagnostics.

Returns:: Diagnostic information.
Return type:: dict

property phi0_current: float: Current phi0 estimate in degrees.

homodyne.optimization.nlsq.shear_weighting.create_shear_weighting(phi_angles, n_physical, config=None, physical_param_names=None)[source]

Factory function to create shear weighting if enabled.

Parameters:

phi_angles (ndarray) – Phi angles in degrees.
n_physical (int) – Number of physical parameters.
config (Mapping | None) – Configuration dictionary.

Returns:

Weighting object if enabled, None otherwise.

Return type:

ShearSensitivityWeighting | None

Parameter Index Mapper¶

Centralized index mapping for anti-degeneracy layers.

This module provides the ParameterIndexMapper class which ensures consistent index ranges regardless of whether Fourier reparameterization is active. This is the single source of truth for parameter group boundaries.

Created: 2025-12-31 Feature: 001-fix-nlsq-anti-degeneracy

class homodyne.optimization.nlsq.parameter_index_mapper.ParameterIndexMapper[source]

Bases: object

Centralized index mapping for anti-degeneracy layers.

Provides consistent index ranges regardless of whether Fourier reparameterization or constant scaling is active. This class is the single source of truth for parameter group boundaries.

Parameters:

n_phi (int) – Number of unique phi angles.
n_physical (int) – Number of physical parameters (typically 7 for laminar_flow mode).
fourier (FourierReparameterizer | None) – Reference to Fourier reparameterizer if Layer 1 is active.
use_constant (bool) – Whether constant scaling mode is active (single contrast/offset shared across all angles).

n_per_angle_total

Total number of per-angle parameters (Fourier coefficients, raw, or 2).

Type:: int

n_per_group

Number of parameters per group (contrast or offset).

Type:: int

use_fourier

Whether Fourier reparameterization is active.

Type:: bool

use_constant

Whether constant scaling mode is active.

Type:: bool

total_params

Total number of parameters.

Type:: int

mode_name

Human-readable name of current mode (“constant”, “fourier”, or “individual”).

Type:: str

Examples

>>> # Constant mode (23 phi angles)
>>> mapper = ParameterIndexMapper(n_phi=23, n_physical=7, use_constant=True)
>>> mapper.get_group_indices()
[(0, 1), (1, 2)]
>>> mapper.n_per_angle_total
2
>>> mapper.mode_name
'constant'

>>> # Non-Fourier mode (23 phi angles)
>>> mapper = ParameterIndexMapper(n_phi=23, n_physical=7, fourier=None)
>>> mapper.get_group_indices()
[(0, 23), (23, 46)]
>>> mapper.n_per_angle_total
46

>>> # Fourier mode (23 phi angles, order=2)
>>> mapper = ParameterIndexMapper(n_phi=23, n_physical=7, fourier=fourier_obj)
>>> mapper.get_group_indices()
[(0, 5), (5, 10)]
>>> mapper.n_per_angle_total
10

n_phi: int

n_physical: int

fourier: FourierReparameterizer | None = None

use_constant: bool = False

__post_init__()[source]

Validate inputs and cache computed values.

Return type:: None

property use_fourier: bool: Check if Fourier reparameterization is active.

property n_per_group: int

Get number of parameters per group (contrast or offset).

Returns:: 1 for constant mode, n_coeffs for Fourier, n_phi for individual.
Return type:: int

property mode_name: str

Get human-readable name of current mode.

Returns:: “constant”, “fourier”, or “individual”
Return type:: str

property n_per_angle_total: int: Get total number of per-angle parameters (scaling params).

property total_params: int: Get total number of parameters.

get_group_indices()[source]

Get (start, end) tuples for contrast and offset parameter groups.

Returns:: Two tuples: [(contrast_start, contrast_end), (offset_start, offset_end)]
Return type:: list[tuple[int, int]]

Notes

Contrast group: indices [0, n_per_group)
Offset group: indices [n_per_group, 2*n_per_group)

get_physical_indices()[source]

Get indices of physical parameters.

Returns:: Indices of physical parameters in the full parameter vector.
Return type:: list[int]

get_per_angle_indices()[source]

Get indices of all per-angle parameters.

Returns:: Indices of per-angle parameters (contrast + offset).
Return type:: list[int]

validate_indices(params)[source]

Validate that group indices are within parameter vector bounds.

Parameters:: params (ndarray) – Full parameter vector.
Returns:: True if all indices are valid, False otherwise.
Return type:: bool
Raises:: ValueError – If indices are out of bounds (with descriptive message).

get_diagnostics()[source]

Get diagnostic information for logging.

Returns:: Diagnostic information including mode, counts, and indices.
Return type:: dict

get_covariance_slice_indices()[source]

Get slice indices for covariance matrix transformation.

Returns slices for extracting per-angle and physical parameter blocks from a covariance matrix.

Returns:: (per_angle_slice, physical_slice) for indexing covariance matrices.
Return type:: tuple[slice, slice]

__init__(n_phi, n_physical, fourier=None, use_constant=False)

Anti-Degeneracy Defense System¶

Introduction¶

Theoretical Background¶

Understanding the Root Causes¶

Root Cause 1: Gradient Cancellation¶

Root Cause 2: Structural Degeneracy¶

Root Cause 3: Ineffective Regularization¶

Layer 1: Fourier Reparameterization¶

Mathematical Foundation¶

Per-Angle Mode Selection Guide¶

Mode Comparison¶

Decision Flowchart¶

Parameter Count Reduction¶

Quantile-Based Scaling Estimation¶

Why Constant Mode Works Best¶

Why Fourier Mode is Still Useful¶

Implementation¶

Basis Matrix Computation¶

Covariance Transformation¶

ParameterIndexMapper Reference¶

Layer 2: Hierarchical Optimization¶

Algorithmic Approach¶

Why This Works¶

Implementation¶

Layer 3: Adaptive CV Regularization¶

From Absolute to Relative Regularization¶

Auto-Tuned Lambda¶

Why CV-Based Regularization¶

Implementation¶

Layer 4: Gradient Collapse Detection¶

Runtime Monitoring¶

Response Actions¶

Implementation¶

Layer 5: Shear-Sensitivity Weighting¶

Mathematical Foundation¶

Why This Works¶

Data Flow¶

Dynamic Weight Updates¶

Implementation¶

Cross-Path Uniformity¶

Configuration Reference¶

YAML Configuration¶

Configuration Options Summary¶

Usage Tutorial¶

Basic Usage¶

Advanced Usage: Fine-Tuning¶

Diagnosing Problems¶

When to Disable¶

API Reference¶

Parameter Utilities¶

Fourier Reparameterization¶

Mathematical Formulation¶

Hierarchical Optimization¶

Why It Works¶

Adaptive Regularization¶

Gradient Monitoring¶

Response Actions¶

Shear-Sensitivity Weighting¶

The Problem¶

The Solution¶

Configuration¶

Parameter Index Mapper¶

See Also¶