Function update_hessian 

pub fn update_hessian(
    h_inv: &DMatrix<f64>,
    sk: &DVector<f64>,
    yk: &DVector<f64>,
) -> DMatrix<f64>

Updates the inverse Hessian matrix using the BFGS formula.

This implements the BFGS update for the inverse Hessian (not Hessian), matching the Fortran MECP implementation exactly.

Fortran BFGS Formula (from UpdateX subroutine)

fac = 1 / (DelG · DelX)
fad = 1 / (DelG · H_inv · DelG)
w = fac * DelX - fad * H_inv · DelG
H_inv_new = H_inv + fac * DelX * DelX^T - fad * (H_inv·DelG) * (H_inv·DelG)^T + fae * w * w^T

where:

• DelX = X_new - X_old (step vector, in A)
• DelG = G_new - G_old (gradient difference, in Ha/A)
• fae = DelG · H_inv · DelG

Arguments

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

Returns

Returns the updated inverse Hessian matrix in Ų/Ha. If the update would be unstable, returns the original inverse Hessian.

Units

• Input h_inv: Ų/Ha
• Input sk: A (step vector)
• Input yk: Ha/A (gradient difference)
• Output: Ų/Ha (maintains inverse Hessian units)

Unit Analysis

The BFGS update preserves units:

• fac = 1 / (yk · sk) = 1 / (Ha/A × A) = 1/Ha
• fac * sk * sk^T = 1/Ha × A × A = Ų/Ha ✓
• fad = 1 / (yk · H_inv · yk) = 1 / (Ha/A × Ų/Ha × Ha/A) = A/Ha
• fad * (H_inv·yk) * (H_inv·yk)^T = A/Ha × A × A = A³/Ha (needs fae correction)
• fae * w * w^T corrects to maintain Ų/Ha

Requirements

Validates: Requirements 3.2, 3.3, 3.4