In this talk, we present Gaussian process (GP) frameworks for PDE-related problems. We first discuss bilevel hyperparameter learning for GP solvers in PDEs and inverse problems, since the accuracy, stability, and generalization of kernel- and neural network-based methods depend strongly on hyperparameter choices. We evaluate the approach on nonlinear PDEs and PDE inverse problems, where the results indicate improved accuracy and robustness compared with random initialization. We then turn to GP operator learning for non-perturbative functional renormalization group equations, which are integro-differential equations defined on functionals. The method learns operators directly on function space, producing a functional representation that does not depend on a specific discretization, while still allowing physical priors to be incorporated through the prior mean or kernel design. We illustrate the approach on examples including the Wetterich and Wilson–Polchinski equations.