Numerical computation of deflating subspaces of skew-Hamiltonian/Hamiltonian pencils

Working title was “Numerical Computation of Deflating Subspaces of Embedded Hamiltonian and Symplectic Pencils"We discuss the numerical solution of structured generalized eigenvalue problems that arise from linear- quadratic optimal control problems, H infinity optimization, multibody systems, and many other areas of applied mathematics, physics, and chemistry. The classical approach for these problems requires computing invariant and deflating subspaces of matrices and matrix pencils with Hamiltonian and/ or skew- Hamiltonian structure. We extend the recently developed methods for Hamiltonian matrices to the general case of skew- Hamiltonian/ Hamiltonian pencils. The algorithms circumvent problems with skew- Hamiltonian/ Hamiltonian matrix pencils that lack structured Schur forms by embedding them into matrix pencils that always admit a structured Schur form. The rounding error analysis of the resulting algorithms is favorable. For the embedded matrix pencils, the algorithms use structure- preserving unitary matrix computations and are strongly backwards stable, i. e., they compute the exact structured Schur form of a nearby matrix pencil with the same structure.All authors were partially supported by Deutsche Forschungsgemeinschaft, Research Grant Me 790/7-2 and Sonderforschungsbereich 393, “Numerische Simulation auf massiv parallelen Rechnern”

Solving polynomial eigenvalue problems by means of the Ehrlich-Aberth method

Given the $n\times n$ matrix polynomial $P(x)=\sum_{i=0}^kP_i x^i$, we consider the associated polynomial eigenvalue problem. This problem, viewed in terms of computing the roots of the scalar polynomial $\det P(x)$, is treated in polynomial form rather than in matrix form by means of the Ehrlich-Aberth iteration. The main computational issues are discussed, namely, the choice of the starting approximations needed to start the Ehrlich-Aberth iteration, the computation of the Newton correction, the halting criterion, and the treatment of eigenvalues at infinity. We arrive at an effective implementation which provides more accurate approximations to the eigenvalues with respect to the methods based on the QZ algorithm. The case of polynomials having special structures, like palindromic, Hamiltonian, symplectic, etc., where the eigenvalues have special symmetries in the complex plane, is considered. A general way to adapt the Ehrlich-Aberth iteration to structured matrix polynomial is introduced. Numerical experiments which confirm the effectiveness of this approach are reported.Comment: Submitted to Linear Algebra App

Perturbation, extraction and refinement of invariant pairs for matrix polynomials

Generalizing the notion of an eigenvector, invariant subspaces are frequently used in the context of linear eigenvalue problems, leading to conceptually elegant and numerically stable formulations in applications that require the computation of several eigenvalues and/or eigenvectors. Similar benefits can be expected for polynomial eigenvalue problems, for which the concept of an invariant subspace needs to be replaced by the concept of an invariant pair. Little has been known so far about numerical aspects of such invariant pairs. The aim of this paper is to fill this gap. The behavior of invariant pairs under perturbations of the matrix polynomial is studied and a first-order perturbation expansion is given. From a computational point of view, we investigate how to best extract invariant pairs from a linearization of the matrix polynomial. Moreover, we describe efficient refinement procedures directly based on the polynomial formulation. Numerical experiments with matrix polynomials from a number of applications demonstrate the effectiveness of our extraction and refinement procedures

Block Circulant and Toeplitz Structures in the Linearized Hartree–Fock Equation on Finite Lattices: Tensor Approach

This paper introduces and analyses the new grid-based tensor approach to approximate solution of the elliptic eigenvalue problem for the 3D lattice-structured systems. We consider the linearized Hartree-Fock equation over a spatial $L_1\times L_2\times L_3$ lattice for both periodic and non-periodic problem setting, discretized in the localized Gaussian-type orbitals basis. In the periodic case, the Galerkin system matrix obeys a three-level block-circulant structure that allows the FFT-based diagonalization, while for the finite extended systems in a box (Dirichlet boundary conditions) we arrive at the perturbed block-Toeplitz representation providing fast matrix-vector multiplication and low storage size. The proposed grid-based tensor techniques manifest the twofold benefits: (a) the entries of the Fock matrix are computed by 1D operations using low-rank tensors represented on a 3D grid, (b) in the periodic case the low-rank tensor structure in the diagonal blocks of the Fock matrix in the Fourier space reduces the conventional 3D FFT to the product of 1D FFTs. Lattice type systems in a box with Dirichlet boundary conditions are treated numerically by our previous tensor solver for single molecules, which makes possible calculations on rather large $L_1\times L_2\times L_3$ lattices due to reduced numerical cost for 3D problems. The numerical simulations for both box-type and periodic $L\times 1\times 1$ lattice chain in a 3D rectangular "tube" with $L$ up to several hundred confirm the theoretical complexity bounds for the block-structured eigenvalue solvers in the limit of large $L$.Comment: 30 pages, 12 figures. arXiv admin note: substantial text overlap with arXiv:1408.383

Designing structured tight frames via an alternating projection method

Tight frames, also known as general Welch-bound- equality sequences, generalize orthonormal systems. Numerous applications - including communications, coding, and sparse approximation- require finite-dimensional tight frames that possess additional structural properties. This paper proposes an alternating projection method that is versatile enough to solve a huge class of inverse eigenvalue problems (IEPs), which includes the frame design problem. To apply this method, one needs only to solve a matrix nearness problem that arises naturally from the design specifications. Therefore, it is the fast and easy to develop versions of the algorithm that target new design problems. Alternating projection will often succeed even if algebraic constructions are unavailable. To demonstrate that alternating projection is an effective tool for frame design, the paper studies some important structural properties in detail. First, it addresses the most basic design problem: constructing tight frames with prescribed vector norms. Then, it discusses equiangular tight frames, which are natural dictionaries for sparse approximation. Finally, it examines tight frames whose individual vectors have low peak-to-average-power ratio (PAR), which is a valuable property for code-division multiple-access (CDMA) applications. Numerical experiments show that the proposed algorithm succeeds in each of these three cases. The appendices investigate the convergence properties of the algorithm