Extreme values of a reproducing kernel for polynomials

For every nonnegative integer {n} there exists a (unique) polynomial {K_n(x, y)} of degree {n} in {x} and {y} separately with the following reproducing property:

{\displaystyle p(x) = \int_{-1}^1 K_n(x, y)p(y)\,dy}

for every polynomial {p} of degree at most {n}, and for every {x}. For example, {K_1(x, y)= (3xy+1)/2}; other examples are found in the post Polynomial delta function.

This fact gives an explicit pointwise bound on a polynomial in terms of its integral on an interval:

{\displaystyle |p(x)| \le M_n(x) \int_{-1}^1 |p(y)|\,dy}

where {M_n(x) = \sup\{|K(x, y)| \colon y\in [-1, 1]\}}. For example, {M_1(x) = (3|x|+1)/2}.

Although in principle {x} could be any real or complex number, it makes sense to restrict attention to {x\in [-1, 1]}, where integration takes place. This leads to the search for extreme values of {K} on the square {Q=[-1, 1]\times [-1, 1]}. Here is how this function looks for {n=1, 2, 3}:

K1
Degree 1
K2
Degree 2
K3
Degree 3

The symmetries {K(x, y)=K(-x, -y) = K(y, x)} are evident here.

Explicitly,

{\displaystyle K_n(x, y) = \sum_{k=0}^n \frac{2k+1}{2} P_k(x)P_k(y)}

where {P_k} is the Legendre polynomial of degree {k} and the factor {(2k+1)/2} is included to make the polynomials an orthonormal set in {L^2(-1, 1)}. Since {P_k} oscillates between {-1} and {1}, it follows that

{\displaystyle |K_n(x, y)|\le \sum_{k=0}^n \frac{2k+1}{2} = \frac{(n+1)^2}{2}}

and this bound is attained at {K(1, 1)=K(-1,-1)=(n+1)^2/2} because {P_k(1)=1} and {P_k(-1)=(-1)^k}.

Is

{\displaystyle K_n(-1, 1) =\sum_{k=0}^n (-1)^k\frac{2k+1}{2} = (-1)^n \frac{n+1}{2}}

the minimum value of {K} on the square {Q}? Certainly not for even {n}. Indeed, differentiating the sum

{\displaystyle S_n(x) = K_n(x, 1) = \sum_{k=0}^n \frac{2k+1}{2} P_k(x)}

with respect to {x} and using {P_k'(-1) =(-1)^{k-1}k(k+1)/2}, we arrive at

{\displaystyle S_n'(-1) = (-1)^{n-1} \frac{n(n^2+3n+2)}{4}}

which is negative if {n} is even, ruling out this point as a minimum.

What about odd {n}, then: is it true that {K_n \ge -(n+1)/2} on the square {Q}?

{n=1}: yes, {K_1(x, y) = (3xy+1)/2 \ge (-3+1)/2 = -1} is clear enough.

{n=3}: the inequality {K_3\ge -2} is also true… at least numerically. It can be simplified to {35(xy)^3 + 9(xy)^2 + 15xy \ge (21x+21y+3)(x^2+y^2)} but I do not see a way forward from there. At least on the boundary of {Q} it can be shown without much work:

{\displaystyle K_3(x, 1) + 2 = \frac{5}{4}(x+1)(7x^2-4x+1)}

The quadratic term has no real roots, which is easy to check.

{n=5}: similar story, the inequality {K_5\ge -3} appears to be true but I can only prove it on the boundary, using

{\displaystyle K_5(x, 1)+3 = \frac{21}{16}(x + 1)(33 x^4 - 18x^3 - 12x^2 + 2x + 3)}

The quartic term has no real roots, which is not so easy to check.

{n=7}: surprisingly, {K_7(4/5, 1) = -2229959/500000} which is about {-4.46}, disproving the conjectural bound {K_7\ge -4}. This fact is not at all obvious from the plot.

K7
Degree 7 kernel

Questions:

  • Is {K_n \ge -Cn} on the square {Q = [-1, 1]\times [-1, 1]} with some universal constant {C}?
  • Is the minimum of {K_n} on {Q} always attained on the boundary of {Q}?
  • Can {M_n(x) = \sup\{|K(x, y)| \colon y\in [-1, 1]\}} be expressed in closed form, at least for small degrees {n}?

Numerical integration visualized

scipy.integrate.quad is a popular method of numerical integration with Python. Let’s see how it chooses the points at which to evaluate the function being integrated. Begin with a simple example, the exponential function.

exp1
exp(x) on [-1, 1]
The blue dots indicate the evaluation points, their y-coordinate being the order of evaluation. So, it begins with x=0, continues with something close to -1, etc. The function, drawn in red, is not to its true vertical scale.

The placement of dots gives away the method of integration: it is the Gauss-Kronrod quadrature with 10 Gauss nodes and 21 Kronrod nodes, abbreviated (G10, K21).  The Gauss nodes are included in the Kronrod nodes, and of course the function is not evaluated there again. The evaluation process is slightly out of order in that 0 is a Kronrod node but comes first, followed by 10 Gauss nodes, followed by 10 remaining Kronrod nodes. The process ends there, as the comparison of G10 and K21 results shows the necessary precision was reached.

The square root on [0, 1] is not such a nice function, so (G10, K21) does not reach the required precision at once. The interval is bisected again and again until it does.

sqrt
sqrt(x) on [0, 4]
  Surely the cube root is even worse. Or is it?

cbrt
cbrt(x) on [-1, 1]
The nodes are symmetric about the midpoint of the interval of integration. So for any odd function on an interval symmetric about 0 we get G10 = 0 and K21 = 0, and the process stops at once. To see the effect of the cube root singularity, one has to use a non-symmetric interval such as [-1, 2].

cbrt12
cbrt(x) on [-1, 2]
That’s still fewer subdivisions than for the square root: cancellation between the left and right neighborhoods of 0 still helps. Let’s look at smooth functions next.

sinxsq
sin(x^2) on [0, 7]
Rapid oscillations forces subdivisions here. There are also other reasons for subdivision, such as the loss of analyticity:

exp_1x
exp(-1/x) on [0, 1]
Although exp(-1/x) is infinitely differentiable on this interval, the fact that it is not analytic at 0 makes it a more difficult integrand that exp(x). Finally, an analytic function with no oscillation which still needs a bunch of subintervals:

runge
1/(1+x^2) on [-5, 5]
This is the standard example used to illustrate the Runge phenomenon. Although we are not interpolating here, numerical integration is also influenced by the function having a small radius of convergence of its Taylor series. After all, G10 can be thought of as degree-9 interpolation of the given function at the Gauss nodes (the zeros of a Legendre polynomial),  with the formula returning the integral of the interpolating polynomial.

The code used to plot these things:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
global eval_points
f = lambda x: np.exp(x)
def f_int(x):
    eval_points.append(x)
    return f(x)

eval_points = []
a, b = -1, 1
quad(f_int, a, b)
n = len(eval_points)
t = np.linspace(a, b, 1000)
y = f(t)
yy = n*(y - y.min())/(y.max() - y.min())
plt.plot(t, yy, 'r') 
plt.plot(eval_points, np.arange(0, n), '.')
plt.show()

Laguerre polynomials under 1

Laguerre polynomials have many neat definitions; I am partial to {\displaystyle L_n(x) = \left(\frac{d}{dx} - 1\right)^n  \frac{x^n}{n!}} because it’s so easy to remember:

  1. Begin with {x^n}
  2. Repeat “subtract the derivative” {n} times
  3. Normalize so the constant term is 1.

For example, for n=3 this process goes as {x^3} to {x^3-3x^2} to {x^3 -6x^2 + 6x} to {x^3-9x^2+18x -6}, which normalizes to {-\frac16x^3+\frac{3}{2}x^2 -3x +1}. This would make a mean exercise on differentiating polynomials: every mistake snowballs throughout the computation.

What would happen if we added the derivative instead? Nothing really new: this change is equivalent to reversing the direction of the x-axis, so we’d end up with {L_n(-x)}. Incidentally, this shows that the polynomial {L_n(-x)} has positive coefficients, which means the behavior of {L_n} for negative {x} is boring: the values go up as {x} becomes more negative. Laguerre polynomials are all about the interval {[0,\infty)} on which they are orthogonal with respect to the weight {\exp(-x)} and therefore change sign often.

L20.png
20 Laguerre polynomials

But when I look at the plot shown above, it’s not the zeros that draw my attention (perhaps because the x-axis is not shown) but the horizontal line {y=1}, the zero-degree polynomial. The coefficients of {L_n} have alternating signs; in particular, {L_n(0)=1} and {L_n'(0)=-n}. So, nonconstant Laguerre polynomials start off with the value of 1 and immediately dive below it. All except the linear one, {L_1(x)=1-x}, eventually recover and reach 1 again (or so it seems; I don’t have a proof).

L30.png
30 Laguerre polynomials

The yellow curve that is the first to cross the blue line is the 5th degree Laguerre polynomial. Let’s see if any of the other polynomials rises about 1 sooner…

L100.png
100 Laguerre polynomials

Still, nobody beats {L_5} (and the second place is held by {L_4}). By the way, the visible expansion of oscillations is approximately exponential; multiplying the polynomials by {\exp(-x/2)} turns the envelopes into horizontal lines:

lag_functions30.png
30 Laguerre polynomials times exp(-x/2)

Back to the crossing of y=1 line. The quantity to study is the smallest positive root of {L_n - 1}, denoted  {r(n)} from now on. (It is the second smallest root overall; as discussed above, this polynomial has a root at x=0 and no negative roots.) For n=2,3,4,5,6, the value of {r(n)} is  {4, 3, 6-2\sqrt{3}, (15-\sqrt{105})/2, 6} which evaluates to 4, 3, 2.536…, 2.377…, and 6 respectively. I got these with Python / SymPy:

from sympy import *
x = Symbol('x')
[Poly(laguerre(n, x) - 1).all_roots()[1] for n in range(2, 7)]

For higher degrees we need numerics. SymPy can still help (applying .evalf() to the roots), but the process gets slow. Switching to NumPy’s roots method speeds things up, but when it indicated than {r(88)}  and a few others are in double digits, I became suspicious…  a closer check showed this was a numerical artifact.

Conjecture: {r(5) \le r(n) \le r(6)} for all {n}. Moreover, {3 < r(n) < 6} when {n \ge 7}.

Here is a closer look at the exceptional polynomials of degrees 3, 4, 5 and 6, with 1 subtracted from each:

l3456

The first local maximum of {L_n} shifts down and to the left as the degree n increases. The degree n=5 is the last for which {L_n} exceeds 1 on the first attempt, so it becomes the quickest to do so. On the other hand, n=6 fails on its first attempt to clear the bar, and its second attempt is later than for any subsequent Laguerre polynomial; so it sets the record for maximal {r(n)}.

Evaluating high-degree Laguerre polynomials is a numerical challenge: adding large terms of alternating signs can reduce accuracy dramatically. Here is a plot of the degree 98 polynomial (minus 1): where is its first positive root?

L98.png
L(98, x) – 1

Fortunately, SymPy can evaluate Laguerre polynomials at rational points using exact arithmetics since the coefficients are rational. For example, when it evaluates the expression laguerre(98, 5) > 1 to True, that’s a (computer-assisted) proof that {r(98) < 5}, which one could in principle "confirm" by computing the same rational value of {L_{98}(5) } by hand (of course, in this situation a human is far less trustworthy than a computer) . Evaluation at the 13 rational points 3, 3.25, 3.5, … , 5.75, 6 is enough to certify that {r(n) < 6} for {n} up to 200 (with the aforementioned exception of {r(6) = 6}).

The lower bounds call for Sturm’s theorem which is more computationally expensive than sampling at a few rational points. SymPy offers a root-counting routine based on this theorem (it counts the roots within the given closed interval):

for n in range(2, 101):
    num_roots = count_roots((laguerre(n,x)-1)/x, 0, 3)
    print('{} roots for n = {}'.format(num_roots, n))

Division by x eliminates the root at 0, so we are left with the root count on (0,3] — which is 1 for n=3,4 and 2 for n=5. The count is zero for all other degrees up to 100, confirming that {r(n) > 3} for {n \ge 6}.

So, the conjecture looks solid. I don’t have a clue to its proof (nor do I know if it’s something known). The only upper bound on {L_n} that I know is Szegő’s {|L_n(x)|\le \exp(x/2)} for {x\ge 0}, which is not helping here.