Symmetric Indefinite Matrix Research Articles

Conditioning analysis for discrete Helmholtz problems

In this paper, we examine conditioning of the discretization of the Helmholtz problem. Although the discrete Helmholtz problem has been studied from different perspectives, to the best of our knowledge, there is no conditioning analysis for it. We aim to fill this gap in the literature. We propose a novel method in 1D to observe the near-zero eigenvalues of a symmetric indefinite matrix. Standard classification of ill-conditioning based on the matrix condition number is not true for the discrete Helmholtz problem. We relate the ill-conditioning of the discretization of the Helmholtz problem with the condition number of the matrix. We carry out analytical conditioning analysis in 1D and extend our observations to 2D with numerical observations. We examine several discretizations. We find different regions in which the condition number of the problem shows different characteristics. We also explain the general behavior of the solutions in these regions.

An efficient relaxed shift-splitting preconditioner for a class of complex symmetric indefinite linear systems

<abstract><p>In this work, by introducing a scalar matrix $ \alpha I $, we transform the complex symmetric indefinite linear systems $ (W+i T)x = b $ into a block two-by-two complex equations equivalently, and propose an efficient relaxed shift-splitting (ERSS) preconditioner. By adopting the relaxation technique, the ERSS preconditioner is not only a computational advantage but also closer to the original two-by-two of complex coefficient matrix. The eigenvalue distributions of the preconditioned matrix are analysed. An efficient and practical formula for computing the parameter value $ \alpha $ is also derived by computing the Frobenius norm of symmetric indefinite matrix $ T $. Numerical examples on a few model problems are illustrated to verify the performances of the ERSS preconditioner.</p></abstract>

Modified projected Newton scheme for non-convex function with simple constraints

In this paper, a descent line search scheme is proposed to find a local minimum point of a non-convex optimization problem with simple constraints. The idea ensures that the scheme escapes the saddle points and finally settles for a local minimum point of the non-convex optimization problem. A positive definite scaling matrix for the proposed scheme is formed through symmetric indefinite matrix factorization of the Hessian matrix of the objective function at each iteration. A numerical illustration is provided, and the global convergence of the scheme is also justified.

Preconditioned iterative methods for the convective FitzHugh–Nagumo equations

We consider preconditioned iterative methods for numerical solution of the convective FitzHugh–Nagumo equations. By eliminating the symmetric indefinite matrix derived from the single nonlinear term of the convective FitzHugh–Nagumo equations, a preconditioner that is robust with respect to the involved model parameters is proposed to solve the arising discretized linear system. Spectral properties of the corresponding preconditioned system are analyzed. Numerical experiments are presented to illustrate the effectiveness of the proposed preconditioner.

Efficient Implementation of Second-Order Stochastic Approximation Algorithms in High-Dimensional Problems.

Stochastic approximation (SA) algorithms have been widely applied in minimization problems when the loss functions and/or the gradient information are only accessible through noisy evaluations. Stochastic gradient (SG) descent-a first-order algorithm and a workhorse of much machine learning-is perhaps the most famous form of SA. Among all SA algorithms, the second-order simultaneous perturbation stochastic approximation (2SPSA) and the second-order stochastic gradient (2SG) are particularly efficient in handling high-dimensional problems, covering both gradient-free and gradient-based scenarios. However, due to the necessary matrix operations, the per-iteration floating-point-operations (FLOPs) cost of the standard 2SPSA/2SG is O(p3) , where p is the dimension of the underlying parameter. Note that the O(p3) FLOPs cost is distinct from the classical SPSA-based per-iteration O(1) cost in terms of the number of noisy function evaluations. In this work, we propose a technique to efficiently implement the 2SPSA/2SG algorithms via the symmetric indefinite matrix factorization and show that the FLOPs cost is reduced from O(p3) to O(p2) . The formal almost sure convergence and rate of convergence for the newly proposed approach are directly inherited from the standard 2SPSA/2SG. The improvement in efficiency and numerical stability is demonstrated in two numerical studies.

On the growth factor upper bound for Aasen’s algorithm

Aasen’s algorithm factorizes a symmetric indefinite matrix A as A=PTLTLTP, where P is a permutation matrix, L is unit lower triangular with its first column being the first column of the identity matrix, and T is tridiagonal. In this note, we provide a growth factor upper bound for Aasen’s algorithm which is much smaller than that given by Higham. We also show that the upper bound we have given is not tight when the matrix dimension is greater than or equal to 6.

Non‐GPU‐resident symmetric indefinite factorization

SummaryWe study various algorithms to factorize a symmetric indefinite matrix that does not fit in the core memory of a computer. There are two sources of the data movement into the memory: one needed for selecting and applying pivots and the other needed to update each column of the matrix for the factorization. It is a challenge to obtain high performance of such an algorithm when the pivoting is required to ensure the numerical stability of the factorization. For example, when factorizing each column of the matrix, a diagonal entry, which ensures the stability, may need to be selected as a pivot among the remaining diagonals, and moved to the leading diagonal by swapping both the corresponding rows and columns of the matrix. If the pivot is not in the core memory, then it must be loaded into the core memory. For updating the matrix, the data locality may be improved by partitioning the matrix. For example, a right‐looking partitioned algorithm first factorizes the leading columns, called panel, and then uses the factorized panel to update the trailing submatrix. This algorithm only accesses the trailing submatrix after each panel factorization (instead of after each column factorization) and performs most of its floating‐point operations (flops) using BLAS‐3, which can take advantage of the memory hierarchy. However, because the pivots cannot be predetermined, the whole trailing submatrix must be updated before the next panel factorization can start. When the whole submatrix does not fit in the core memory all at once, loading the block columns into the memory can become the performance bottleneck. Similarly, the left‐looking variant of the algorithm would require to update each panel with all of the previously factorized columns. This makes it a much greater challenge to implement an efficient out‐of‐core symmetric indefinite factorization compared with an out‐of‐core nonsymmetric LU factorization with partial pivoting, which only requires to swap the rows of the matrix and accesses the trailing submatrix after each in‐core factorization (instead of after each panel factorization by the symmetric factorization). To reduce the amount of the data transfer, in this paper we uses the recently proposed left‐looking communication‐avoiding variant of the symmetric factorization algorithm to factorize the columns in the core memory, and then perform the partitioned right‐looking out‐of‐core trailing submatrix updates. This combination may still require to load the pivots into the core memory, but it only updates the trailing submatrix after each in‐core factorization, while the previous algorithm updates it after each panel factorization.Although these in‐core and out‐of‐core algorithms can be applied at any level of the memory hierarchy, we apply our designs to the GPU and CPU memory, respectively. We call this specific implementation of the algorithm a non–GPU‐resident implementation. Our performance results on the current hybrid CPU/GPU architecture demonstrate that when the matrix is much larger than the GPU memory, the proposed algorithm can obtain significant speedups over the communication‐hiding implementations of the previous algorithms.

Cholesky-Like Factorization of Symmetric Indefinite Matrices and Orthogonalization with Respect to Bilinear Forms

It is well known that orthogonalization of column vectors in a rectangular matrix $B$ with respect to the bilinear form induced by a nonsingular symmetric indefinite matrix $A$ can be eventually seen as its factorization $B=QR$ that is equivalent to the Cholesky-like factorization in the form $B^TAB=R^T \Omega R$, where $R$ is upper triangular and $\Omega$ is a signature matrix. Under the assumption of nonzero principal minors of the matrix $M=B^T A B$ we give bounds for the conditioning of the triangular factor $R$ in terms of extremal singular values of $M$ and of only those principal submatrices of $M$ where there is a change of sign in $\Omega$. Using these results we study the numerical behavior of two types of orthogonalization schemes and we give the worst-case bounds for quantities computed in finite precision arithmetic. In particular, we analyze the implementation based on the Cholesky-like factorization of $M$ and the Gram--Schmidt process with respect to the bilinear form induced by the matrix...

The Antitriangular Factorization of Symmetric Matrices

Indefinite symmetric matrices occur in many applications, such as optimization, least squares problems, partial differential equations, and variational problems. In these applications one is often interested in computing a factorization of the indefinite matrix that puts into evidence the inertia of the matrix or possibly provides an estimate of its eigenvalues. In this paper we propose an algorithm that provides this information for any symmetric indefinite matrix by transforming it to a block antitriangular form using orthogonal similarity transformations. We also show that the algorithm is backward stable and has a complexity that is comparable to existing matrix decompositions for dense indefinite matrices.

Matrix pencils and existence conditions for quadratic programming with a sign-indefinite quadratic equality constraint

We consider minimization of a quadratic objective function subject to a sign-indefinite quadratic equality constraint. We derive necessary and sufficient conditions for the existence of solutions to the constrained minimization problem. These conditions involve a generalized eigenvalue of the matrix pencil consisting of a symmetric positive-semidefinite matrix and a symmetric indefinite matrix. A complete characterization of the solution set to the constrained minimization problem in terms of the eigenspace of the matrix pencil is provided.

Tridiagonal-Diagonal Reduction of Symmetric Indefinite Pairs

We consider the reduction of a symmetric indefinite matrix pair (A,B), with B nonsingular, to tridiagonal-diagonal form by congruence transformations. This is an important reduction in solving polynomial eigenvalue problems with symmetric coefficient matrices and in frequency response computations. The pair is first reduced to symmetric-diagonal form. We describe three methods for reducing the symmetric-diagonal pair to tridiagonal-diagonal form. Two of them employ more stable versions of Brebner and Grad's pseudosymmetric Givens and pseudosymmetric Householder reductions, while the third is new and based on a combination of Householder reflectors and hyperbolic rotations. We prove an optimality condition for the transformations used in the third reduction. We present numerical experiments that compare the different approaches and show improvements over Brebner and Grad's reductions.

Perturbation analysis for the generalized Cholesky factorization

Let K be a symmetric indefinite matrix. Suppose that K= LJL T is the generalized Cholesky factorization of K. In this paper we present perturbation analysis for the generalized Cholesky factorization. We obtain the first-order bound on the norm of the perturbation in the generalized Cholesky factor. Also, we give rigorous perturbation bounds.

Inertia-controlling factorizations for optimization algorithms

When solving smooth optimization problems by sequential quadratic programming methods or interior methods using exact second-derivative information, often a linear system of equations involving a large sparse symmetric indefinite matrix has to be solved at each iteration. This matrix is often referred to as the barrier KKT matrix or the KKT matrix, depending on if an interior approach is used or not. For simplicity, we refer to it as the KKT matrix below. The inertia of the KKT matrix reveals if the problem is locally convex at the current iterate. In the convex case, the KKT matrix always has the correct inertia, and in this situation an off-the shelf sparse factorization routine can be used, where the pivot selection is based on sparsity and numerical stability. In the nonconvex case, it is of interest to find out the inertia during the factorization process so that the matrix may be modified a posteriori. This can be done by imposing an additional restriction on the pivot selection, based on controlling the inertia of the submatrix factorized so far. We describe the issues of an inertia-controlling factorization and describe our experiments with a particular sparse implementation.

A Note on the LDLT Decomposition of Matrices from Saddle-Point Problems

Sparse linear systems Kx=b are considered, where K is a specially structured symmetric indefinite matrix. These systems arise frequently, e.g., from mixed finite element discretizations of PDE problems. The LDLT factorization of K with diagonal D and unit lower triangular L is known to exist for natural ordering of K, but the resulting triangular factors can be rather dense. On the other hand, for a given permutation matrix P, the LDLT factorization of PTKP may not exist. In this paper a new way to obtain a fill-in minimizing permutation based on initial fill-in minimizing ordering is introduced. For an important subclass of matrices arising from mixed and hybrid finite element discretizations, the existence of the LDLT factorization of the permuted matrix is proved. Experimental results on practical problems indicate that the amount of computational savings can be substantial when compared with the approach based on Schur complement.

A hybrid iterative method for symmetric indefinite linear systems

Hybrid iterative methods that combine a conjugate direction method with a simpler iteration scheme, such as Chebyshev or Richardson iteration, were first proposed in the 1950s. The ease with which Chebyshev and Richardson iteration can be implemented efficiently on a large variety of computer architectures has in recent years lead to renewed interest in iterative methods that use Chebyshev or Richardson iteration. This paper presents a new hybrid iterative method for the solution of linear systems of equations with a symmetric indefinite matrix. Our method combines the conjugate residual method with Richardson iteration. Special attention is paid to the determination of two real intervals, one on each side of the origin, that contain most of the eigenvalues of the matrix. These intervals are used to compute suitable iteration parameters for Richardson iteration. We also discuss when to switch between the methods. The hybrid scheme typically uses the Richardson method for most iterations, and this reduces the number of arithmetic vector operations significantly compared with the number of arithmetic vector operations required when only the conjugate residual method is used. Computed examples illustrate the competitiveness of the hybrid scheme.

The Potential Fluid Flow Problem and the Convergence Rate of the Minimal Residual Method

In the paper the potential fluic flow problem in porous media using Darcy's law and the continuity equation is solved. Mixed-hybrid finite element formulation based on general trilateral prismatic elements is considered. Spectral properties of resulting symmetric indefinite system of linear equations are examined. Minimal residual method for the solution of systems with a symmetric indefinite matrix is applied. The rate of convergence and the asymptotic convergence factor which depend on the eigenvalue distribution of the system matrix are estimated.

On the Stability of Cholesky Factorization for Symmetric Quasidefinite Systems

Sparse linear equations $Kd = r$ are considered, where K is a specially structured symmetric indefinite matrix that arises in numerical optimization and elsewhere. Under certain conditions, K is quasidefinite. The Cholesky factorization $PKP^T = LDL^T $ is then known to exist for any permutation P, even though D is indefinite. Quasidefinite matrices have been used successfully by Vanderbei within barrier methods for linear and quadratic programming. An advantage is that for a sequence of K’s, P may be chosen once and for all to optimize the sparsity of L, as in the positive-definite case. A preliminary stability analysis is developed here. It is observed that a quasidefinite matrix is closely related to an unsymmetric positive-definite matrix, for which an $LDM^T $ factorization exists. Using the Golub and Van Loan analysis of the latter, conditions are derived under which Cholesky factorization is stable for quasidefinite systems. Some numerical results confirm the predictions.

Newton Methods for Large-Scale Linear Equality-Constrained Minimization

Newton methods for large-scale minimization subject to linear equality constraints are discussed. For large-scale problems, it may be prohibitively expensive to reduce the problem to an unconstrained problem in the null space of the constraint matrix. We investigate computational schemes that enable the computation of descent directions and directions of negative curvature without the need to know the null-space matrix. The schemes are based on factorizing a sparse symmetric indefinite matrix. Three different methods are proposed based on the schemes described for computing the search directions. Convergence properties for the methods are established.

The Lanczos algorithm applied to Kron's method

AbstractIn this paper Kron's primitive composite system matrix is shown to be reducible to a symmetric indefinite matrix. This matrix, although never formed explicitly, can then be transformed to a tridiagonal matrix of reduced order by the Lanczos algorithm. Eigenvalue solutions of this partially tridiagonalized matrix give very good approximations to the eigenvalues of the composite system. The method has multiplication counts which are over 90 per cent lower than when the (condensed) Kron matrix is solved by the Newton‐Raphson iteration applied to the Kron scalar equation. Numerical examples illustrating the Kron‐Lanczos method when solving for the natural frequencies of frames with large numbers of degrees of freedom are presented.

On the use of directions of negative curvature in a modified newton method

We present a modified Newton method for the unconstrained minimization problem. The modification occurs in non-convex regions where the information contained in the negative eigenvalues of the Hessian is taken into account by performing a line search along a path which is initially tangent to a direction of negative curvature. We give termination criteria for the line search and prove that the resulting iterates are guaranteed to converge, under reasonable conditions, to a critical point at which the Hessian is positive semidefinite. We also show how the Bunch and Parlett decomposition of a symmetric indefinite matrix can be used to give entirely adequate directions of negative curvature.