Fast and stable two‐way algorithm for diagonal plus semi‐separable systems of linear equations

An algorithm for the solution of linear systems of equations where the coefficient matrix is diagonal plus a semi‐separable matrix is considered. The algorithm is stable with linear complexity. Furthermore, it is suitable for an implementation on a system of two processors. Copyright © 2001 John Wiley & Sons, Ltd.     

The retraction algorithm for factoring banded symmetric matrices

Let A be an n × n symmetric matrix of bandwidth 2m + 1. The matrix need not be positive definite. In this paper we will present an algorithm for factoring A which preserves symmetry and the band structure and limits element growth in the factorization. With this factorization one may solve a linear system with A as the coefficient matrix and determine the inertia of A, the number of positive, negative, and zero eigenvalues of A. The algorithm requires between 1/2nm² and 5/4nm² multiplications and at most (2m + 1)n locations compared to non‐symmetric Gaussian elimination which requires between nm² and 2nm² multiplications and at most (3m + 1)n locations. Our algorithm reduces A to block diagonal form with 1 × 1 and 2 × 2 blocks on the diagonal. When pivoting for stability and subsequent transformations produce non‐zero elements outside the original band, column/row transformations are used to retract the bandwidth. To decrease the operation count and the necessary storage, we use the fact that the correction outside the band is rank‐1 and invert the process, applying the transformations that would restore the bandwidth first, followed by a modified correction. This paper contains an element growth analysis and a computational comparison with LAPACKs non‐symmetric band routines and the Snap‐back code of Irony and Toledo. Copyright © 2007 John Wiley & Sons, Ltd.     

An implicit QR algorithm for symmetric semiseparable matrices

The QR algorithm is one of the classical methods to compute the eigendecomposition of a matrix. If it is applied on a dense n × n matrix, this algorithm requires O(n³) operations per iteration step. To reduce this complexity for a symmetric matrix to O(n), the original matrix is first reduced to tridiagonal form using orthogonal similarity transformations. In the report (Report TW360, May 2003) a reduction from a symmetric matrix into a similar semiseparable one is described. In this paper a QR algorithm to compute the eigenvalues of semiseparable matrices is designed where each iteration step requires O(n) operations. Hence, combined with the reduction to semiseparable form, the eigenvalues of symmetric matrices can be computed via intermediate semiseparable matrices, instead of tridiagonal ones. The eigenvectors of the intermediate semiseparable matrix will be computed by applying inverse iteration to this matrix. This will be achieved by using an O(n) system solver, for semiseparable matrices. A combination of the previous steps leads to an algorithm for computing the eigenvalue decompositions of semiseparable matrices. Combined with the reduction of a symmetric matrix towards semiseparable form, this algorithm can also be used to calculate the eigenvalue decomposition of symmetric matrices. The presented algorithm has the same order of complexity as the tridiagonal approach, but has larger lower order terms. Numerical experiments illustrate the complexity and the numerical accuracy of the proposed method. Copyright © 2005 John Wiley & Sons, Ltd.     

New fast divide‐and‐conquer algorithms for the symmetric tridiagonal eigenvalue problem

In this paper, two accelerated divide‐and‐conquer (ADC) algorithms are proposed for the symmetric tridiagonal eigenvalue problem, which cost O(N²r) flops in the worst case, where N is the dimension of the matrix and r is a modest number depending on the distribution of eigenvalues. Both of these algorithms use hierarchically semiseparable (HSS) matrices to approximate some intermediate eigenvector matrices, which are Cauchy‐like matrices and are off‐diagonally low‐rank. The difference of these two versions lies in using different HSS construction algorithms, one (denoted by ADC1) uses a structured low‐rank approximation method and the other (ADC2) uses a randomized HSS construction algorithm. For the ADC2 algorithm, a method is proposed to estimate the off‐diagonal rank. Numerous experiments have been carried out to show their stability and efficiency. These algorithms are implemented in parallel in a shared memory environment, and some parallel implementation details are included. Comparing the ADCs with highly optimized multithreaded libraries such as Intel MKL, we find that ADCs could be more than six times faster for some large matrices with few deflations. Copyright © 2016 John Wiley & Sons, Ltd.     

A Low-Complexity Divide-and-Conquer Method for Computing Eigenvalues and Eigenvectors of Symmetric Band Matrices

A framework for an efficient low-complexity divide-and-conquer algorithm for computing eigenvalues and eigenvectors of an n × n symmetric band matrix with semibandwidth b n is proposed and its arithmetic complexity analyzed. The distinctive feature of the algorithm—after subdivision of the original problem into p subproblems and their solution—is a separation of the eigenvalue and eigenvector computations in the central synthesis problem. The eigenvalues are computed recursively by representing the corresponding symmetric rank b(p–1) modification of a diagonal matrix as a series of rank-one modifications. Each rank-one modifications problem can be solved using techniques developed for the tridiagonal divide-and-conquer algorithm. Once the eigenvalues are known, the corresponding eigenvectors can be computed efficiently using modified QR factorizations with restricted column pivoting. It is shown that the complexity of the resulting divide-and-conquer algorithm is O (n ² b ²) (in exact arithmetic).This revised version was published online in October 2005 with corrections to the Cover Date.     

Parallel Dixon Matrices by Bracket

It is known that the Dixon matrix can be constructed in parallel either by entry or by diagonal. This paper presents another parallel matrix construction, this time by bracket. The parallel by bracket algorithm is the fastest among the three, but not surprisingly it requires the highest number of processors. The method also shows analytically that the Dixon matrix has a total of m(m+1)²(m+2)n(n+1)²(n+2)/36 brackets but only mn(m+1)(n+1)(mn+2m+2n+1)/6 of them are distinct.     

A fast algorithm for Gaussian type quadrature formulae with mixed boundary conditions and some lumped mass spectral approximations

After studying Gaussian type quadrature formulae with mixed boundary conditions, we suggest a fast algorithm for computing their nodes and weights. It is shown that the latter are computed in the same manner as in the theory of the classical Gauss quadrature formulae. In fact, all nodes and weights are again computed as eigenvalues and eigenvectors of a real symmetric tridiagonal matrix. Hence, we can adapt existing procedures for generating such quadrature formulae. Comparative results with various methods now in use are given. In the second part of this paper, new algorithms for spectral approximations for second-order elliptic problems are derived. The key to the efficiency of our algorithms is to find an appropriate spectral approximation by using the most accurate quadrature formula, which takes the boundary conditions into account in such a way that the resulting discrete system has a diagonal mass matrix. Hence, our algorithms can be used to introduce explicit resolutions for the time-dependent problems. This is the so-called lumped mass method. The performance of the approach is illustrated with several numerical examples in one and two space dimensions.

     

Symmetrizing a Hessenberg matrix: Designs for VLSI parallel processor arrays

A symmetrizer of a nonsymmetric matrix A is the symmetric matrixX that satisfies the equationXA =A ^tX, wheret indicates the transpose. A symmetrizer is useful in converting a nonsymmetric eigenvalue problem into a symmetric one which is relatively easy to solve and finds applications in stability problems in control theory and in the study of general matrices. Three designs based on VLSI parallel processor arrays are presented to compute a symmetrizer of a lower Hessenberg matrix. Their scope is discussed. The first one is the Leiserson systolic design while the remaining two, viz., the double pipe design and the fitted diagonal design are the derived versions of the first design with improved performance.     

A divide-and-conquer algorithm for the eigendecomposition of symmetric block-diagonal plus semiseparable matrices

Summary. We present a fast and numerically stable algorithm for computing the eigendecomposition of a symmetric block diagonal plus semiseparable matrix. We report numerical experiments that indicate that our algorithm is significantly faster than the standard method which treats the given matrix as a general symmetric dense matrix. Mathematics Subject Classification (1991):15A09, 15A23, 65F05, 65L10, 65R20This research was supported in part by NSF Career Award CCR-9734290.This research was supported in part by NSF Career Award CCR-9702866 and by Alfred Sloan Research Fellowship BR-3720.Received: 10, September 2001     

A Jacobi eigenreduction algorithm for definite matrix pairs

Summary We propose a Jacobi eigenreduction algorithm for symmetric definite matrix pairsA, J of small to medium-size real symmetric matrices withJ ²=I,J diagonal (neitherJ norA itself need be definite). Our Jacobi reduction works only on one matrix and usesJ-orthogonal elementary congruences which include both trigonometric and hyperbolic rotations and preserve the symmetry throughout the process. For the rotation parameters only the pivotal elements of the current matrix are needed which facilitates parallelization. We prove the global convergence of the method; the quadratic convergence was observed in all experiments. We apply our method in two situations: (i) eigenreducing a single real symmetric matrix and (ii) eigenreducing an overdamped quadratic matrix pencil. In both cases our method is preceded by a symmetric indefinite decomposition and performed in its one-sided variant on the thus obtained factors. Our method outdoes the standard methods like standard Jacobi orqr/ql in accuracy in spite of the use of hyperbolic transformations which are not orthogonal (a theoretical justification of this behaviour is made elsewhere). The accuracy advantage of our method can be particularly drastic if the eigenvalues are of different order. In addition, in working with quadratic pencils our method is shown to either converge or to detect non-overdampedness.     

Finding off‐diagonal entries of the inverse of a large symmetric sparse matrix

The method fast inverse using nested dissection (FIND) was proposed to calculate the diagonal entries of the inverse of a large sparse symmetric matrix. In this paper, we show how the FIND algorithm can be generalized to calculate off‐diagonal entries of the inverse that correspond to ‘short’ geometric distances within the computational mesh of the original matrix. The idea is to extend the downward pass in FIND that eliminates all nodes outside of each node cluster. In our advanced downwards pass, it eliminates all nodes outside of each ‘node cluster pair’ from a subset of all node cluster pairs. The complexity depends on how far (i,j) is from the main diagonal. In the extension of the algorithm, all entries of the inverse that correspond to vertex pairs that are geometrically closer than a predefined length limit l will be calculated. More precisely, let α be the total number of nodes in a two‐dimensional square mesh. We will show that our algorithm can compute O(α^{3 ∕ 2 + 2ε}) entries of the inverse in O(α^{3 ∕ 2 + 2ε}) time where l = O(α^{1 ∕ 4 + ε}) and 0 ≤ ε ≤1 ∕ 4. Numerical examples are given to illustrate the efficiency of the proposed method. Copyright © 2012 John Wiley & Sons, Ltd.     

Orthogonal similarity transformation into block‐semiseparable matrices of semiseparability rank k

Very recently, an algorithm, which reduces any symmetric matrix into a semiseparable one of semi‐ separability rank 1 by similar orthogonality transformations, has been proposed by Vandebril, Van Barel and Mastronardi. Partial execution of this algorithm computes a semiseparable matrix whose eigenvalues are the Ritz‐values obtained by the Lanczos' process applied to the original matrix. Also a kind of nested subspace iteration is performed at each step. In this paper, we generalize the above results and propose an algorithm to reduce any symmetric matrix into a similar block‐semiseparable one of semiseparability rank k, with k ∈ ?, by orthogonal similarity transformations. Also in this case partial execution of the algorithm computes a block‐semiseparable matrix whose eigenvalues are the Ritz‐values obtained by the block‐Lanczos' process with k starting vectors, applied to the original matrix. Subspace iteration is performed at each step as well. Copyright © 2005 John Wiley & Sons, Ltd.     

Efficient Ritz-Galerkin Discretization of PDEs with Variable Coefficients in arbitrary Dimensions using Sparse Grids

Sparse grids can be used to discretize elliptic differential equations of second order on a d-dimensional cube. Using the Ritz-Galerkin discretization, one obtains a linear equation system with 𝒪 (N (log N)^d−1) unknowns. The corresponding discretization error is 𝒪 (N⁻¹ (log N)^d−1) in the H¹-norm. A major difficulty in using this sparse grid discretization is the complexity of the related stiffness matrix. To reduce the complexity of the sparse grid discretization matrix, we apply prewavelets and a discretization with semi-orthogonality. Furthermore, a recursive algorithm is used, which performs a matrix vector multiplication with the stiffness matrix by 𝒪 (N (log N)^d−1) operations. Simulation results up to level 10 are presented for a 6-dimensional Helmholtz problem with variable coefficients. (© 2017 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Compatible wavelet-Galerkin method to solve stokes interior problem with circular boundary

This article addresses Neumann boundary value interior problem of Stokes equations with circular boundary. By using natural boundary element method, the Stokes interior problem is reduced into an equivalent natural integral equation with a hyper-singular kernel, which is viewed as Hadamard finite part. Based on trigonometric wavelet functions, the compatible wavelet space is constructed so that it can serve as Galerkin trial function space. In proposed compatible wavelet-Galerkin method, the simple and accurate computational formulae of the entries in stiffness matrix are obtained by singularity removal technique. It is also proved that the stiffness matrix is almost a block diagonal matrix, and its diagonal sub-blocks all are both symmetric and circulant submatrices. These good properties indicate that a 2 ^J+3 × 2 ^J+3 stiffness matrix can be determined only by its 2 ^J + 3J + 1 entries. It greatly decreases the computational complexity. Some error estimates for the compatible wavelet-Galerkin projection solutions are established. Finally, several numerical examples are given to demonstrate the validity of the proposed approach.     

Generalized qd algorithm for block band matrices

The generalized qd algorithm for block band matrices is an extension of the block qd algorithm applied to a block tridiagonal matrix. This algorithm is applied to a positive definite symmetric block band matrix. The result concerning the behavior of the eigenvalues of the first and the last diagonal block of the matrix containing the entries q ^(k) which was obtained in the tridiagonal case is still valid for positive definite symmetric block band matrices. The eigenvalues of the first block constitute strictly increasing sequences and those of the last block constitute strictly decreasing sequences. The theorem of convergence, given in Draux and Sadik (Appl Numer Math 60:1300?C1308, 2010), also remains valid in this more general case.     

Orthogonal similarity transformation of a symmetric matrix into a diagonal-plus-semiseparable one with free choice of the diagonal

In this paper we describe an orthogonal similarity transformation for transforming arbitrary symmetric matrices into a diagonal-plus-semiseparable matrix, where we can freely choose the diagonal. Very recently an algorithm was proposed for transforming arbitrary symmetric matrices into similar semiseparable ones. This reduction is strongly connected to the reduction to tridiagonal form. The class of semiseparable matrices can be considered as a subclass of the diagonalplus- semiseparable matrices. Therefore we can interpret the proposed algorithm here as an extension of the reduction to semiseparable form. A numerical experiment is performed comparing thereby the accuracy of this reduction algorithm with respect to the accuracy of the traditional reduction to tridiagonal form, and the reduction to semiseparable form. The experiment indicates that all three reduction algorithms are equally accurate. Moreover it is shown in the experiments that asymptotically all the three approaches have the same complexity, i.e. that they have the same factor preceding the n³ term in the computational complexity. Finally we illustrate that special choices of the diagonal create a specific convergence behavior. The research was partially supported by the Research Council K.U.Leuven, project OT/05/40 (Large rank structured matrix computations), by the Fund for Scientific Research–Flanders (Belgium), projects G.0078.01 (SMA: Structured Matrices and their Applications), G.0176.02 (ANCILA: Asymptotic aNalysis of the Convergence behavior of Iterative methods in numerical Linear Algebra), G.0184.02 (CORFU: Constructive study of Orthogonal Functions) and G.0455.0 (RHPH: Riemann-Hilbert problems, random matrices and Padé-Hermite approximation), and by the Belgian Programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister's Office for Science, Technology and Culture, project IUAP V-22 (Dynamical Systems and Control: Computation, Identification & Modelling). The scientific responsibility rests with the authors.     

Complexity reduction of least squares problems involving special vandermonde matrices

This paper presents a new QRD factorization of a rectangular Vandermonde matrix for a special point distribution, including the symmetric case, based on ak-dimensional block decomposition of the matrix and some properties of the Kronecker product. The computational reduction factor with respect to any QR method isk ², in the general case, and 4 in the symmetric one. By the resulting matrix factorization, new formulas are devised for the least squares system solution, whose implementation produces an algorithm of reduced computational cost and computer storage. Finally the perturbation bounds of this new factorization are devised.     

The Inertia Set of Nonnegative Symmetric Sign Pattern with Zero Diagonal

The inertia set of a symmetric sign pattern A is the set i(A) = {i(B) | B = B ^T ∈ Q(A)}, where i(B) denotes the inertia of real symmetric matrix B, and Q(A) denotes the sign pattern class of A. In this paper, a complete characterization on the inertia set of the nonnegative symmetric sign pattern A in which each diagonal entry is zero and all off-diagonal entries are positive is obtained. Further, we also consider the bound for the numbers of nonzero entries in the nonnegative symmetric sign patterns A with zero diagonal that require unique inertia.