[Paper Review] Sparse Cholesky Factorization for Solving Nonlinear PDEs via Gaussian Processes

Yifan Chen, Houman Owhadi+1 more|arXiv (Cornell University)|Apr 3, 2023

Gaussian Processes and Bayesian Inference4 citations

TL;DR

This paper introduces a sparse Cholesky factorization algorithm for Gaussian process (GP)-based solution of nonlinear PDEs, leveraging a novel ordering of Dirac and derivative measurements to induce near-sparsity in the inverse kernel matrix. The method achieves $O(N/log^{d}(N/ ho))$ space and $O(N/log^{2d}(N/ ho))$ time complexity, enabling scalable, high-accuracy solution of nonlinear PDEs like elliptic, Burgers’, and Monge-Ampère equations via efficient optimization with preconditioned conjugate gradient methods.

ABSTRACT

In recent years, there has been widespread adoption of machine learning-based approaches to automate the solving of partial differential equations (PDEs). Among these approaches, Gaussian processes (GPs) and kernel methods have garnered considerable interest due to their flexibility, robust theoretical guarantees, and close ties to traditional methods. They can transform the solving of general nonlinear PDEs into solving quadratic optimization problems with nonlinear, PDE-induced constraints. However, the complexity bottleneck lies in computing with dense kernel matrices obtained from pointwise evaluations of the covariance kernel, and its extit{partial derivatives}, a result of the PDE constraint and for which fast algorithms are scarce. The primary goal of this paper is to provide a near-linear complexity algorithm for working with such kernel matrices. We present a sparse Cholesky factorization algorithm for these matrices based on the near-sparsity of the Cholesky factor under a novel ordering of pointwise and derivative measurements. The near-sparsity is rigorously justified by directly connecting the factor to GP regression and exponential decay of basis functions in numerical homogenization. We then employ the Vecchia approximation of GPs, which is optimal in the Kullback-Leibler divergence, to compute the approximate factor. This enables us to compute $ε$-approximate inverse Cholesky factors of the kernel matrices with complexity $O(N\log^d(N/ε))$ in space and $O(N\log^{2d}(N/ε))$ in time. We integrate sparse Cholesky factorizations into optimization algorithms to obtain fast solvers of the nonlinear PDE. We numerically illustrate our algorithm's near-linear space/time complexity for a broad class of nonlinear PDEs such as the nonlinear elliptic, Burgers, and Monge-Ampère equations.

Motivation & Objective

• To address the computational bottleneck of dense kernel matrices in GP-based PDE solvers, which scale as $O(N^3)$ for $N$ collocation points.
• To develop a scalable, accurate, and efficient algorithm for solving general nonlinear PDEs using Gaussian processes.
• To rigorously analyze the sparsity and approximation accuracy of the Cholesky factor of the inverse kernel matrix under a new measurement ordering.
• To extend the framework to Gauss-Newton methods for nonlinear PDEs, where reduced kernel matrices lack sufficient Dirac measurements.
• To demonstrate near-linear complexity and high accuracy across diverse nonlinear PDEs, including fully nonlinear ones like Monge-Ampère.

Proposed method

• Propose a novel ordering of Dirac and derivative measurements to induce near-sparsity in the Cholesky factor of the inverse kernel matrix.
• Apply Vecchia approximation to the GP, proven optimal in Kullback-Leibler divergence, to compute sparse approximations of the inverse Cholesky factor.
• Use KL minimization to compute the sparse Cholesky factor, with rigorous error bounds showing exponential convergence in the Frobenius norm.
• For reduced kernel matrices in Gauss-Newton iterations, apply the same sparse factor as a preconditioner in conjugate gradient (pCG) solvers.
• Leverage auto-differentiation for kernel derivatives in complex PDEs like Monge-Ampère, ensuring flexibility and accuracy.
• Combine the sparse factorization with sequential quadratic programming and Gauss-Newton methods to achieve scalable, high-accuracy PDE solutions.

Experimental results

Research questions

• RQ1 Can a novel measurement ordering induce near-sparsity in the Cholesky factor of the inverse kernel matrix for GP-based PDE solvers?
• RQ2 What is the theoretical accuracy of the resulting Vecchia approximation, and how does it relate to the Kullback-Leibler divergence?
• RQ3 Can the sparse Cholesky factorization achieve $O(N\log^{d}(N/\epsilon))$ space and time complexity for general nonlinear PDEs?
• RQ4 How does the method perform when applied to reduced kernel matrices in Gauss-Newton iterations, where sparsity decay is weaker?
• RQ5 Does the algorithm maintain near-linear complexity and high accuracy across diverse nonlinear PDEs, including fully nonlinear ones like Monge-Ampère?

Key findings

• The proposed sparse Cholesky factorization achieves $O(N\log^{d}(N/\epsilon))$ space and $O(N\log^{2d}(N/\epsilon))$ time complexity, with $\epsilon$-accuracy in the Frobenius norm.
• The Vecchia approximation based on the new ordering is optimal in Kullback-Leibler divergence and exhibits exponentially convergent accuracy.
• For the Burgers’ equation, the method achieves $L^2$ errors of $6.111 \times 10^{-5}$ at $N=2000$ and $7.453 \times 10^{-5}$ at $N=4000$, indicating high accuracy with fixed parameters.
• In the Monge-Ampère equation example, $L^2$ errors decrease with increasing $N_{\text{domain}}$, confirming convergence despite complex derivative structure.
• The Gauss-Newton method converges in $O(1)$ steps, and pCG with the sparse factor as preconditioner requires nearly constant iterations, enabling near-linear overall complexity.
• The algorithm demonstrates near-linear scaling in CPU time across all tested PDEs, including those requiring auto-differentiation of kernel derivatives.