[jessegrabowski]

Kernel is now a graph on row placeholders, and its leaves are typed KernelOp nodes that take their parameters as explicit inputs. Sums and products are ordinary graph arithmetic, and calling a kernel replaces the placeholders. SumKernel, ProductKernel, _eval, and the hand-written diag methods are removed, and Kernel.diag is pt.diagonal(self(X)).

Everything in kernel construction is now a graph node, so kernels can be changed with graph manipulations alone. graph_replace, vectorize_graph, pm.do, pm.observe, and pm.set_data all reach kernel parameters and rows without a kernel-specific path. None of this needs a pm.Potential, because pm.do also reaches a kernel inside a valued pm.MvNormal. Efficient forms are rewrites on the Op types, for example diagonal_of_kernel_to_kernel_diagonal.

This is the kernel counterpart to the symbolic likelihood work in #30, #50, and #68.

The VFE gradient depends on the PyTensor rewrites in pymc-devs/pytensor#2613, which aren't released yet. That PR will move to the index write versions of these rewrites, and the timings below use those versions. On PyTensor 3.3.7, loss plus gradients at N=3200 takes 152 ms against 1.56 ms on main, and it scales quadratically in N. With those rewrites it takes 1.53 ms and scales linearly. Five tests fail on 3.3.7, and all of them pass with those rewrites. Four need the rewrites, and test_dims_only_rv_compiles_under_jax needs a constant arange argument under JAX.

K(X, X) is annotated symmetric and positive definite by declaration, for composites too, as it is on main. -1.0 * K is therefore lowered to Cholesky as if it were positive definite.

---

📚 Documentation preview 📚: https://ptgp--76.org.readthedocs.build/en/76/