This work studies the invariant measure of numerical approximation to a class of superlinear stochastic differential equations (SDEs) with periodic coefficients by the truncated Euler–Maruyama (EM) method. The existing studies on numerical invariant measures mostly focus on autonomous SDEs. To the best of our knowledge, this paper is the first devoted to the case of non-autonomous SDEs which allow superlinear growth. Technical challenges including superlinearity, periodicity and time-inhomogeneity of the system make this a non-trivial work. Not only the existence and uniqueness of a numerical invariant measure but also the convergence of the numerical invariant measure to the underlying one are deduced in the Wasserstein metric. Consequently, a case study is carried out to demonstrate the main results.
In this paper, we consider the low rank matrix feasibility problem arising in quantum marginal computing. Based on the alternating projection method and Newton’s method, we propose the Alternating-Newton algorithm to solve this problem. The convergence analysis of the Alternating-Newton algorithm is given and its quadratic convergence rate is also derived. Some numerical experiments demonstrate the feasibility and efficiency of the proposed algorithm.
In this paper, we propose three unisolvent weighted quadratic enrichments to extend classical linear histopolation on 3D tetrahedral meshes. The first combines face and interior weighted moments, the second relies exclusively on volumetric quadratic moments, and the third employs edge supported probabilistic moments. In all three cases, the degrees of freedom are defined by weighted integral functionals associated with suitable probability densities and quadratic trial spaces. We establish the unisolvence of the enriched schemes and derive explicit conditions under which this property holds. Representative density families, including two-parameter symmetric Dirichlet laws and volumetric families defined by convex combinations of densities, are examined in detail, and a general procedure for constructing the associated quadratic basis functions is presented. For the considered parametric density families, the corresponding parameters are selected via a grid search procedure. Numerical experiments show improved accuracy compared to the classical linear histopolation scheme.
The quaternion CUR (qCUR) matrix decomposition constructs low-rank approximations of quaternion matrices by selecting a few key rows and columns, yet existing qCUR approximation methods typically fail to provide tight error bounds. This work focuses on constructing qCUR via leverage score sampling, and derives novel gap error bounds for the rank-k qCUR decomposition: we leverage real representations of quaternion random matrices (which connect these matrices to probability theory) to synthesize the combined impacts of column and row selection strategies, and the derived bound intrinsically ties the decomposition’s approximation accuracy to the decay rate of the matrix’s singular values. Specifically, when sampling is performed on the p ( p≥ k ) exact dominant singular vectors, the rank-k approximation exhibits a relative error (relative to the best rank-k approximation error) of O( εσ _p+1^2/σ _k+1^2) , where σ _j denotes the j-th largest singular value and the number of selected columns and rows satisfies c, r = O( plog p/ε) . The approximation error is further analyzed for cases where the leverage score is computed based on approximate dominant singular vectors from randomized quaternion singular value decomposition. Empirical results on various classes of synthetic data and real-world color image completion tasks show that our algorithm outperforms the qCUR variants with squared-length sampling or the discrete empirical interpolation method.
We investigate numerical aspects of Riemannian interpolation on the Grassmann manifold. Instead of relying on the Riemannian normal coordinates, i.e. the Riemannian exponential and logarithm maps, we approach the interpolation problem with an alternative set of local coordinates and corresponding parameterizations. We show that these coordinates define a second-order retraction. Numerical evaluation does not formally require matrix decompositions. This is an advantage over Riemannian normal coordinates and many other retractions on the Grassmann manifold, especially when derivative data are to be treated. To estimate the interpolation error, we examine the conditioning of the coordinate mappings and state explicit bounds. It turns out that the parameterizations are well-conditioned, but the coordinate charts are generally not. As a remedy, we introduce canonically centered coordinates based on a homogeneous transition to a canonical Stiefel representative of a Grassmann data point chosen to act as the coordinate center. We show that the order of magnitude of the asymptotic interpolation error on Gr(n,p) is the same as in the Euclidean space. Numerical experiments illustrate the findings. The first is academic, where we interpolate a parametric orthogonal projector QQ^T . The Q–factor stems from a parametric compact QR–decomposition. Moreover, we assess the computation time, present an example where closed Riemannian normal coordinates struggle, and conduct an experiment in the context of parametric model reduction of dynamical systems, where we interpolate reduced subspaces that are obtained by proper orthogonal decomposition.
This paper presents stability and accuracy analysis of a high-order explicit time stepping scheme introduced by [4, Section 2.2], which exhibits superior stability compared to classical Adams-Bashforth. A conjecture that is supported by several numerical phenomena in [3, Figure 2.5], the method appears to remain stable when the accuracy approaches infinity, although it is not yet proven. We have disproven this conjecture from the perspective of harmonic analysis in this work. Notwithstanding the aforementioned, this method displays considerably enhanced stability in comparison to conventional explicit schemes. Furthermore, we present a criterion for ascertaining the maximum permissible accuracy for a given specific parabolic stability radius. Conversely, the original method will lose one order associated with the expected accuracy, which can be explored theoretically. Consequently, a unified analysis strategy for the L^2 -stability will be presented for extensional PDEs under the CFL condition. Finally, a selection of representative numerical examples will be shown in order to substantiate the theoretical analysis.
In this paper, we study the convergence behavior of the diffuse domain method (DDM) for solving a class of second-order parabolic partial differential equations with Neumann boundary condition posed on general irregular domains. The DDM employs a phase-field function to extend the original parabolic problem to a similar but slightly modified problem defined over a larger rectangular domain that contains the target physical domain. Based on the weighted Sobolev spaces, we rigorously establish the convergence of the diffuse domain solution to the original solution as the interface thickness parameter goes to zero, together with the corresponding optimal error estimates under the weighted L^2 and H^1 norms. Numerical experiments are also presented to validate the theoretical results.
In this paper, we analyze stability properties of the two-derivative strong stability preserving schemes presented in [Gottlieb et al., SIAM Journal on Numerical Analysis 60, 2022]. Stability analysis shows that the diagonally implicit two-derivative two-stage third-order strong stability preserving scheme can never be A-stable. We provide a detailed investigation of the third-order schemes and discuss stabilizing strategies. The stabilizing techniques are applicable to tune any general implicit two-derivative scheme. We implement the two-derivative strong stability preserving schemes for partial differential equations with a discontinuous Galerkin spectral element spatial discretization. We use Newton’s method for non-linear stage equations and the generalized minimal residual method with a matrix-free approach for solving linear algebraic equations under suitable preconditioning. The method is applied for compressible Euler and Navier-Stokes equations with orders up to four. Numerical results show that the second and fourth-order strong stability preserving schemes attain their desired order of convergence for relatively large timesteps. In contrast, third-order schemes require smaller timesteps to exhibit convergence. Nevertheless, the improved adaptive third-order scheme yields stable solutions.
In recent years, the strong convergence analysis of positivity-preserving numerical methods and the weak convergence analysis of numerical methods for stochastic differential equations in general have garnered widespread attention. However, research focusing on the weak convergence analysis of positivity-preserving methods remains limited. In particular, the weak convergence analysis of these methods in the multi-dimensional setting remains unexplored, which constitutes the core objective of the present work. In this paper, we first prove the existence and uniqueness of positive strong solution of the underlying equations under certain conditions. Then, we investigate the weak convergence of the positivity-preserving truncated Euler–Maruyama method and demonstrate that, under additional conditions, its convergence order can be made arbitrarily close to 1. Lastly, numerical experiments are conducted to validate the theoretical findings.
We consider four product integration rules, two for the Chebyshev weight of the first-kind based on the Chebyshev abscissae of the third or fourth-kind, and another two for the Chebyshev weight of the second-kind based again on the Chebyshev abscissae of the third or fourth-kind. The new rules have positive weights given by explicit formulae, while the rules for the Chebyshev weight of the second-kind have the best possible degree of exactness for an interpolatory formula not of Gauss type. On certain spaces of analytic functions the error term of these rules is a continuous linear functional. By means of a new approach, we compute explicitly the norm of the error functional, which leads to efficient error bounds.
Equivalent preservation of asymptotic mean square stability and instability of balanced midpoint Milstein methods (BMMMs) applied to stochastic differential equations (SDEs) driven by standard Wiener processes is shown whenever the underlying SDE has an asymptotically mean square stable equilibrium or not in ℂ^1 , respectively (called mean square E-stability by this article). These are certain numerical methods built up by the class of implicit Milstein methods combined with midpoint drift-implicitness and additional balanced terms. The paper verifies that it is indeed possible to construct such higher order numerical methods for SDEs, which are asymptotically mean square stable for all possible step sizes h>0 if and only if the standard test class of bi-linear SDEs (the stochastic 1D test equation in ℂ^1 in the sense of Dahlquist) with multiplicative noise has an asymptotically mean square stable trivial solution. This investigation goes far beyond the common requirement of mean square A-stability of stochastic-numerical methods. Previously, it has been shown that this is the case for the drift-implicit midpoint method, which represents a lower mean square order stochastic Theta method for SDEs with parameter θ =0.5 .
The paper deals with bounds for Krylov methods which are insensitive in low rank perturbations. In finite dimensional cases resolvents are meromorphic in the whole plane and robust bounds have been constructed using special growth functions created for operator valued meromorphic functions. In this paper such bounds are derived without use of those special tools. In particular, convergence in generic hermitean problems and highly non-normal problems are effectively analysed with the same technique based on representing the resolvent using spectral polynomials and thus for example the conditioning of eigenvector bases does not show up at all.
Principal Component Analysis (PCA) is a foundational technique in machine learning for reducing the dimensionality of high-dimensional datasets. However, PCA can lead to biased representations that disadvantage certain subgroups within the data. To address this issue, a Fair PCA (FPCA) model was introduced to equalize the reconstruction loss between subgroups, but the existing semidefinite relaxation (SDR) based approach is computationally expensive even for a suboptimal solution. Although several alternative FPCA variants have been developed to improve efficiency, they often shift attention away from equalizing the reconstruction loss – the central goal of FPCA. In this paper, we identify a hidden convexity in FPCA and introduce a new algorithm that solves the resulting convex optimization via an eigenvalue optimization. Our approach achieves the desired fairness in reconstruction loss without sacrificing performance. Experiments on real-world datasets show that the proposed FPCA algorithm is approximately 8× faster than the SDR-based algorithm while being at most 85
We present a hybrid a–priori/a–posteriori goal–oriented error estimator for a combination of dynamic iteration based solution of linear ordinary differential equations discretized by finite elements. Our novel error estimator combines estimates from classical dynamic iteration methods, usually used to enable splitting–based distributed simulation, and from the dual weighted residual method to be able to evaluate and balance both, the dynamic iteration error and the discretization error in desired quantities of interest. The obtained error estimators are used to conduct refinements of the computational mesh and as a stopping criterion for the dynamic iteration. In particular, we allow for an adaptive and flexible discretization of the time domain, where variables can be discretized differently to match both goal and solution requirements, e.g. in view of multiple time scales. We endow the scheme with efficient solvers from numerical linear algebra to ensure its applicability to complex problems. Numerical experiments compare the adaptive approach to a uniform refinement.
We study Kaczmarz type methods to solve consistent linear matrix equations. We first present a block Kaczmarz (BK) method that employs a deterministic cyclic row selection strategy. Assuming that the associated coefficient matrix has full column or row rank, we derive matrix formulas for a cycle of this BK method. Moreover, we propose a greedy randomized block Kaczmarz (GRBK) method and further extend it to a relaxed variant (RGRBK) and a deterministic counterpart (MWRBK). We establish the convergence properties of the proposed methods. Numerical tests verify the theoretical findings, and we apply the proposed methods to color image restoration problems.
A novel class of bivariate α -fractal functions is constructed on rectangular grids in this paper, with sufficient conditions established for these α -fractal functions to become α -fractal interpolation functions ( α -FIFs). The box-counting dimension of the resulting α -FIFs is subsequently estimated. Three numerical examples are provided to illustrate the effectiveness and practical applicability of the proposed construction method. Finally, it is demonstrated that the Riemann-Liouville fractional integrals of the constructed bivariate α -fractal functions remain α -fractal functions.
Sketching techniques have gained popularity in numerical linear algebra to accelerate the solution of least squares problems. The so-called epsilon-subspace embedding property of a sketching matrix S has been largely used to characterize the problem residual norm, since the procedure is no longer optimal in terms of the (classical) Frobenius or Euclidean norm. By building on available results on the SVD of the sketched matrix SA derived by Gilbert, Park, and Wakin (Proc. of SPARS-2013), a novel decomposition of A, the S inverted perpendicular S-SVD, is proposed, which holds with high probability, and in which the left singular vectors are orthonormal with respect to a (semi-)norm defined by the sketching matrix S. The new decomposition is less expensive to compute than the standard SVD, while preserving the singular values with probabilistic confidence. The S inverted perpendicular S-SVD appears to be the right tool to analyze the quality of several sketching-based techniques in the literature, for which examples are reported. For instance, it is possible to simply bound the distance from (standard) orthogonality of sketching-based orthogonal matrices in state-of-the-art randomized algorithms for QR factorizations. As an application, the classical problem of the nearest orthogonal matrix is generalized to the new S inverted perpendicular S-orthogonality, and the S inverted perpendicular S-SVD is used to solve it. Probabilistic bounds on the quality of the solution are also derived.
In this paper, a class of two-scale mixed finite element methods, including the two-scale mixed finite element method and the two-scale postprocessed mixed finite element method, is proposed for solving the Stokes problem. Both methods are grounded in the Bernardi-Raugel element. The essence of the two-scale mixed finite element solution lies in combining the mixed finite element solution from a coarse grid with those from several univariate fine grids. The underlying principle is that the low-frequency components of the mixed finite element solution can be adequately captured on the coarse grid, while the high-frequency components are more efficiently handled by the univariate fine grids. Both theoretical and numerical analyses demonstrate that the two-scale mixed finite element solution achieves the same order of accuracy as the traditional mixed finite element solution, yet at substantially lower computational costs. Additionally, we incorporate a postprocessing technique to develop and analyze the two-scale postprocessed mixed finite element method. This method is even more computationally efficient surpassing both the postprocessed mixed finite element method and the two-scale mixed finite element method.
The Kaczmarz method is successfully used for solving discretizations of linear inverse problems, especially in computed tomography where it is known as ART. Practitioners often observe and appreciate its fast convergence in the first few iterations, leading to the same favorable semi-convergence that we observe for simultaneous iterative reconstruction methods. While the latter methods have symmetric and positive definite iteration operators that facilitate their analysis, the operator in Kaczmarz's method is nonsymmetric and it has been an open question so far to understand this fast initial convergence. We perform a spectral analysis of Kaczmarz's method that gives new insight into its (often fast) initial behavior. We also carry out a statistical analysis of how the data noise enters the iteration vectors, which sheds new light on the semi-convergence. Our results are illustrated with several numerical examples.
Matrix joint block-diagonalization (jbd) frequently arises from diverse applications such as independent component analysis, blind source separation, and common principal component analysis (CPCA), among others. Particularly, CPCA aims at joint diagonalization, i.e., each block size being 1-by-1. This paper is concerned with principal joint block-diagonalization (pjbd), which aims to achieve two goals: 1) partial joint block-diagonalization, and 2) identification of dominant common block-diagonal parts for all involved matrices. This is in contrast to most existing methods, especially the popular ones based on Givens rotation, which focus on full joint diagonalization and quickly become impractical for matrices of even moderate size (300-by-300 or larger). An NPDo approach, directly aiming at dominant common block-diagonal parts, is proposed and it is built on the nonlinear polar decomposition with orthonormal polar factor dependency that characterizes the solutions of the optimization problem designed to achieve pjbd, and it is shown the associated SCF iteration is globally convergent to a stationary point while the objective function increases monotonically during the iterative process. Numerical experiments, including practical applications to multi-view subspace clustering and independent subspace analysis, are presented to illustrate the effectiveness of the NPDo approach and its superiority to Givens rotation-based methods.
