"Gaussian processes (GPs) and kernel methods are promising automatic approaches for solving PDEs, as they combine the theoretical rigor of traditional numerical algorithms with the flexible design of machine learning solvers. The complexity bottleneck of GP-based methods lies in computing with dense kernel matrices. In the case of PDE problems, these matrices may also involve partial derivatives of the kernels, and fast algorithms for such matrices are less developed compared to the derivative-free case. In this talk, we will discuss a rigorous sparse Cholesky factorization algorithm to make GP-based PDE solvers scalable. The algorithm relies on the near-sparsity of the Cholesky factor under a multiscale ordering of the pointwise and derivative-type entries of the matrices. It enables us to compute $\epsilon$-approximate inverse Cholesky factors of the kernel matrices with complexity $O(N \log^d (N/\epsilon))$ in space and $O(N\log^{2d}(N/\epsilon))$ in time. As a result, this leads to a near-linear space/time complexity method for solving general PDEs with GPs."