Skip to content

AFT and PO share the PH optimiser's two omissions — measure optimise_nm_tnc before changing it #331

Description

@derrynknife

Split out of #328, which fixed the same two omissions for PH only.

What is shared

optimise_nm_tnc in surpyval/univariate/regression/_fit_skeleton.py serves the AFT and PO families:

def optimise_nm_tnc(fun: Callable, init_t):
    """AFT/PO's historical ladder: Nelder-Mead, then TNC kept only on
    success."""
    res = minimize(
        fun, init_t, method="Nelder-Mead", options={"maxiter": 1000}
    )
    res2 = minimize(fun, res.x, method="TNC")
    return res2 if res2.success else res

Like the old optimise_ph, it passes no jac even though fun closes over regression_neg_ll, which is written in autograd.numpy and is fully traceable. And like the old optimise_ph, the search is not preconditioned, so it inherits the scale sensitivity fixed in #323 and #328: the gradient shrinks like 1/theta with the data magnitude and grows like n with the sample size, and no absolute gtol serves every fit.

What is not shared, and why this is not a copy of #328

It does already guard the TNC rung on success — that was the third defect in optimise_ph and it does not apply here.

More importantly the first rung is Nelder-Mead, which is derivative-free. The specific failure in #328 was BFGS meeting scipy's absolute gtol well short of the optimum on large-magnitude data and reporting success; a derivative-free method does not stop for that reason. So the plausible outcomes here are different:

  • it may be slow but correct, where PH was fast and wrong — Nelder-Mead pays for its robustness in function evaluations, and every evaluation is O(n);
  • it may fail somewhere else entirely, since TNC still polishes the answer and TNC is gradient-based;
  • or the missing gradient may cost little because Nelder-Mead never asks for one.

I have not measured any of this. That is the point of this issue.

Do this first

Run the same A/B I ran for PH before touching the code:

  • Grid over n and p, and a separate sweep over data magnitude (1e-3 to 1e9), which is where PH broke.
  • Compare against lifelines' WeibullAFTFitter and LogNormalAFTFitter. This is a direct comparison — no reparameterisation needed, unlike the PH case which had to go through beta_aft = -gamma_ph / shape.
  • Score both packages' answers on an independently written AFT log-likelihood rather than trusting either one's reported figure. Scoring on surpyval's own objective is what caught the Truncation correction manufactures likelihood when the survival probability underflows #326 misdiagnosis.
  • Proportional odds has no lifelines counterpart, so it needs a different oracle — a fit at unit scale transported to large scale by the exact change of units, as in test_weibull_ph_at_large_scale_reaches_the_optimum.

Only then decide the change. If it looks like PH did, the fix is the same: preconditioned_bfgs (now in fitters/__init__.py) on the autograd gradient, then TNC, then Nelder-Mead as the derivative-free fallback, stopping at the first rung that converges.

Reference

For PH the equivalent change gave 1.4–6x on wall clock and turned a 1.468-nat shortfall at data scale 1e6 into 1.1e-08. There is no reason yet to expect the same numbers here — the ladder is differently shaped.

Test pattern to reuse: test_weibull_ph_fit_is_invariant_to_the_units_of_x and test_weibull_ph_at_large_scale_reaches_the_optimum in surpyval/tests/univariate/regression/test_proportional_hazards.py. Both fail on the pre-#328 PH ladder, so they are known to discriminate.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions