Create a gist now

Instantly share code, notes, and snippets.

Christoffel-Darboux recurrence relation for orthogonal polynomials using numpy
# this is a translation of orthpoly.ado from Stata 11. Their license likely applies.
def orthpoly(X, deg, weights=None):
"""
Christoffel-Darboux recurrence relation for orthogonal polynomials.
"""
nobs = len(X)
orth_poly = np.ones((nobs, deg+1))
for i in range(1, deg+1):
t = X*orth_poly[:,i-1]**2
b = np.average(t, weights=weights)
t = X*t
a = np.average(t, weights=weights)
if i > 1:
k = i - 2
t = X*orth_poly[:,i-1]*orth_poly[:,i-2]
c = np.average(t, weights=weights)
else:
k = 0
c = 0
a = 1/np.sqrt(a - b**2 - c**2)
b *= -a
c *= a
orth_poly[:,i] = (a*X + b)*orth_poly[:,i-1] - c * orth_poly[:,k]
return orth_poly[:,1:]
@jseabold
Owner

A one-liner free of any license issues, but not numerically sound given the Vandermonde expansion.

poly = lambda x, degree : np.linalg.qr(np.vander(x, degree + 1)[:, ::-1])[0][:, 1:]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment