#344 imhof(): "Result too large" (error 32) beyond about 2,500-3,000 eigenvalues; suggest computing on the log scale.

[Kirill Naralenkov]

gretl 2026d, Windows 11 Pro 64-bit.

imhof() stops when the quadratic form has a few thousand eigenvalues: it prints "imhof: Result too large" and returns error 32, "The convergence criterion was not met". gretl reports the failure, so nothing goes wrong silently; but exact Durbin–Watson bounds for a panel need about N(T − 1) eigenvalues, for example 9,000 at N = 1000, T = 10.

The attached imhof-many-eigenvalues.inp computes P(d < 1.97) for the Durbin–Watson statistic of a regression on a constant alone, whose T − 1 nonzero roots are 2 − 2cos(πj/T), as imhof(roots - 1.97, 0):

roots   imhof()          imhof_log()       gap
 1000   0.317491598180   0.317491598180    9.8e-15
 2000   0.251075254070   0.251075254070    5.0e-15
 2800   0.213602724094   0.213602724094    3.9e-16
 2900   0.209536884316   0.209536884316    7.7e-15
 3000   error 32         0.205585239916
 5000   error 32         0.144372445977
 9000   error 32         0.077335670637

The limit depends on the spectrum. Bhargava, Franzini & Narendranathan's (1982, Table I) lower bound for T = 6, H = 500, n = 3, printed as 1.9393, uses 2,497 roots: 4 sin²(πq/12) for q = 1, ..., 5, each 500 times, less 3 of the largest. With them imhof(roots - 1.94, 0) already fails; the answer is 0.0527.

The message suggests an overflow. In Imhof's (1961) integral the integrand has the factor 1/ρ(u), with ρ(u) = Π_j (1 + λ_j² u²)^(1/4). With the 3,000 roots above, ρ(u) exceeds the largest double from u ≈ 0.99 on, where 1/(u ρ(u)) is about 10⁻³⁰⁸, while log ρ(u) stays small. The attached imhof_log() is Imhof's integral in hansl with ρ(u) kept on the log scale, by Simpson's rule up to the point where log u + log ρ(u) = 40. It matches imhof() to about 10⁻¹⁴ wherever imhof() works, and gives the values above where it does not. An adaptive quadrature in Python (scipy's quad) on the same integrand gives the same values to the 12 decimals shown.

The suggestion: inside imhof(), compute ρ(u), and anything else that multiplies over the roots, on the log scale. imhof_log() is only a demonstration: its fixed grid suits many roots, where the integrand dies out fast, not a handful.

[Allin Cottrell]

Thanks for the report, Kirill. It turns out that the overflow was basically spurious, in that we were already cumulating rho(u) on the log scale and it was just that the final exponentiation was misplaced. Fixed in git. Snapshots to follow.