Add rational P(x)/Q(x) fitting + 4 rough-edge fixes#16
Conversation
|
Hi! I really appreciate all the work you put into this, first of all! Here are my honest thoughts Horner'sYou are correct that this needs to be there. Luckily the upcoming changes included a new method in the Basis trait for this, so I have added implementations for Horner's and Clenshaw (chebyshev) Because Horner's is Monomial specific, adding a module just for that which bypasses all my traits is a massive footgun, so it will need to live in the trait. Chebyshev warningThis is a great point and needs to be documented however, please keep in mind the target audience for polyfit are devs without knowlege of stats, so the explanation will need to be simplified, with the specifics in a 'technical details' block as specified in the contribution guidelines Rational moduleIn theory I love the addition of rational polynomials to polyfit, but unfortunately I cannot merge this implementation. While very solid, and detailed, the API surface is not compatible with the vision I have for polyfit as a whole - a library for non-statisticians to leverage the power of stats. 1. The traitsLike the horners addition, one of my bigger issues is that is isn't leveraging any of polyfit's infrastructure - this could easily be generic over 2. Magic numbersIn my experience, if you ask for a random number, you get a random number. A great example is this tolerance parameter. I would have zero idea how to find a good value for this. I'd prefer that rational did what the rest of polyfit does, provide a custom option for expert use, but wrap the rest in opinionated options, look at DegreeBound for example of what I mean, or better yet, just derive a good value from what info you have (Like I do for the cutoff in Polynomial::spectral_energy_filter 3. The examplesMy examples are designed for non-experts, I try to keep them short, specific, simple to understand, and critically, over documented. Most of them have more comments than code. The main reason is that I like to turn them into tutorial posts for the Polyfit website. The examples around rationals would need to be simplified and split up, named according to use case (for example I called mine 'outlier detection', not 'using a covariance matrix - try to seperate the how from the why) And they need a lot more comments - I was having trouble understanding what they all do |
5 commits ahead of master:
The rational module is the bulk (Sanathanan-Koerner + Levenberg-Marquardt, constraints, 5 examples, ~2,647 LOC in src/rational.rs).
I've been using polyfit for magetypes (SIMD primitive types for the archmage crate, vector-width polynomial evaluators for transcendentals) plus various AI/numerical experiments. The four small follow-up commits address rough edges hit while driving the library — most recently fitting f(x) = x^(1/2.4) on [0, 12] for a SIMD hot loop in zentone:
polyfit::eval module — every caller shipping a polyfit-derived polynomial ends up re-implementing the same Horner-form evaluator. Exposes
horner_f32,horner_f64,horner_f32_via_f64, andrational_horner_f32_via_f64(the last is bit-identical to the internal scorer inRationalPolynomial::optimize_f32, verified by a parity test at 200 sample points).Constraint::soft alias for with_weight — the
with_weight(0.0, 0.0, 1e6)API name suggests a near-hard pin, but the weight is a soft LS penalty and SK iteration drifts (boundary ulp ~8.5e8 even with1e6set). AddsConstraint::softas the clearer name and strengthens the type-level docs to spell out that all weights are soft.with_weightstays as an alias; no API break.as_monomial precision warning —
ChebyshevBasis::as_monomial()silently loses precision past ~degree 12 on wide domains (Vandermonde-style condition number explodes). On[0, 12]√x-substituted fit ofx^(1/2.4): direct Chebyshev evaluation stable past degree 20, monomial afteras_monomialgives 8.0e-3 at deg 12, 1.9e-1 at deg 14. Adds a doc warning with a degree/domain-half-width stable/risky table.F32 ULP metric warning —
F32SearchResult::max_ulpis technically correct but useless as a quality signal when the reference output approaches zero. Onf(x) = x^(1/2.4)every candidate fit reportedmax_ulp ≈ 1e9even whenmax_abs_err < 5e-4. Adds doc warnings onF32ScoreMetric::Ulpand the relevant fields listing trap-case reference functions and a rule of thumb for switching toAbsoluteError.Known limitation
SK iteration is unstable on functions with a singular tangent at the fit-domain boundary (
f'(x_0) → ∞for somex_0 ∈ {a, b}). Concretely, fittingf(x) = x^(1/2.4)on[0, 12]produces non-monotonicmax_absacross(p_deg, q_deg):(3,3) → 2.2e-2,(4,4) → 4.1e-2,(5,4) → 3.7e-2,(6,5) → 5.8e-2. Even withrestarts=8andn_samples=40000the fitter consistently lands on bad local minima. Workaround: sqrt substitution (g(z) = z^(2/2.4)forz = √x) flattens the singular tangent and produces clean monomial fits (3.44e-4 max abs error with a 2-piece monomial split at√0.07). Could be addressed in a follow-up PR via either aRationalFit::from_function_with_transformconstructor or a sqrt-substitution heuristic in the SK fitter when a boundary singular tangent is detected.All gates pass: build, test (123 lib + 118 doc + 2 ignored), fmt, clippy, doc, bench, rdme. No new clippy warnings.
Draft because the rational module is large enough to want a real look before merge.