Search

We design two classes of ultra-fast meta-solvers for linear systems arising after discretizing PDEs by combining neural operators with either simple iterative solvers, e.g., Jacobi and Gauss-Seidel, or with Krylov methods, e.g., GMRES and BiCGStab, using the trunk basis of DeepONet as a coarse preconditioner. The idea is to leverage the spectral bias of neural networks to account for the lower part of the spectrum in the error distribution while the upper part is handled easily and inexpensively using relaxation methods or fine-scale preconditioners. We create a pareto front of optimal meta-solvers using a plurarilty of metrics, and we introduce a preference function to select the best solver most suitable for a specific scenario. This automation for finding optimal solvers can be extended to nonlinear systems and other setups, e.g. finding the best meta-solver for space-time in time-dependent PDEs.

The interplay between stochastic processes and optimal control has been extensively explored in the literature. With the recent surge in the use of diffusion models, stochastic processes have increasingly been applied to sample generation. This paper builds on the log transform, known as the Cole-Hopf transform in Brownian motion contexts, and extends it within a more abstract framework that includes a linear operator. Within this framework, we found that the well-known relationship between the Cole-Hopf transform and optimal transport is a particular instance where the linear operator acts as the infinitesimal generator of a stochastic process. We also introduce a novel scenario where the linear operator is the adjoint of the generator, linking to Bayesian inference under specific initial and terminal conditions. Leveraging this theoretical foundation, we develop a new algorithm, named the HJ-sampler, for Bayesian inference for the inverse problem of a stochastic differential equation with given terminal observations. The HJ-sampler involves two stages: (1) solving the viscous Hamilton-Jacobi partial differential equations, and (2) sampling from the associated stochastic optimal control problem. Our proposed algorithm naturally allows for flexibility in selecting the numerical solver for viscous HJ PDEs. We introduce two variants of the solver: the Riccati-HJ-sampler, based on the Riccati method, and the SGM-HJ-sampler, which utilizes diffusion models. We demonstrate the effectiveness and flexibility of the proposed methods by applying them to solve Bayesian inverse problems involving various stochastic processes and prior distributions, including applications that address model misspecifications and quantifying model uncertainty.

We propose a novel neural network architecture (TSympOCNet) to address high--dimensional optimal control problems with linear and nonlinear dynamics. An important application of this method is to solve the path planning problem of multi-agent vehicles in real time. The new method extends our previous SympOCNet framework by introducing a time-dependent symplectic network into the architecture. In addition, we propose a more general latent representation, which greatly improves model expressivity based on the universal approximation theorem. We demonstrate the efficacy of TSympOCNet in path planning problems with obstacle and collision avoidance, including systems with Newtonian dynamics and non-convex environments, up to dimension 512. Our method shows significant promise in handling efficiently both complex dynamics and constraints.

Uncertainty quantification (UQ) in scientific machine learning (SciML) combines the powerful predictive power of SciML with methods for quantifying the reliability of the learned models. However, two major challenges remain: limited interpretability and expensive training procedures. We provide a new interpretation for UQ problems by establishing a new theoretical connection between some Bayesian inference problems arising in SciML and viscous Hamilton-Jacobi partial differential equations (HJ PDEs). Namely, we show that the posterior mean and covariance can be recovered from the spatial gradient and Hessian of the solution to a viscous HJ PDE. As a first exploration of this connection, we specialize to Bayesian inference problems with linear models, Gaussian likelihoods, and Gaussian priors. In this case, the associated viscous HJ PDEs can be solved using Riccati ODEs, and we develop a new Riccati-based methodology that provides computational advantages when continuously updating the model predictions. Specifically, our Riccati-based approach can efficiently add or remove data points to the training set invariant to the order of the data and continuously tune hyperparameters. Moreover, neither update requires retraining on or access to previously incorporated data. We provide several examples from SciML involving noisy data and \textit{epistemic uncertainty} to illustrate the potential advantages of our approach. In particular, this approach's amenability to data streaming applications demonstrates its potential for real-time inferences, which, in turn, allows for applications in which the predicted uncertainty is used to dynamically alter the learning process.

Maximum entropy (Maxent) models are a class of statistical models that use the maximum entropy principle to estimate probability distributions from data. Due to the size of modern data sets, Maxent models need efficient optimization algorithms to scale well for big data applications. State-of-the-art algorithms for Maxent models, however, were not originally designed to handle big data sets; these algorithms either rely on technical devices that may yield unreliable numerical results, scale poorly, or require smoothness assumptions that many practical Maxent models lack. In this paper, we present novel optimization algorithms that overcome the shortcomings of state-of-the-art algorithms for training large-scale, non-smooth Maxent models. Our proposed first-order algorithms leverage the Kullback-Leibler divergence to train large-scale and non-smooth Maxent models efficiently. For Maxent models with discrete probability distribution of $n$ elements built from samples, each containing $m$ features, the stepsize parameters estimation and iterations in our algorithms scale on the order of $O(mn)$ operations and can be trivially parallelized. Moreover, the strong $\ell_{1}$ convexity of the Kullback--Leibler divergence allows for larger stepsize parameters, thereby speeding up the convergence rate of our algorithms. To illustrate the efficiency of our novel algorithms, we consider the problem of estimating probabilities of fire occurrences as a function of ecological features in the Western US MTBS-Interagency wildfire data set. Our numerical results show that our algorithms outperform the state of the arts by one order of magnitude and yield results that agree with physical models of wildfire occurrence and previous statistical analyses of wildfire drivers., Comment: The main text of our manuscript is 20 pages long, the appendices are 4 pages long, and the references are 4 pages long,for a total of 28 pages

We address two major challenges in scientific machine learning (SciML): interpretability and computational efficiency. We increase the interpretability of certain learning processes by establishing a new theoretical connection between optimization problems arising from SciML and a generalized Hopf formula, which represents the viscosity solution to a Hamilton-Jacobi partial differential equation (HJ PDE) with time-dependent Hamiltonian. Namely, we show that when we solve certain regularized learning problems with integral-type losses, we actually solve an optimal control problem and its associated HJ PDE with time-dependent Hamiltonian. This connection allows us to reinterpret incremental updates to learned models as the evolution of an associated HJ PDE and optimal control problem in time, where all of the previous information is intrinsically encoded in the solution to the HJ PDE. As a result, existing HJ PDE solvers and optimal control algorithms can be reused to design new efficient training approaches for SciML that naturally coincide with the continual learning framework, while avoiding catastrophic forgetting. As a first exploration of this connection, we consider the special case of linear regression and leverage our connection to develop a new Riccati-based methodology for solving these learning problems that is amenable to continual learning applications. We also provide some corresponding numerical examples that demonstrate the potential computational and memory advantages our Riccati-based approach can provide.

Hamilton-Jacobi partial differential equations (HJ PDEs) have deep connections with a wide range of fields, including optimal control, differential games, and imaging sciences. By considering the time variable to be a higher dimensional quantity, HJ PDEs can be extended to the multi-time case. In this paper, we establish a novel theoretical connection between specific optimization problems arising in machine learning and the multi-time Hopf formula, which corresponds to a representation of the solution to certain multi-time HJ PDEs. Through this connection, we increase the interpretability of the training process of certain machine learning applications by showing that when we solve these learning problems, we also solve a multi-time HJ PDE and, by extension, its corresponding optimal control problem. As a first exploration of this connection, we develop the relation between the regularized linear regression problem and the Linear Quadratic Regulator (LQR). We then leverage our theoretical connection to adapt standard LQR solvers (namely, those based on the Riccati ordinary differential equations) to design new training approaches for machine learning. Finally, we provide some numerical examples that demonstrate the versatility and possible computational advantages of our Riccati-based approach in the context of continual learning, post-training calibration, transfer learning, and sparse dynamics identification.

Solving high-dimensional optimal control problems in real-time is an important but challenging problem, with applications to multi-agent path planning problems, which have drawn increased attention given the growing popularity of drones in recent years. In this paper, we propose a novel neural network method called SympOCnet that applies the Symplectic network to solve high-dimensional optimal control problems with state constraints. We present several numerical results on path planning problems in two-dimensional and three-dimensional spaces. Specifically, we demonstrate that our SympOCnet can solve a problem with more than 500 dimensions in 1.5 hours on a single GPU, which shows the effectiveness and efficiency of SympOCnet. The proposed method is scalable and has the potential to solve truly high-dimensional path planning problems in real-time.

Logistic regression is a widely used statistical model to describe the relationship between a binary response variable and predictor variables in data sets. It is often used in machine learning to identify important predictor variables. This task, variable selection, typically amounts to fitting a logistic regression model regularized by a convex combination of $\ell_1$ and $\ell_{2}^{2}$ penalties. Since modern big data sets can contain hundreds of thousands to billions of predictor variables, variable selection methods depend on efficient and robust optimization algorithms to perform well. State-of-the-art algorithms for variable selection, however, were not traditionally designed to handle big data sets; they either scale poorly in size or are prone to produce unreliable numerical results. It therefore remains challenging to perform variable selection on big data sets without access to adequate and costly computational resources. In this paper, we propose a nonlinear primal-dual algorithm that addresses these shortcomings. Specifically, we propose an iterative algorithm that provably computes a solution to a logistic regression problem regularized by an elastic net penalty in $O(T(m,n)\log(1/\epsilon))$ operations, where $\epsilon \in (0,1)$ denotes the tolerance and $T(m,n)$ denotes the number of arithmetic operations required to perform matrix-vector multiplication on a data set with $m$ samples each comprising $n$ features. This result improves on the known complexity bound of $O(\min(m^2n,mn^2)\log(1/\epsilon))$ for first-order optimization methods such as the classic primal-dual hybrid gradient or forward-backward splitting methods., Comment: 15 pages

Two key challenges in optimal control include efficiently solving high-dimensional problems and handling optimal control problems with state-dependent running costs. In this paper, we consider a class of optimal control problems whose running costs consist of a quadratic on the control variable and a convex, non-negative, piecewise affine function on the state variable. We provide the analytical solution for this class of optimal control problems as well as a Hopf-type representation formula for the corresponding Hamilton-Jacobi partial differential equations. Finally, we propose efficient numerical algorithms based on our Hopf-type representation formula, convex optimization algorithms, and min-plus techniques. We present several high-dimensional numerical examples, which demonstrate that our algorithms overcome the curse of dimensionality. We also describe a field-programmable gate array (FPGA) implementation of our numerical solver whose latency scales linearly in the spatial dimension and that achieves approximately a 40 times speedup compared to a parallelized central processing unit (CPU) implementation. Thus, our numerical results demonstrate the promising performance boosts that FPGAs are able to achieve over CPUs. As such, our proposed methods have the potential to serve as a building block for solving more complicated high-dimensional optimal control problems in real-time., Comment: arXiv admin note: text overlap with arXiv:2109.14849

Two of the main challenges in optimal control are solving problems with state-dependent running costs and developing efficient numerical solvers that are computationally tractable in high dimension. In this paper, we provide analytical solutions to certain optimal control problems whose running cost depends on the state variable and with constraints on the control. We also provide Lax-Oleinik-type representation formulas for the corresponding Hamilton-Jacobi partial differential equations with state-dependent Hamiltonians. Additionally, we present an efficient, grid-free numerical solver based on our representation formulas, which is shown to scale linearly with the state dimension, and thus, to overcome the curse of dimensionality. Using existing optimization methods and the min-plus technique, we extend our numerical solvers to address more general classes of convex and nonconvex initial costs. We demonstrate the capabilities of our numerical solvers using implementations on a central processing unit (CPU) and a field-programmable gate array (FPGA). In several cases, our FPGA implementation obtains over a 10 times speedup compared to the CPU, which demonstrates the promising performance boosts FPGAs can achieve. Our numerical results show that our solvers have the potential to serve as a building block for solving broader classes of high-dimensional optimal control problems in real-time.

The linear primal-dual hybrid gradient (PDHG) method is a first-order method that splits convex optimization problems with saddle-point structure into smaller subproblems. Unlike those obtained in most splitting methods, these subproblems can generally be solved efficiently because they involve simple operations such as matrix-vector multiplications or proximal mappings that are fast to evaluate numerically. This advantage comes at the price that the linear PDHG method requires precise stepsize parameters for the problem at hand to achieve an optimal convergence rate. Unfortunately, these stepsize parameters are often prohibitively expensive to compute for large-scale optimization problems, such as those in machine learning. This issue makes the otherwise simple linear PDHG method unsuitable for such problems, and it is also shared by most first-order optimization methods as well. To address this issue, we introduce accelerated nonlinear PDHG methods that achieve an optimal convergence rate with stepsize parameters that are simple and efficient to compute. We prove rigorous convergence results, including results for strongly convex or smooth problems posed on infinite-dimensional reflexive Banach spaces. We illustrate the efficiency of our methods on $\ell_{1}$-constrained logistic regression and entropy-regularized matrix games. Our numerical experiments show that the nonlinear PDHG methods are considerably faster than competing methods., Comment: 52 pages, no figures

We consider image denoising problems formulated as variational problems. It is known that Hamilton-Jacobi PDEs govern the solution of such optimization problems when the noise model is additive. In this work, we address certain non-additive noise models and show that they are also related to Hamilton-Jacobi PDEs. These findings allow us to establish new connections between additive and non-additive noise imaging models. Specifically, we study how the solutions to these optimization problems depend on the parameters and the observed images. We show that the optimal values are ruled by some Hamilton-Jacobi PDEs, while the optimizers are characterized by the spatial gradient of the solution to the Hamilton-Jacobi PDEs. Moreover, we use these relations to investigate the asymptotic behavior of the variational model as the parameter goes to infinity, that is, when the influence of the noise vanishes. With these connections, some non-convex models for non-additive noise can be solved by applying convex optimization algorithms to the equivalent convex models for additive noise. Several numerical results are provided for denoising problems with Poisson noise or multiplicative noise.

Solving high dimensional optimal control problems and corresponding Hamilton-Jacobi PDEs are important but challenging problems in control engineering. In this paper, we propose two abstract neural network architectures which are respectively used to compute the value function and the optimal control for certain class of high dimensional optimal control problems. We provide the mathematical analysis for the two abstract architectures. We also show several numerical results computed using the deep neural network implementations of these abstract architectures. A preliminary implementation of our proposed neural network architecture on FPGAs shows promising speed up compared to CPUs. This work paves the way to leverage efficient dedicated hardware designed for neural networks to solve high dimensional optimal control problems and Hamilton-Jacobi PDEs.

Many imaging problems can be formulated as inverse problems expressed as finite-dimensional optimization problems. These optimization problems generally consist of minimizing the sum of a data fidelity and regularization terms. In [23,26], connections between these optimization problems and (multi-time) Hamilton--Jacobi partial differential equations have been proposed under the convexity assumptions of both the data fidelity and regularization terms. In particular, under these convexity assumptions, some representation formulas for a minimizer can be obtained. From a Bayesian perspective, such a minimizer can be seen as a maximum a posteriori estimator. In this chapter, we consider a certain class of non-convex regularizations and show that similar representation formulas for the minimizer can also be obtained. This is achieved by leveraging min-plus algebra techniques that have been originally developed for solving certain Hamilton--Jacobi partial differential equations arising in optimal control. Note that connections between viscous Hamilton--Jacobi partial differential equations and Bayesian posterior mean estimators with Gaussian data fidelity terms and log-concave priors have been highlighted in [25]. We also present similar results for certain Bayesian posterior mean estimators with Gaussian data fidelity and certain non-log-concave priors using an analogue of min-plus algebra techniques.

We propose a novel Caputo fractional derivative-based optimization algorithm. Upon defining the Caputo fractional gradient with respect to the Cartesian coordinate, we present a generic Caputo fractional gradient descent (CFGD) method. We prove that the CFGD yields the steepest descent direction of a locally smoothed objective function. The generic CFGD requires three parameters to be specified, and a choice of the parameters yields a version of CFGD. We propose three versions -- non-adaptive, adaptive terminal and adaptive order. By focusing on quadratic objective functions, we provide a convergence analysis. We prove that the non-adaptive CFGD converges to a Tikhonov regularized solution. For the two adaptive versions, we derive error bounds, which show convergence to integer-order stationary point under some conditions. We derive an explicit formula of CFGD for quadratic functions. We computationally found that the adaptive terminal (AT) CFGD mitigates the dependence on the condition number in the rate of convergence and results in significant acceleration over gradient descent (GD). For non-quadratic functions, we develop an efficient implementation of CFGD using the Gauss-Jacobi quadrature, whose computational cost is approximately proportional to the number of the quadrature points and the cost of GD. Our numerical examples show that AT-CFGD results in acceleration over GD, even when a small number of the Gauss-Jacobi quadrature points (including a single point) is used.

We consider the problem of optimizing the trajectory of an Unmanned Aerial Vehicle (UAV). Assuming a traffic intensity map of users to be served, the UAV must travel from a given initial location to a final position within a given duration and serves the traffic on its way. The problem consists in finding the optimal trajectory that minimizes a certain cost depending on the velocity and on the amount of served traffic. We formulate the problem using the framework of Lagrangian mechanics. We derive closed-form formulas for the optimal trajectory when the traffic intensity is quadratic (single-phase) using Hamilton-Jacobi equations. When the traffic intensity is bi-phase, i.e. made of two quadratics, we provide necessary conditions of optimality that allow us to propose a gradient-based algorithm and a new algorithm based on the linear control properties of the quadratic model. These two solutions are of very low complexity because they rely on fast convergence numerical schemes and closed form formulas. These two approaches return a trajectory satisfying the necessary conditions of optimality. At last, we propose a data processing procedure based on a modified K-means algorithm to derive a bi-phase model and an optimal trajectory simulation from real traffic data., Comment: 30 pages, 10 figures, 2 tables. arXiv admin note: substantial text overlap with arXiv:1812.08759

Physics informed neural networks (PINNs) are deep learning based techniques for solving partial differential equations (PDEs) encounted in computational science and engineering. Guided by data and physical laws, PINNs find a neural network that approximates the solution to a system of PDEs. Such a neural network is obtained by minimizing a loss function in which any prior knowledge of PDEs and data are encoded. Despite its remarkable empirical success in one, two or three dimensional problems, there is little theoretical justification for PINNs. As the number of data grows, PINNs generate a sequence of minimizers which correspond to a sequence of neural networks. We want to answer the question: Does the sequence of minimizers converge to the solution to the PDE? We consider two classes of PDEs: linear second-order elliptic and parabolic. By adapting the Schauder approach and the maximum principle, we show that the sequence of minimizers strongly converges to the PDE solution in $C^0$. Furthermore, we show that if each minimizer satisfies the initial/boundary conditions, the convergence mode becomes $H^1$. Computational examples are provided to illustrate our theoretical findings. To the best of our knowledge, this is the first theoretical work that shows the consistency of PINNs.

Variational and Bayesian methods are two approaches that have been widely used to solve image reconstruction problems. In this paper, we propose original connections between Hamilton--Jacobi (HJ) partial differential equations and a broad class of Bayesian methods and posterior mean estimators with Gaussian data fidelity term and log-concave prior. Whereas solutions to certain first-order HJ PDEs with initial data describe maximum a posteriori estimators in a Bayesian setting, here we show that solutions to some viscous HJ PDEs with initial data describe a broad class of posterior mean estimators. These connections allow us to establish several representation formulas and optimal bounds involving the posterior mean estimate. In particular, we use these connections to HJ PDEs to show that some Bayesian posterior mean estimators can be expressed as proximal mappings of twice continuously differentiable functions, and furthermore we derive a representation formula for these functions., Comment: 33 pages, 2 figures

We propose novel connections between several neural network architectures and viscosity solutions of some Hamilton--Jacobi (HJ) partial differential equations (PDEs) whose Hamiltonian is convex and only depends on the spatial gradient of the solution. To be specific, we prove that under certain assumptions, the two neural network architectures we proposed represent viscosity solutions to two sets of HJ PDEs with zero error. We also implement our proposed neural network architectures using Tensorflow and provide several examples and illustrations. Note that these neural network representations can avoid curve of dimensionality for certain HJ PDEs, since they do not involve neither grids nor discretization. Our results suggest that efficient dedicated hardware implementation for neural networks can be leveraged to evaluate viscosity solutions of certain HJ PDEs.

We propose new and original mathematical connections between Hamilton-Jacobi (HJ) partial differential equations (PDEs) with initial data and neural network architectures. Specifically, we prove that some classes of neural networks correspond to representation formulas of HJ PDE solutions whose Hamiltonians and initial data are obtained from the parameters of the neural networks. These results do not rely on universal approximation properties of neural networks; rather, our results show that some classes of neural network architectures naturally encode the physics contained in some HJ PDEs. Our results naturally yield efficient neural network-based methods for evaluating solutions of some HJ PDEs in high dimension without using grids or numerical approximations. We also present some numerical results for solving some inverse problems involving HJ PDEs using our proposed architectures., Comment: 44 pages, 9 figures

This paper provides new theoretical connections between multi-time Hamilton-Jacobi partial differential equations and variational image decomposition models in imaging sciences. We show that the minimal values of these optimization problems are governed by multi-time Hamilton-Jacobi partial differential equations. The minimizers of these optimization problems can be represented using the momentum in the corresponding Hamilton-Jacobi partial differential equation. Moreover, variational behaviors of both the minimizers and the momentum are investigated as the regularization parameters approach zero. In addition, we provide a new perspective from convex analysis to prove the uniqueness of convex solutions to Hamilton-Jacobi equations. Finally we consider image decomposition models that do not have unique minimizers and we propose a regularization approach to perform the analysis using multi-time Hamilton-Jacobi partial differential equations.

In this paper, we consider the problem of optimizing the trajectory of an Unmanned Aerial Vehicle (UAV) Base Station (BS). We consider a map characterized by a traffic intensity of users to be served. The UAV BS must travel from a given initial location at an initial time to a final position within a given duration and serves the traffic on its way. The problem consists in finding the optimal trajectory that minimizes a certain cost depending on the velocity and on the amount of served traffic. We formulate the problem using the framework of Lagrangian mechanics. When the traffic intensity is quadratic (single-phase), we derive closed-form formulas for the optimal trajectory. When the traffic intensity is bi-phase, we provide necessary conditions of optimality and propose an Alternating Optimization Algorithm that returns a trajectory satisfying these conditions. The Algorithm is initialized with a Model Predictive Control (MPC) online algorithm. Numerical results show how we improve the trajectory with respect to the MPC solution.