Function update_hessian_psb 

pub fn update_hessian_psb(
    hessian: &DMatrix<f64>,
    sk: &DVector<f64>,
    yk: &DVector<f64>,
) -> DMatrix<f64>

Updates the Hessian matrix using the PSB (Powell-Symmetric-Broyden) formula.

This is a direct port of HessianUpdator() function. PSB is more appropriate than BFGS for MECP optimization because MECPs have saddle-point-like character on the difference PES.

PSB Formula

v = yk - B·sk
B_new = B + (v·sk^T + sk·v^T) / (sk^T·sk)
      - (sk^T·v) · (sk·sk^T) / (sk^T·sk)²

where:

• sk = x_new - x_old (step vector, in A)
• yk = g_new - g_old (gradient difference, in Ha/A)

Arguments

• hessian - Current Hessian approximation (Ha/A²)
• sk - Step vector (x_new - x_old) in A
• yk - Gradient difference (g_new - g_old) in Ha/A

Returns

Returns the updated Hessian matrix in Ha/A².

Unit Analysis

• v = yk - B·sk → Ha/A - (Ha/A²)(A) = Ha/A ✓
• v·sk^T → (Ha/A)(A) = Ha (matrix outer product) → / A² → Ha/A² ✓
• sk^T·v → A·(Ha/A) = Ha (scalar)
• sk·sk^T / (sk^T·sk)² → A²/A⁴ = 1/A² → × Ha → Ha/A² ✓