We study the problem of the nonparametric estimation for the density $\pi$ of the stationary distribution of a $d$-dimensional stochastic differential equation $(X_t)_{t \in [0, T]}$. From the continuous observation of the sampling path on $[0, T]$, we study the estimation of $\pi(x)$ as $T$ goes to infinity. For $d\ge2$, we characterize the minimax rate for the $\mathbf{L}^2$-risk in pointwise estimation over a class of anisotropic H\"older functions $\pi$ with regularity $\beta = (\beta_1, ... , \beta_d)$. For $d \ge 3$, our finding is that, having ordered the smoothness such that $\beta_1 \le ... \le \beta_d$, the minimax rate depends on whether $\beta_2 < \beta_3$ or $\beta_2 = \beta_3$. In the first case, this rate is $(\frac{\log T}{T})^\gamma$, and in the second case, it is $(\frac{1}{T})^\gamma$, where $\gamma$ is an explicit exponent dependent on the dimension and $\bar{\beta}_3$, the harmonic mean of smoothness over the $d$ directions after excluding $\beta_1$ and $\beta_2$, the smallest ones. We also demonstrate that kernel-based estimators achieve the optimal minimax rate. Furthermore, we propose an adaptive procedure for both $L^2$ integrated and pointwise risk. In the two-dimensional case, we show that kernel density estimators achieve the rate $\frac{\log T}{T}$, which is optimal in the minimax sense. Finally we illustrate the validity of our theoretical findings by proposing numerical results.
The problem of recovering a moment-determinate multivariate function $f$ via its moment sequence is studied. Under mild conditions on $f$, the point-wise and $L_1$-rates of convergence for the proposed constructions are established. The cases where $f$ is the indicator function of a set, and represents a discrete probability mass function are also investigated. Calculations of the approximants and simulation studies are conducted to graphically illustrate the behavior of the approximations in several simple examples. Analytical and simulated errors of proposed approximations are recorded in Tables 1-3.
We consider the coupled system of the Landau--Lifshitz--Gilbert equation and the conservation of linear momentum law to describe magnetic processes in ferromagnetic materials including magnetoelastic effects in the small-strain regime. For this nonlinear system of time-dependent partial differential equations, we present a decoupled integrator based on first-order finite elements in space and an implicit one-step method in time. We prove unconditional convergence of the sequence of discrete approximations towards a weak solution of the system as the mesh size and the time-step size go to zero. Compared to previous numerical works on this problem, for our method, we prove a discrete energy law that mimics that of the continuous problem and, passing to the limit, yields an energy inequality satisfied by weak solutions. Moreover, our method does not employ a nodal projection to impose the unit length constraint on the discrete magnetisation, so that the stability of the method does not require weakly acute meshes. Furthermore, our integrator and its analysis hold for a more general setting, including body forces and traction, as well as a more general representation of the magnetostrain. Numerical experiments underpin the theory and showcase the applicability of the scheme for the simulation of the dynamical processes involving magnetoelastic materials at submicrometer length scales.
The Riesz maps of the $L^2$ de Rham complex frequently arise as subproblems in the construction of fast preconditioners for more complicated problems. In this work we present multigrid solvers for high-order finite element discretizations of these Riesz maps with the same time and space complexity as sum-factorized operator application, i.e.~with optimal complexity in polynomial degree in the context of Krylov methods. The key idea of our approach is to build new finite elements for each space in the de Rham complex with orthogonality properties in both the $L^2$- and $H(\mathrm{d})$-inner products ($\mathrm{d} \in \{\mathrm{grad}, \mathrm{curl}, \mathrm{div}\})$ on the reference hexahedron. The resulting sparsity enables the fast solution of the patch problems arising in the Pavarino, Arnold--Falk--Winther and Hiptmair space decompositions, in the separable case. In the non-separable case, the method can be applied to an auxiliary operator that is sparse by construction. With exact Cholesky factorizations of the sparse patch problems, the application complexity is optimal but the setup costs and storage are not. We overcome this with the finer Hiptmair space decomposition and the use of incomplete Cholesky factorizations imposing the sparsity pattern arising from static condensation, which applies whether static condensation is used for the solver or not. This yields multigrid relaxations with time and space complexity that are both optimal in the polynomial degree.
This work puts forth low-complexity Riemannian subspace descent algorithms for the minimization of functions over the symmetric positive definite (SPD) manifold. Different from the existing Riemannian gradient descent variants, the proposed approach utilizes carefully chosen subspaces that allow the update to be written as a product of the Cholesky factor of the iterate and a sparse matrix. The resulting updates avoid the costly matrix operations like matrix exponentiation and dense matrix multiplication, which are generally required in almost all other Riemannian optimization algorithms on SPD manifold. We further identify a broad class of functions, arising in diverse applications, such as kernel matrix learning, covariance estimation of Gaussian distributions, maximum likelihood parameter estimation of elliptically contoured distributions, and parameter estimation in Gaussian mixture model problems, over which the Riemannian gradients can be calculated efficiently. The proposed uni-directional and multi-directional Riemannian subspace descent variants incur per-iteration complexities of $\O(n)$ and $\O(n^2)$ respectively, as compared to the $\O(n^3)$ or higher complexity incurred by all existing Riemannian gradient descent variants. The superior runtime and low per-iteration complexity of the proposed algorithms is also demonstrated via numerical tests on large-scale covariance estimation and matrix square root problems.
We extend the use of piecewise orthogonal collocation to computing periodic solutions of renewal equations, which are particularly important in modeling population dynamics. We prove convergence through a rigorous error analysis. Finally, we show some numerical experiments confirming the theoretical results, and a couple of applications in view of bifurcation analysis.
Neufeld and Wu (arXiv:2310.12545) developed a multilevel Picard (MLP) algorithm which can approximately solve general semilinear parabolic PDEs with gradient-dependent nonlinearities, allowing also for coefficient functions of the corresponding PDE to be non-constant. By introducing a particular stochastic fixed-point equation (SFPE) motivated by the Feynman-Kac representation and the Bismut-Elworthy-Li formula and identifying the first and second component of the unique fixed-point of the SFPE with the unique viscosity solution of the PDE and its gradient, they proved convergence of their algorithm. However, it remained an open question whether the proposed MLP schema in arXiv:2310.12545 does not suffer from the curse of dimensionality. In this paper, we prove that the MLP algorithm in arXiv:2310.12545 indeed can overcome the curse of dimensionality, i.e. that its computational complexity only grows polynomially in the dimension $d\in \mathbb{N}$ and the reciprocal of the accuracy $\varepsilon$, under some suitable assumptions on the nonlinear part of the corresponding PDE.
We introduce a family of identities that express general linear non-unitary evolution operators as a linear combination of unitary evolution operators, each solving a Hamiltonian simulation problem. This formulation can exponentially enhance the accuracy of the recently introduced linear combination of Hamiltonian simulation (LCHS) method [An, Liu, and Lin, Physical Review Letters, 2023]. For the first time, this approach enables quantum algorithms to solve linear differential equations with both optimal state preparation cost and near-optimal scaling in matrix queries on all parameters.
We study the stability and sensitivity of an absorbing layer for the Boltzmann equation by examining the Bhatnagar-Gross-Krook (BGK) approximation and using the perfectly matched layer (PML) technique. To ensure stability, we discard some parameters in the model and calculate the total sensitivity indices of the remaining parameters using the ANOVA expansion of multivariate functions. We conduct extensive numerical experiments to study stability and compute the total sensitivity indices, which allow us to identify the essential parameters of the model.
Quadratization of polynomial and nonpolynomial systems of ordinary differential equations is advantageous in a variety of disciplines, such as systems theory, fluid mechanics, chemical reaction modeling and mathematical analysis. A quadratization reveals new variables and structures of a model, which may be easier to analyze, simulate, control, and provides a convenient parametrization for learning. This paper presents novel theory, algorithms and software capabilities for quadratization of non-autonomous ODEs. We provide existence results, depending on the regularity of the input function, for cases when a quadratic-bilinear system can be obtained through quadratization. We further develop existence results and an algorithm that generalizes the process of quadratization for systems with arbitrary dimension that retain the nonlinear structure when the dimension grows. For such systems, we provide dimension-agnostic quadratization. An example is semi-discretized PDEs, where the nonlinear terms remain symbolically identical when the discretization size increases. As an important aspect for practical adoption of this research, we extended the capabilities of the QBee software towards both non-autonomous systems of ODEs and ODEs with arbitrary dimension. We present several examples of ODEs that were previously reported in the literature, and where our new algorithms find quadratized ODE systems with lower dimension than the previously reported lifting transformations. We further highlight an important area of quadratization: reduced-order model learning. This area can benefit significantly from working in the optimal lifting variables, where quadratic models provide a direct parametrization of the model that also avoids additional hyperreduction for the nonlinear terms. A solar wind example highlights these advantages.
While many phenomena in physics and engineering are formally high-dimensional, their long-time dynamics often live on a lower-dimensional manifold. The present work introduces an autoencoder framework that combines implicit regularization with internal linear layers and $L_2$ regularization (weight decay) to automatically estimate the underlying dimensionality of a data set, produce an orthogonal manifold coordinate system, and provide the mapping functions between the ambient space and manifold space, allowing for out-of-sample projections. We validate our framework's ability to estimate the manifold dimension for a series of datasets from dynamical systems of varying complexities and compare to other state-of-the-art estimators. We analyze the training dynamics of the network to glean insight into the mechanism of low-rank learning and find that collectively each of the implicit regularizing layers compound the low-rank representation and even self-correct during training. Analysis of gradient descent dynamics for this architecture in the linear case reveals the role of the internal linear layers in leading to faster decay of a "collective weight variable" incorporating all layers, and the role of weight decay in breaking degeneracies and thus driving convergence along directions in which no decay would occur in its absence. We show that this framework can be naturally extended for applications of state-space modeling and forecasting by generating a data-driven dynamic model of a spatiotemporally chaotic partial differential equation using only the manifold coordinates. Finally, we demonstrate that our framework is robust to hyperparameter choices.