We investigate the performance of two approximation algorithms for the Hafnian of a nonnegative square matrix, namely the Barvinok and Godsil-Gutman estimators. We observe that, while there are examples of matrices for which these algorithms fail to provide a good approximation, the algorithms perform surprisingly well for adjacency matrices of random graphs. In most cases, the Godsil-Gutman estimator provides a far superior accuracy. For dense graphs, however, both estimators demonstrate a slow growth of the variance. For complete graphs, we show analytically that the relative variance $\sigma / \mu$ grows as a square root of the size of the graph. Finally, we simulate a Gaussian Boson Sampling experiment using the Godsil-Gutman estimator and show that the technique used can successfully reproduce low-order correlation functions.
Differential equation models are crucial to scientific processes. The values of model parameters are important for analyzing the behaviour of solutions. A parameter is called globally identifiable if its value can be uniquely determined from the input and output functions. To determine if a parameter estimation problem is well-posed for a given model, one must check if the model parameters are globally identifiable. This problem has been intensively studied for ordinary differential equation models, with theory and several efficient algorithms and software packages developed. A comprehensive theory of algebraic identifiability for PDEs has hitherto not been developed due to the complexity of initial and boundary conditions. Here, we provide theory and algorithms, based on differential algebra, for testing identifiability of polynomial PDE models. We showcase this approach on PDE models arising in the sciences.
We introduce free probability analogues of the stochastic theta methods for free stochastic differential equations, which generalize the free Euler-Maruyama method introduced by Schl\"{u}chtermann and Wibmer [27]. Under some mild conditions, we prove the strong convergence and exponential stability in mean square of the numerical solution. The free stochastic theta method with $\theta=1$ can inherit the exponential stability of original equations for any given step size. Our method can offer better stability and efficiency than the free Euler-Maruyama method. Moreover, numerical results are reported to confirm these theoretical findings.
Nonlinear conservation laws such as the system of ideal magnetohydrodynamics (MHD) equations may develop singularities over time. In these situations, viscous regularization is a common approach to regain regularity of the solution. In this paper, we present a new viscous flux to regularize the MHD equations which holds many attractive properties. In particular, we prove that the proposed viscous flux preserves positivity of density and internal energy, satisfies the minimum entropy principle, is consistent with all generalized entropies, and is Galilean and rotationally invariant. We also provide a variation of the viscous flux that conserves angular momentum. To make the analysis more useful for numerical schemes, the divergence of the magnetic field is not assumed to be zero. Using continuous finite elements, we show several numerical experiments including contact waves and magnetic reconnection.
This paper is devoted to the theoretical and numerical investigation of the initial boundary value problem for a system of equations used for the description of waves in coastal areas, namely, the Boussinesq-Abbott system in the presence of topography. We propose a procedure that allows one to handle very general linear or nonlinear boundary conditions. It consists in reducing the problem to a system of conservation laws with nonlocal fluxes and coupled to an ODE. This reformulation is used to propose two hybrid finite volumes/finite differences schemes of first and second order respectively. The possibility to use many kinds of boundary conditions is used to investigate numerically the asymptotic stability of the boundary conditions, which is an issue of practical relevance in coastal oceanography since asymptotically stable boundary conditions would allow one to reconstruct a wave field from the knowledge of the boundary data only, even if the initial data is not known.
We study the approximation by a Voronoi finite-volume scheme of the Gross-Pitaevskii equation with time-dependent potential in two and three dimensions. We perform an explicit splitting scheme for the time integration alongside a two-point flux approximation scheme in space. We rigorously analyze the error bounds relying on discrete uniform Sobolev inequalities. We also prove the convergence of the pseudo-vorticity of the wave function. We finally perform some numerical simulations to illustrate our theoretical results.
We consider the estimation of the cumulative hazard function, and equivalently the distribution function, with censored data under a setup that preserves the privacy of the survival database. This is done through a $\alpha$-locally differentially private mechanism for the failure indicators and by proposing a non-parametric kernel estimator for the cumulative hazard function that remains consistent under the privatization. Under mild conditions, we also prove lowers bounds for the minimax rates of convergence and show that estimator is minimax optimal under a well-chosen bandwidth.
The joint bidiagonalization (JBD) process iteratively reduces a matrix pair $\{A,L\}$ to two bidiagonal forms simultaneously, which can be used for computing a partial generalized singular value decomposition (GSVD) of $\{A,L\}$. The process has a nested inner-outer iteration structure, where the inner iteration usually can not be computed exactly. In this paper, we study the inaccurately computed inner iterations of JBD by first investigating influence of computational error of the inner iteration on the outer iteration, and then proposing a reorthogonalized JBD (rJBD) process to keep orthogonality of a part of Lanczos vectors. An error analysis of the rJBD is carried out to build up connections with Lanczos bidiagonalizations. The results are then used to investigate convergence and accuracy of the rJBD based GSVD computation. It is shown that the accuracy of computed GSVD components depend on the computing accuracy of inner iterations and condition number of $(A^T,L^T)^T$ while the convergence rate is not affected very much. For practical JBD based GSVD computations, our results can provide a guideline for choosing a proper computing accuracy of inner iterations in order to obtain approximate GSVD components with a desired accuracy. Numerical experiments are made to confirm our theoretical results.
An efficient approximate version of implicit Taylor methods for initial-value problems of systems of ordinary differential equations (ODEs) is introduced. The approach, based on an approximate formulation of Taylor methods, produces a method that requires less evaluations of the function that defines the ODE and its derivatives than the usual version. On the other hand, an efficient numerical solution of the equation that arises from the discretization by means of Newton's method is introduced for an implicit scheme of any order. Numerical experiments illustrate that the resulting algorithm is simpler to implement and has better performance than its exact counterpart.
We develop a sparse spectral method for a class of fractional differential equations, posed on $\mathbb{R}$, in one dimension. These equations can include sqrt-Laplacian, Hilbert, derivative and identity terms. The numerical method utilizes a basis consisting of weighted Chebyshev polynomials of the second kind in conjunction with their Hilbert transforms. The former functions are supported on $[-1,1]$ whereas the latter have global support. The global approximation space can contain different affine transformations of the basis, mapping $[-1,1]$ to other intervals. Remarkably, not only are the induced linear systems sparse, but the operator decouples across the different affine transformations. Hence, the solve reduces to solving $K$ independent sparse linear systems of size $\mathcal{O}(n)\times \mathcal{O}(n)$, with $\mathcal{O}(n)$ nonzero entries, where $K$ is the number of different intervals and $n$ is the highest polynomial degree contained in the sum space. This results in an $\mathcal{O}(n)$ complexity solve. Applications to fractional heat and wave equations are considered.
We provide a Lyapunov convergence analysis for time-inhomogeneous variable coefficient stochastic differential equations (SDEs). Three typical examples include overdamped, irreversible drift, and underdamped Langevin dynamics. We first formula the probability transition equation of Langevin dynamics as a modified gradient flow of the Kullback-Leibler divergence in the probability space with respect to time-dependent optimal transport metrics. This formulation contains both gradient and non-gradient directions depending on a class of time-dependent target distribution. We then select a time-dependent relative Fisher information functional as a Lyapunov functional. We develop a time-dependent Hessian matrix condition, which guarantees the convergence of the probability density function of the SDE. We verify the proposed conditions for several time-inhomogeneous Langevin dynamics. For the overdamped Langevin dynamics, we prove the $O(t^{-1/2})$ convergence in $L^1$ distance for the simulated annealing dynamics with a strongly convex potential function. For the irreversible drift Langevin dynamics, we prove an improved convergence towards the target distribution in an asymptotic regime. We also verify the convergence condition for the underdamped Langevin dynamics. Numerical examples demonstrate the convergence results for the time-dependent Langevin dynamics.