Proof of the Christoffel–Darboux formula without induction

In this post I prove the Christoffel–Darboux formula without using induction. It seems that often the Christoffel–Darboux formula is proved with induction. However, I find that the proof with induction does not give insight why the Christoffel–Darboux formula is correct. I found the proof below in a paper by Barry Simon.
The Christoffel–Darboux formula is a well-known formula in the theory of orthogonal polynomials. If \( p_k(x) \) with \(k = 0,1, \ldots \) is a set of orthogonal polynomials, normalized such that \begin{equation*} \int_{-\infty}^{+\infty}\!\! dx\ w(x)\ p_k(x) p_l(x) = \delta_{kl} \end{equation*} then the kernel \begin{equation*} K_n(x,y) = \sum_{k=0}^n p_k(x) p_k(y) \end{equation*} has the form \begin{equation}\label{eq:20151111a} K_n(x,y) = \frac{k_n}{k_{n+1}}\frac{p_{n+1} (x) p_n(y) -p_n(x) p_{n+1} (y)}{x-y} \end{equation} with \( k_n \) the leading coefficient of \( p_n (x) \). Formula \eqref{eq:20151111a} is the Christoffel–Darboux formula.

Proof
To make the proof more transparent, I use the bracket notation which is common in quantum mechanics: \( p_n(x) \) is represented by the ket \( | n \rangle \). The convolution with \( K_n \) is by definition \begin{equation*} \int_{-\infty}^{+\infty}\!\! dy\ w(y)\ K_n(x,y) f(y) = \sum_{k=0}^n p_k(x) \int_{-\infty}^{+\infty}\!\! dy\ w(y) \ p_k(y) f(y) \end{equation*} Thus the operator \( K_n: f \mapsto \int_{-\infty}^{+\infty}\!\! dy\ w(y)\ K_n(x,y) f(y) \) has the form \begin{equation}\label{eq:20151111b} K_n = \sum_{k =0}^n | k \rangle \langle k | \end{equation} The essential point in the proof is to calculate the commutator \( [ x, K_n ] \) \begin{align*} [ x, K_n ] f(x) & = x K_n f(x) - K_n x f(x) \\ & = x \int_{-\infty}^{+\infty}\!\! dy\ w(y)\ K_n(x,y) f(y) - \int_{-\infty}^{+\infty}\!\! dy\ w(y)\ K_n(x,y) y f(y)\\ & = \int_{-\infty}^{+\infty}\!\! dy\ w(y) (x - y) \ K_n(x,y) f(y) \end{align*} The Christoffel–Darboux formula is thus equivalent with \begin{equation}\label{eq:20151111c} [ x, K_n ] = \frac{k_n}{k_{n+1}} \Big( | n+1 \rangle \langle n | - | n\rangle \langle n +1| \Big) \end{equation} Formula \eqref{eq:20151111c} is easily calculated based on \eqref{eq:20151111b}. First expand the multiplication operator \( x \) in the basis \( | k \rangle \) \begin{equation}\label{eq:20151111d} x = \sum_{k,l} | k \rangle\ \langle k | x | l \rangle\ \langle l | \end{equation} The matrix elements \( \langle k | x | l \rangle \) form the Jacobi matrix; that is why I use the notation \( J_{k,l} = \langle k | x | l \rangle \). Now I calculate \begin{align*} [ x, K_n ] &= \sum_{\substack{k,l \\ m \le n}} J_{k,l} \Big[ |k\rangle \langle l| \ ,\ | m \rangle \langle m| \Big] \\ &= \sum_{\substack{k,l \\ m \le n}} J_{k,l} \Big( |k\rangle \langle l| m \rangle \langle m| - | m \rangle \langle m|k\rangle \langle l| \Big) \end{align*} using orthonormality gives \begin{equation*} [ x, K_n ] = \sum_{\substack{k \\ l \le n}} J_{k,l} |k\rangle \langle l| - \sum_{\substack{l \\ k \le n}} J_{k,l} |k\rangle \langle l| \end{equation*} Now use the fact that the Jacobi matrix \( J_{k,l} \) is tridiagonal, almost all terms cancel, except two. This gives \begin{equation*} [ x, K_n ] = J_{n+1,n} |n+1 \rangle \langle n| - J_{n,n+1} |n \rangle \langle n+ 1| \end{equation*} The Christoffel–Darboux formula is then proved because it is well-known that \begin{equation*} J_{n+1,n} = J_{n,n+1} = \frac{k_n}{k_{n+1}} \end{equation*}