Calculus7

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:

$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:

$|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}$ :

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

Explicitly,

$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

$|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}$ .

$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

$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

$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:

$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

$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.

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}$ ?

Institutions ranked by the number of AMS Fellows, 2019 edition

Only those with the count of 3 or greater are included. Not counting the deceased. Considering CUNY as a single institution. Source of data.

1. Rutgers The State University of New Jersey New Brunswick : 44
2. University of California, Los Angeles : 38
3. Massachusetts Institute of Technology : 34
4. University of California, Berkeley : 34
5. University of Michigan : 32
6. Cornell University : 26
7. Princeton University : 26
8. University of Illinois, Urbana-Champaign : 24
9. University of Wisconsin, Madison : 24
10. Brown University : 23
11. Stanford University : 22
12. University of Texas at Austin : 22
13. New York University, Courant Institute : 21
14. The City University of New York : 21
15. University of California, San Diego : 21
16. University of Illinois at Chicago : 21
17. University of Chicago : 20
18. University of Washington : 20
19. Stony Brook University : 19
20. Texas A&M University : 19
21. University of California, Santa Barbara : 19
22. University of Minnesota-Twin Cities : 18
23. University of Pennsylvania : 17
24. Duke University : 16
25. Indiana University, Bloomington : 16
26. Purdue University : 16
27. University of Maryland : 16
28. Georgia Institute of Technology : 15
29. Northwestern University : 15
30. Ohio State University, Columbus : 15
31. Pennsylvania State University : 15
32. University of California, Irvine : 13
33. University of Utah : 13
34. Johns Hopkins University, Baltimore : 12
35. University of British Columbia : 12
36. Boston University : 11
37. Harvard University : 11
38. University of Notre Dame : 11
39. University of Toronto : 11
40. Eidgenössische Technische Hochschule Zürich (ETH Zürich) : 10
41. University of North Carolina at Chapel Hill : 10
42. University of Virginia : 10
43. Vanderbilt University : 10
44. Brandeis University : 9
45. University of California, Davis : 9
46. University of Georgia : 9
47. Columbia University : 8
48. Institute for Advanced Study : 8
49. Rice University : 8
50. Tel Aviv University : 8
51. University of Oregon : 8
52. California Institute of Technology : 7
53. Ecole Polytechnique Fédérale de Lausanne (EPFL) : 7
54. Michigan State University : 7
55. Microsoft Research : 7
56. North Carolina State University : 7
57. University of Nebraska-Lincoln : 7
58. University of Oxford : 7
59. University of Southern California : 7
60. Williams College : 7
61. Carnegie Mellon University : 6
62. The Hebrew University of Jerusalem : 6
63. Université Pierre et Marie Curie (Paris VI) : 6
64. University of Arizona : 6
65. University of Rochester : 6
66. Harvey Mudd College : 5
67. Northeastern University : 5
68. Temple University : 5
69. Université Paris-Diderot : 5
70. University of California, Riverside : 5
71. University of Colorado, Boulder : 5
72. Virginia Polytechnic Institute and State University : 5
73. Boston College : 4
74. Florida State University : 4
75. Louisiana State University, Baton Rouge : 4
76. NYU Polytechnic School of Engineering : 4
77. Norwegian University of Science and Technology : 4
78. Rutgers The State University of New Jersey Newark : 4
79. University of Cambridge : 4
80. University of Connecticut, Storrs : 4
81. University of Missouri-Columbia : 4
82. University of Tennessee, Knoxville : 4
83. University of Warwick : 4
84. Washington University : 4
85. Weizmann Institute of Science : 4
86. Yale University : 4
87. Alfréd Rényi Institute of Mathematics, Hungarian Academy of Sciences : 3
88. Australian National University : 3
89. Barnard College, Columbia University : 3
90. Emory University : 3
91. Imperial College : 3
92. Københavns Universitet : 3
93. KTH Royal Institute of Technology : 3
94. Mathematical Institute, Oxford University : 3
95. McGill University : 3
96. Pomona College : 3
97. Shanghai Jiao Tong University : 3
98. Tsinghua University : 3
99. Tufts University : 3
100. Università degli Studi di Milano : 3
101. Université Paris-Sud (Paris XI) : 3
102. Université de Montréal : 3
103. University of Iowa : 3
104. University of Melbourne : 3
105. University of Memphis : 3

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.

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.

Surely the cube root is even worse. Or is it?

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].

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.

Rapid oscillations forces subdivisions here. There are also other reasons for subdivision, such as the loss of analyticity:

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:

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()

Nonexpanding Jordan curves

A simple closed curve ${\Gamma}$ on the plane can be parameterized by a homeomorphism ${f\colon \mathbb T\to \Gamma}$ in infinitely many ways. It is natural to look for “nice” parameterizations: smooth ones, for example. I do not want to require smooth ${\Gamma}$ here, so let us try to find ${f}$ that is nonexpanding, that is ${|f(a)-f(b)|\le |a-b|}$ for all ${a, b\in \mathbb T}$ . Note that Euclidean distance is used here, not arclength.

What are some necessary conditions for the existence of a nonexpanding parameterization?

The curve must have length at most ${2\pi}$ , since nonexpanding maps do not increase length. But this is not sufficient: an equilateral triangle of sidelength ${2\pi/3}$ has no nonexpanding parameterization, despite its length being ${2\pi}$ .
The curve must have diameter at most 2 (which the triangle in item 1 fails). Indeed, nonexpanding maps do not increase the diameter either. However, ${\mathrm{diam}\,\Gamma\le 2}$ is not sufficient either: an equilateral triangle of sidelength ${2}$ has no nonexpanding parameterization, despite its diameter being 2 (and length ${6 < 2\pi}$ ).
The curve must be contained in some closed disk of radius 1. This is not as obvious as the previous two items. We need Kirszbraun’s theorem here: any nonexpanding map ${f\colon \mathbb T\to \Gamma}$ extends to a nonexpanding map ${F\colon \mathbb R^2\to\mathbb R^2}$ , and therefore ${f(\mathbb T)}$ is contained in the closed disk of radius 1 centered at ${F(0)}$ . (This property fails for the triangle in item 2.)

The combination of 1 and 3 (with 2 being superseded by 3) still is not sufficient. A counterexample is given by any polygon that has length ${2\pi}$ but is small enough to fit in a unit disk, for example:

crown — Uneasy lies the head that cannot find a nonexpanding parameterization

Indeed, since the length is exactly ${2\pi}$ , a nonexpanding parametrization must have constant speed 1. But mapping a circular arc onto a line segment with speed 1 increases pairwise Euclidean distances, since we are straightening out the arc.

Since I do not see a way to refine the necessary conditions further, let us turn to the sufficient ones.

It is sufficient for ${\Gamma}$ to be a convex curve contained in the unit disk. Indeed, the nearest-point projection onto a convex set is a nonexpanding map, and projecting the unit circle onto the curve in this way gives the desired parameterization.
It is sufficient for the curve to have length at most 4. Indeed, in this case there exists a parameterization ${f\colon \mathbb T\to\Gamma}$ with ${|f'|\le 4/(2\pi) = 2/\pi}$ . For any ${a, b\in \mathbb T}$ the length of the shorter subarc ${\gamma}$ between ${a}$ and ${b}$ is at most ${(\pi/2)|a-b|}$ . Therefore, the length of ${f(\gamma)}$ is at most ${|a-b|}$ , which implies ${|f(a)-f(b)|\le |a-b|}$ .

Clearly, neither of two sufficient conditions is necessary for the existence of a nonexpanding parameterization. But one can consider “length ${\le 2\pi}$ is necessary, length ${\le 4}$ is sufficient” a reasonable resolution of the problem: up to some constant, the length of a curve can decide the existence one way or the other.

Stay tuned for a post on noncontracting parameterizations…

Modes of convergence on abstract (and not so abstract) sets

In how many ways can a series of real-valued functions on an abstract set converge? Having no measure on the domain eliminates the infinitude of modes of convergence based on integral norms. I can think of five modes of convergence of ${\sum f_n}$ where ${f_n\colon X\to \mathbb R}$ :

(P) Pointwise convergence: ${\sum f_n(x)}$ converges for each ${x\in X}$ .
(U) Uniform convergence: the partial sums ${S_n}$ of the series converge to some function ${f\colon X\to\mathbb R}$ uniformly, i.e., ${\sup |S_n-f| \to 0}$ .
(PA) Pointwise convergence of absolute values: ${\sum |f_n(x)|}$ converges for each ${x\in X}$ .
(UA) Uniform convergence of absolute values: like uniform, but for ${\sum |f_n|}$ .
(M) Weierstrass M-test convergence: ${\sum \sup |f_n| }$ converges.

Implications (all easy): (M) implies (UA), which implies both (U) and (PA). Neither (U) nor (PA) implies the other one, but each of them implies (P).

Perhaps (U) and (PA) deserve an illustration, being incomparable. Let ${X = [0, 1]}$ . The constant functions ${f_n(x)=(-1)^n/n}$ form a series that converges uniformly but not in the sense (PA). In the opposite direction, a series of triangles with height 1 and disjoint supports converges (PA) but not (U).

Notably, the sum of the latter series is not a continuous function. This had to happen: by Dini’s theorem, if a series of continuous functions is (PA)-convergent and its sum is continuous, then it is (UA)-convergent. This “self-improving” property of (PA) convergence will comes up again in the second part of this post.

From abstract sets to normed spaces

In functional analysis, the elements of a normed space ${E}$ can often be usefully interpreted as functions on the unit ball of the dual space ${E^*}$ . Indeed, each ${x\in E}$ induces ${f_x(z) = z(x)}$ for ${z\in E^*}$ . Applying the aforementioned modes of convergence to ${\sum x_n}$ with ${x_n\in E}$ , we arrive at

(P) ⇔ Convergence in the weak topology of E.
(U) ⇔ Convergence in the norm topology of E.
(PA) ⇔ Unconditional convergence in the weak topology of E.
(UA) ⇔ Unconditional convergence in the norm topology of E.
(M) ⇔ Absolute convergence, i.e., ${\sum \|x_n\|}$ converges.

The equivalences (P), (U), (M) are straightforward exercises, but the unconditional convergence merits further discussion. For one thing, there are subtly different approaches to formalizing the concept. Following “Normed Linear Spaces” by M. M. Day, let’s say that a series ${\sum x_n}$ is

(RC) Reordered convergent if there exists ${x}$ such that ${\sum x_{\pi(n)} =x}$ for every bijection ${\pi:\mathbb{N}\to\mathbb{N}}$
(UC) Unordered convergent if there exists ${x}$ such that for every neighborhood ${U}$ of ${x}$ there exists a finite set ${E\subset \mathbb{N}}$ with the property that ${\sum_{n\in F}x_n\in U}$ for every finite set ${F}$ containing ${E}$ .
(SC) Subseries convergent if for every increasing sequence of integers ${(n_k)}$ the series ${\sum x_{n_k}}$ converges.
(BC) Bounded-multiplier convergent if for every bounded sequence of scalars ${(a_n)}$ , the series ${\sum a_n x_n}$ converges.

In a general locally convex space, (BC) ⇒ (SC) ⇒ (UC) ⇔ (RC). The issue with reversing the first two implications is that they involve the existence of a sum for some new series, and if the space lacks completeness, the sum might fail to exist for no good reason. All four properties are equivalent in sequentially complete spaces (those where every Cauchy sequence converges).

Let’s prove that interpretation of (PA) stated above, using the (BC) form of unconditional convergence. Suppose ${\sum x_n}$ converges in the sense (PA), that is for every linear functional ${z}$ the series ${\sum |z(x_n)|}$ converges. Then it’s clear that ${\sum a_n x_n}$ has the same property for any bounded scalar sequence ${(a_n)}$ . That is, (PA) implies bounded-multiplier convergence in the weak topology. Conversely, suppose ${\sum x_n}$ enjoys weak bounded-multiplier convergence and let ${z\in E^*}$ . Multiplying each ${x_n}$ by a suitable unimodular factor ${a_n}$ we can get ${z(a_n x_n) > 0}$ for all ${n}$ . Now the weak convergence of ${\sum a_n x_n}$ yields the pointwise convergence of ${\sum |z(x_n)|}$ .

A theorem of Orlicz, proved in the aforementioned book by Day, says that (SC) convergence in the weak topology of a Banach space is equivalent to (SC) convergence in the norm topology. Thanks to completeness, in the norm topology of a Banach space all forms of unconditional convergence are equivalent. The upshot is that (PA) automatically upgrades to (UA) in the context of the elements of a Banach space being considered as functions on the dual unit ball.

Orthonormal bases formed by translation of sequences

The standard orthonormal basis (ONB) in the Hilbert space ${\ell^2(\mathbb N_0)}$ consists of the vectors
(1, 0, 0, 0, …)
(0, 1, 0, 0, …)
(0, 0, 1, 0, …)
…
Let S be the forward shift operator: ${S(x_0, x_1, \dots) = (0, x_0, x_1, \dots)}$ . The aforementioned ONB is precisely the orbit of the first basis vector ${e_0}$ under the iteration of S. Are there any other vectors x whose orbit under S is an ONB?

If one tries to construct such x by hand, taking some finite linear combination of ${e_k}$ , this is not going to work. Indeed, if the coefficient sequence has finitely many nonzero terms, then one of them, say, ${c_m}$ , is the last one. Then ${S^m x}$ is not orthogonal to ${x}$ because the inner product is ${\overline{c_0} c_m}$ and that is not zero.

However, such vectors x exist, and arise naturally in complex analysis. Indeed, to a sequence ${(c_n)\in\ell^2(\mathbb N_0)}$ we can associate a function ${f(z) = \sum_{n=1}^\infty c_n z^n}$ . The series converges in the open unit disk to a holomorphic function which, being in the Hardy space ${H^2}$ , has boundary values represented by an square-integrable function on the unit circle ${\mathbb T}$ . Forward shift corresponds to multiplication by ${z}$ . Thus, the orthogonality requires that for every ${k\in \mathbb N}$ the function ${z^kf}$ be orthogonal to ${f}$ in ${L^2(\mathbb T)}$ . This means that ${|f|^2}$ is orthogonal to all such ${z^k}$ ; and since it’s real, it is orthogonal to ${z^k}$ for all ${k\in \mathbb Z\setminus \{0\}}$ by virtue of conjugation. Conclusion: |f| has to be constant on the boundary; specifically we need |f|=1 a.e. to have a normalized basis. All the steps can be reversed: |f|=1 is also sufficient for orthogonality.

So, all we need is a holomorphic function f on the unit disk such that almost all boundary values are unimodular and f(0) is nonzero; the latter requirement comes from having to span the entire space. In addition to the constant 1, which yields the canonical ONB, one can use

A Möbius transformation ${f(z)=(z-a)/(1-\bar a z)}$ where ${a\ne 0}$ .
A product of those (a Blaschke product), which can be infinite if the numbers ${a_k}$ converge to the boundary at a sufficient rate to ensure the convergence of the series.
The function ${f(z) = \exp((z+1)/(z-1))}$ which is not a Blaschke product (indeed, it has no zeros) yet satisfies ${|f(z)|=1}$ for all ${z\in \mathbb T\setminus \{-1\}}$ .
Most generally, an inner function which is a Blaschke product multiplied by an integral of rotated versions of the aforementioned exponential function.

Arguably the simplest of these is the Möbius transformation with ${a=1/2}$ ; expanding it into the Taylor series we get
$-\frac12 +\sum_{n=1}^\infty \frac{3}{2^{n+1}} z^n$
Thus, the second simplest ONB-by-translations after the canonical one consists of
(-1/2, 3/4, 3/8, 3/16, 3/32, 3/64, …)
(0, -1/2, 3/4, 3/8, 3/16, 3/32, …)
(0, 0, -1/2, 3/4, 3/8, 3/16, …)
and so on. Direct verification of the ONB properties is an exercise in summing geometric series.

What about the exponential one? The Taylor series of ${\exp((z+1)/(z-1))}$ begins with
$\frac{1}{e}\left(1 - 2z +\frac23 z^3 + \frac23 z^4 +\frac25 z^5 + \frac{4}{45} z^6 - \frac{10}{63} z^7 - \frac{32}{105}z^8 - \cdots \right)$
I don’t know if these coefficients in parentheses have any significance. Well perhaps they do because the sum of their squares is ${e^2}$ . But I don’t know anything else about them. For example, are there infinitely many terms of either sign?

Geometrically, a Möbius transform corresponds to traversing the boundary circle once, a Blaschke product of degree n means doing it n times, while the exponential function, as well as infinite Blaschke products, manage to map a circle onto itself so that it winds around infinitely many times.

Finally, is there anything like that for the backward shift ${S^*(x_0, x_1, x_2, \dots) = (x_1, x_2, \dots)}$ ? The vector ${(S^*)^k x}$ is orthogonal to ${x}$ if and only if ${x}$ is orthogonal to ${S^kx}$ , so the condition for orthogonality is the same as above. But the orbit of any ${\ell^2(\mathbb N_0)}$ vector under ${S^*}$ tends to zero, thus cannot be an orthonormal basis.

Critical points of a cubic spline

The choice of piecewise polynomials of degree 3 for interpolation is justifiably popular: even-degree splines are algebraically awkward to construct, degree 1 is simply piecewise linear interpolation (not smooth), and degree 5, while feasible, entails juggling too many coefficients. Besides, a cubic polynomial minimizes the amount of wiggling (the integral of second derivative squared) for given values and slopes at the endpoints of an interval. (Recall Connecting dots naturally.)

But the derivative of a cubic spline is a quadratic spline. And one needs the derivative to find the critical points. This results in an awkward example in SciPy documentation, annotated with “(NB: sproot only works for order 3 splines, so we fit an order 4 spline)”.

Although not implemented in SciPy, the task of computing the roots of a quadratic spline is a simple one. Obtaining the roots from the internal representation of a quadratic spline in SciPy (as a linear combination of B-splines) would take some work and reading. But a quadratic polynomial is determined by three values, so sampling it at three points, such as two consecutive knots and their average, is enough.

Quadratic formula with values instead of coefficients

Suppose we know the values of a quadratic polynomial q at -1, 0, 1, and wish to find if it has roots between -1 and 1. Let’s normalize so that q(0)=1, and let x = q(-1), y = q(1). If either x or y is negative, there is definitely a root on the interval. If they are positive, there is still a chance: we need the parabola to be concave up, have a minimum within [-1, 1], and for the minimum to be negative. All of this is easily determined once we note that the coefficients of the polynomial are a = (x+y)/2 – 1, b = (y-x)/2, and c = 1.

The inequality ${(x-y)^2 \ge 8(x+y-2)}$ ensures the suitable sign of the discriminant. It describes a parabola with vertex (1, 1) and focus (2, 2), contained in the first quadrant and tangent to the axes at (4, 0) and (0, 4). Within the orange region there are no real roots.

The line x+y=2, tangent to the parabola at its vertex, separates convex and concave parabolas. While concavity in conjunction with x, y being positive definitely precludes having roots in [-1, 1], slight convexity is not much better: it results in real roots outside of the interval. Here is the complete picture: green means there is a root in [-1, 1], orange means no real roots, red covers the rest.

all_colors — Green = there is a root in the interval [-1, 1]

Back to splines

Since the derivative of a spline is implemented in SciPy (B-splines have a nice formula for derivatives), all we need is a root-finding routine for quadratic splines. Here it is, based on the above observations but using built-in NumPy polynomial solver np.roots to avoid dealing with various special cases for the coefficients.

def quadratic_spline_roots(spl):
    roots = []
    knots = spl.get_knots()
    for a, b in zip(knots[:-1], knots[1:]):
        u, v, w = spl(a), spl((a+b)/2), spl(b)
        t = np.roots([u+w-2*v, w-u, 2*v])
        t = t[np.isreal(t) & (np.abs(t) <= 1)]
        roots.extend(t*(b-a)/2 + (b+a)/2)
    return np.array(roots)

A demonstration, which plots the spline (blue), its critical points (red), and original data points (black) as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import InterpolatedUnivariateSpline

x = np.arange(7)
y = np.array([3, 1, 1, 2, 2, 4, 3])
f = InterpolatedUnivariateSpline(x, y, k=3)
crit_pts = quadratic_spline_roots(f.derivative())

t = np.linspace(x[0], x[-1], 500)
plt.plot(t, f(t))
plt.plot(x, y, 'kd')
plt.plot(crit_pts, f(crit_pts), 'ro')
plt.show()