We analyze when an arbitrary matrix pencil is equivalent to a dissipative Hamiltonian pencil and show that this heavily restricts the spectral properties. In order to relax the spectral properties, we introduce matrix pencils with coefficients that have positive semidefinite Hermitian parts. We will make a detailed analysis of their spectral properties and their numerical range. In particular, we relate the Kronecker structure of these pencils to that of an underlying skew-Hermitian pencil and discuss their regularity, index, numerical range, and location of eigenvalues. Further, we study matrix polynomials with positive semidefinite Hermitian coefficients and use linearizations with positive semidefinite Hermitian parts to derive sufficient conditions for a spectrum in the left half plane and derive bounds on the index.
In this paper, we consider several efficient data structures for the problem of sampling from a dynamically changing discrete probability distribution, where some prior information is known on the distribution of the rates, in particular the maximum and minimum rate, and where the number of possible outcomes N is large. We consider three basic data structures, the Acceptance-Rejection method, the Complete Binary Tree and the Alias method. These can be used as building blocks in a multi-level data structure, where at each of the levels, one of the basic data structures can be used, with the top level selecting a group of events, and the bottom level selecting an element from a group. Depending on assumptions on the distribution of the rates of outcomes, different combinations of the basic structures can be used. We prove that for particular data structures the expected time of sampling and update is constant when the rate distribution follows certain conditions. We show that for any distribution, combining a tree structure with the Acceptance-Rejection method, we have an expected time of sampling and update of $O\left(\log\log{r_{max}}/{r_{min}}\right)$ is possible, where $r_{max}$ is the maximum rate and $r_{min}$ the minimum rate. We also discuss an implementation of a Two Levels Acceptance-Rejection data structure, that allows expected constant time for sampling, and amortized constant time for updates, assuming that $r_{max}$ and $r_{min}$ are known and the number of events is sufficiently large. We also present an experimental verification, highlighting the limits given by the constraints of a real-life setting.
This paper develops a robust dynamic mode decomposition (RDMD) method endowed with statistical and numerical robustness. Statistical robustness ensures estimation efficiency at the Gaussian and non-Gaussian probability distributions, including heavy-tailed distributions. The proposed RDMD is statistically robust because the outliers in the data set are flagged via projection statistics and suppressed using a Schweppe-type Huber generalized maximum-likelihood estimator that minimizes a convex Huber cost function. The latter is solved using the iteratively reweighted least-squares algorithm that is known to exhibit a better convergence property and numerical stability than the Newton algorithms. Several numerical simulations using canonical models of dynamical systems demonstrate the excellent performance of the proposed RDMD method. The results reveal that it outperforms several other methods proposed in the literature.
We present a new enriched Galerkin (EG) scheme for the Stokes equations based on piecewise linear elements for the velocity unknowns and piecewise constant elements for the pressure. The proposed EG method augments the conforming piecewise linear space for velocity by adding an additional degree of freedom which corresponds to one discontinuous linear basis function per element. Thus, the total number of degrees of freedom is significantly reduced in comparison with standard conforming, non-conforming, and discontinuous Galerkin schemes for the Stokes equation. We show the well-posedness of the new EG approach and prove that the scheme converges optimally. For the solution of the resulting large-scale indefinite linear systems we propose robust block preconditioners, yielding scalable results independent of the discretization and physical parameters. Numerical results confirm the convergence rates of the discretization and also the robustness of the linear solvers for a variety of test problems.
This paper deals with the problem of graph matching or network alignment for Erd\H{o}s--R\'enyi graphs, which can be viewed as a noisy average-case version of the graph isomorphism problem. Let $G$ and $G'$ be $G(n, p)$ Erd\H{o}s--R\'enyi graphs marginally, identified with their adjacency matrices. Assume that $G$ and $G'$ are correlated such that $\mathbb{E}[G_{ij} G'_{ij}] = p(1-\alpha)$. For a permutation $\pi$ representing a latent matching between the vertices of $G$ and $G'$, denote by $G^\pi$ the graph obtained from permuting the vertices of $G$ by $\pi$. Observing $G^\pi$ and $G'$, we aim to recover the matching $\pi$. In this work, we show that for every $\varepsilon \in (0,1]$, there is $n_0>0$ depending on $\varepsilon$ and absolute constants $\alpha_0, R > 0$ with the following property. Let $n \ge n_0$, $(1+\varepsilon) \log n \le np \le n^{\frac{1}{R \log \log n}}$, and $0 < \alpha < \min(\alpha_0,\varepsilon/4)$. There is a polynomial-time algorithm $F$ such that $\mathbb{P}\{F(G^\pi,G')=\pi\}=1-o(1)$. This is the first polynomial-time algorithm that recovers the exact matching between vertices of correlated Erd\H{o}s--R\'enyi graphs with constant correlation with high probability. The algorithm is based on comparison of partition trees associated with the graph vertices.
In this article, we develop and analyse a new spectral method to solve the semi-classical Schr\"odinger equation based on the Gaussian wave-packet transform (GWPT) and Hagedorn's semi-classical wave-packets (HWP). The GWPT equivalently recasts the highly oscillatory wave equation as a much less oscillatory one (the $w$ equation) coupled with a set of ordinary differential equations governing the dynamics of the so-called GWPT parameters. The Hamiltonian of the $ w $ equation consists of a quadratic part and a small non-quadratic perturbation, which is of order $ \mathcal{O}(\sqrt{\varepsilon }) $, where $ \varepsilon\ll 1 $ is the rescaled Planck's constant. By expanding the solution of the $ w $ equation as a superposition of Hagedorn's wave-packets, we construct a spectral method while the $ \mathcal{O}(\sqrt{\varepsilon}) $ perturbation part is treated by the Galerkin approximation. This numerical implementation of the GWPT avoids imposing artificial boundary conditions and facilitates rigorous numerical analysis. For arbitrary dimensional cases, we establish how the error of solving the semi-classical Schr\"odinger equation with the GWPT is determined by the errors of solving the $ w $ equation and the GWPT parameters. We prove that this scheme has the spectral convergence with respect to the number of Hagedorn's wave-packets in one dimension. Extensive numerical tests are provided to demonstrate the properties of the proposed method.
Time-space fractional Bloch-Torrey equations (TSFBTEs) are developed by some researchers to investigate the relationship between diffusion and fractional-order dynamics. In this paper, we first propose a second-order implicit difference scheme for TSFBTEs by employing the recently proposed L2-type formula [A. A. Alikhanov, C. Huang, Appl. Math. Comput. (2021) 126545]. Then, we prove the stability and the convergence of the proposed scheme. Based on such a numerical scheme, an L2-type all-at-once system is derived. In order to solve this system in a parallel-in-time pattern, a bilateral preconditioning technique is designed to accelerate the convergence of Krylov subspace solvers according to the special structure of the coefficient matrix of the system. We theoretically show that the condition number of the preconditioned matrix is uniformly bounded by a constant for the time fractional order $\alpha \in (0,0.3624)$. Numerical results are reported to show the efficiency of our method.
We propose and analyze a stochastic Newton algorithm for homogeneous distributed stochastic convex optimization, where each machine can calculate stochastic gradients of the same population objective, as well as stochastic Hessian-vector products (products of an independent unbiased estimator of the Hessian of the population objective with arbitrary vectors), with many such stochastic computations performed between rounds of communication. We show that our method can reduce the number, and frequency, of required communication rounds compared to existing methods without hurting performance, by proving convergence guarantees for quasi-self-concordant objectives (e.g., logistic regression), alongside empirical evidence.
We study the MARINA method of Gorbunov et al (2021) -- the current state-of-the-art distributed non-convex optimization method in terms of theoretical communication complexity. Theoretical superiority of this method can be largely attributed to two sources: the use of a carefully engineered biased stochastic gradient estimator, which leads to a reduction in the number of communication rounds, and the reliance on {\em independent} stochastic communication compression operators, which leads to a reduction in the number of transmitted bits within each communication round. In this paper we i) extend the theory of MARINA to support a much wider class of potentially {\em correlated} compressors, extending the reach of the method beyond the classical independent compressors setting, ii) show that a new quantity, for which we coin the name {\em Hessian variance}, allows us to significantly refine the original analysis of MARINA without any additional assumptions, and iii) identify a special class of correlated compressors based on the idea of {\em random permutations}, for which we coin the term Perm$K$, the use of which leads to $O(\sqrt{n})$ (resp. $O(1 + d/\sqrt{n})$) improvement in the theoretical communication complexity of MARINA in the low Hessian variance regime when $d\geq n$ (resp. $d \leq n$), where $n$ is the number of workers and $d$ is the number of parameters describing the model we are learning. We corroborate our theoretical results with carefully engineered synthetic experiments with minimizing the average of nonconvex quadratics, and on autoencoder training with the MNIST dataset.
In this paper, we consider the asymptotical regularization with convex constraints for nonlinear ill-posed problems. The method allows to use non-smooth penalty terms, including the L1-like and the total variation-like penalty functionals, which are significant in reconstructing special features of solutions such as sparsity and piecewise constancy. Under certain conditions we give convergence properties of the methods. Moreover, we propose Runge-Kutta type methods to discrete the initial value problems to construct new type iterative regularization methods.
In this paper a two-sided, parallel Kogbetliantz-type algorithm for the hyperbolic singular value decomposition (HSVD) of real and complex square matrices is developed, with a single assumption that the input matrix, of order $n$, admits such a decomposition into the product of a unitary, a non-negative diagonal, and a $J$-unitary matrix, where $J$ is a given diagonal matrix of positive and negative signs. When $J=\pm I$, the proposed algorithm computes the ordinary SVD. The paper's most important contribution -- a derivation of formulas for the HSVD of $2\times 2$ matrices -- is presented first, followed by the details of their implementation in floating-point arithmetic. Next, the effects of the hyperbolic transformations on the columns of the iteration matrix are discussed. These effects then guide a redesign of the dynamic pivot ordering, being already a well-established pivot strategy for the ordinary Kogbetliantz algorithm, for the general, $n\times n$ HSVD. A heuristic but sound convergence criterion is then proposed, which contributes to high accuracy demonstrated in the numerical testing results. Such a $J$-Kogbetliantz algorithm as presented here is intrinsically slow, but is nevertheless usable for matrices of small orders.