Function apply_chisep_operator 

fn apply_chisep_operator(
    ws: &mut ChiSepWorkspace,
    d_kernel: &[f32],
    mx: &[f32],
    my: &[f32],
    mz: &[f32],
    vr_pos: &[f32],
    vr_neg: &[f32],
    vr_sum: &[f32],
    lambda_para: f32,
    lambda_dia: f32,
    lambda_sum: f32,
    lambda_cpl: f32,
    dr_p: f32,
    dr_q: f32,
    mask: &[u8],
    dx: &[f32],
    out: &mut [f32],
)

Apply chi-sep Hessian operator A to doubled vector dx = [d_pos; d_neg].

A_pos = λ_para * TV_hess(d_pos) + λ_sum * TV_hess_sum(d_pos+d_neg) + λ_cpl * D²(d_pos + d_neg) + λ_cpl * dr_p * (dr_p * d_pos - dr_q * d_neg) A_neg = λ_dia * TV_hess(d_neg) + λ_sum * TV_hess_sum(d_pos+d_neg) + λ_cpl * D²(d_pos + d_neg) - λ_cpl * dr_q * (dr_p * d_pos - dr_q * d_neg)

TV_hessian uses IRLS: bdiv(wGVrwG*fgrad(.)), positive semi-definite. Field fidelity: D², same for both, positive semi-definite. R2’: rank-1 structure [dr_p, -dr_q]^T * [dr_p, -dr_q], positive semi-definite.