OrthogonalPolynomial is unreasonably slow

[nvtroshkin]

OrthogonalPolynomial uses spectral Lanczos algorithm that operates with matrices of ncap size (calculates ncap coefficients) while we are interested only in first n polynomials. Because usually ncap >> n, the overhead is huge. Also the Lanczos algorithm may not converge (#6).

The Steiltjes algorithm for the fejer quadrature may by preferable to use Gautschi, 1993, because it is stable and fast.

Benchmark

_orthpol.dlancz

import numpy as np
import orthpol._orthpol as orth
import timeit

n=20
ncap = 100000
quad = orth.dfejer(ncap)

print(timeit.timeit(lambda: orth.dlancz(n,quad[0], np.array(quad[1]) * 1/ncap), number=5))

_orthpol.dsti

print(timeit.timeit(lambda: orth.dsti(n,quad[0], np.array(quad[1]) * 1/ncap), number=5))

Results

_orthpol.dlancz:

89.07238823200169

_orthpol.dsti:

0.07561909200012451

Discussion

The execution time of the Stieljes algorithm depends linearly on n. For the Lanczos algorithm there is no dependency on n, and we could think that for n >> 1 the Lanczos algorithm is preferable. Actually, it is not the case because of the recurrence coefficients undergo saturation (i.e. alpha[n+1] is almost equal to alpha[n]) and the underlying Fortran code exits with an error, indicating the maximum recurrence index, before the advantage could realize.