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