{'fraction_non_alphanumeric': 0.036801208720051805, 'fraction_numerical': 0.015648607813511763, 'mean_word_length': 4.574293162221776, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Robust Iterative Scheme for Symmetric Indefinite Systems

We propose a two-level nested preconditioned iterative scheme for solving sparse linear systems of equations in which the coefficient matrix is symmetric and indefinite with relatively small number of negative eigenvalues. The proposed scheme consists of an outer Minimum Residual (MINRES) iteration, preconditioned by an inner Conjugate Gradient (CG) iteration in which CG can be further preconditioned. The robustness of the proposed scheme is illustrated by solving indefinite linear systems that arise in the solution of quadratic eigenvalue problems in the context of model reduction methods for finite element models of disk brakes as well as on other problems that arise in a variety of applications.


Introduction. Symmetric indefinite linear systems
arise in many applications ranging from optimization problems to problems in computational physics, see e.g. [2,17]. In this paper we assume that A ∈ R n×n is a sparse, full-rank, symmetric and indefinite matrix with only few negative eigenvalues. Our motivation to develop a new preconditioned iterative method arises from an application in the automotive industry. In order to control brake squeal, large scale eigenvalue problems are solved via a shift-and-invert Arnoldi method to obtain a reduced model that can be used for parameter studies and optimization, see [10] and Section 3.1. We propose the use of a two-level preconditioned iterative method with a positive definite preconditioner for the solution of the arising linear systems. The basic idea of such a preconditioner iteration is well-known. In the context of optimization problems, see [9], a sparse Bunch-Parlett factorization (2) P AP T = LDL T is suggested as a solver for the systems involving the indefinite blocks of various preconditioners. Where P is a permutation matrix (with P P T = I), L is a sparse lower triangular matrix (typically with some fill-in compared to the sparsity pattern of A), and D is a block-diagonal matrix that contains either 1 × 1 or 2 × 2 blocks. Given such a factorization, one can modify the diagonal matrix D to obtain a positive definiteD such that the eigenvalues ofD are the absolute values of the eigenvalues of D, so that also M := LDL T is positive definite. If a diagonal block of D is 1 × 1 and negative, then one replaces it with its absolute value. Otherwise, it is a symmetric 2 × 2 block, and one computes the spectral decomposition where c, s ∈ R satisfy c 2 + s 2 = 1, and one replaces the 2 × 2 block with The matrix M , if easily available, is a good preconditioner for a preconditioned Krylov subspace method, such as the Minimum Residual method (MINRES) [14], since due to the fact that the spectrum of M −1 A has only the values +1, −1, it would converge in at most 2 iterations in exact arithmetic if the factorization is exact. However, this preconditioner is, in general, not practical for large problems due to fill-in and large storage requirements. In [12], therefore, an incomplete LDL T factorization (ILDL T ) based preconditioner for MINRES is proposed. Another suggestion for a preconditioner of MINRES, proposed in [21], is the positive definite absolute value of A, defined as |A| := V |Λ|V T in which A = V ΛV T is the spectral decomposition of A, however, to avoid the high computational complexity of the spectral decomposition, in [21] it is suggested to use a geometric multigrid method instead of the absolute value preconditioner and it is illustrated via a model problem that this approach is very effective when the system matrix arises from elliptic partial differential equations.
In our motivating problem, the indefinite matrix arises from a perturbed wave equation where the resulting linear system depends on parameters and has the extra property that the number of negative eigenvalues is much smaller than the number of positive eigenvalues. For this class of problems we propose a new two-level iterative scheme that combines the absolute value preconditioner approach with a deflation procedure and we show that this method is also very effective for a large class of indefinite problems arising in other applications.

2.
A two-level iterative scheme. In this section we describe a new two-level preconditioned iterative scheme for symmetric indefinite linear systems where the coefficient matrix has only very few negative eigenvalues. The method employs MIN-RES together with a modified absolute value preconditioner that is constructed via a deflation procedure which, however, is not carried out explicitly. The linear systems involving the preconditioner are solved again iteratively via the preconditioned Conjugate Gradient (CG) [7] which can be preconditioned via an incomplete LU (ILU ) decomposition, see e.g. [17], of the original coefficient matrix A or any other preconditioner obtained from the original coefficient matrix . These include but are not limited to Sparse Approximate Inverse Algebraic Multigrid based preconditioner as well. We illustrate that this MINRES-CG iterative scheme is very effective and more robust than other preconditioned general Krylov subspace methods, such as the restarted Generalized Minimum Residual (GMRES) [18], the stabilized Bi-Conjugate Gradient method (BiCGStab) [19], inner-outer FGMRES-GMRES [16] or just ILDL T preconditioned MINRES.
As an approximation to the absolute value preconditioner we use whereV is an approximate invariant subspace of A associated with the (say k) negative eigenvalues and |Λ| is the corresponding absolute value of the diagonal matrix of negative eigenvalues. Since we have assumed that k is much smaller than n, the modification (or as it is sometimes called deflation) is of small rank. In each iteration of MINRES applied to (1) a system of the form (7) M mr z = y has to be solved, and again the preconditioned matrix M −1 mr A has only eigenvalues +1 or −1 so that MINRES with the exact preconditioner converges theoretically again in at most 2 iterations. However, since M mr is symmetric and positive definite, we propose to use a preconditioned CG iteration for solving system (7) approximately with an indefinite preconditioner, M cg , which is an approximation of the original coefficient matrix itself. Note that the eigenvalues of the preconditioned matrix for CG, M −1 cg M mr , would again be either +1 or −1 if the exact matrix A −1 was used. Indefinite preconditioning for the CG method is rarely applied with the exception of [15], where CG for indefinite systems with indefinite preconditioner is used but it is assumed that the preconditioned matrix is positive definite. In our case, however, this will not be the case.
The first level preconditioner (M mr ) is symmetric and positive definite, but dense, so it should not be formed explicitly. On the other hand, the second level preconditioner (M cg ) is sparse and symmetric but not positive definite. However, the preconditioned CG is still guaranteed not to break down (see [17, p. 277]) using an indefinite preconditioner which can be seen as follows. It is well-known, see e.g. [17, p. 279], that preconditioned CG with a preconditioner M applied to a system W x = b with symmetric positive definite W can be expressed in an indefinite M -scalar product by replacing the Euclidean inner products in CG by the M -inner products. If W is symmetric positive definite, and M is symmetric indefinite (but invertible), then we can define the indefinite M -inner product as (x, y) M = (M x, y) = y T M x = x T M y = (y, x) M , so M −1 W is positive definite with respect to the M -inner product, since (M −1 W x, x) M > 0 for all x = 0.
Given the system W x = z, an initial guess x 0 , and a preconditioner M , as CG is a projection based Krylov subspace method, the vectors x m must satisfy the orthogonality condition (8) is equivalent to the orthogonality condition of CG without preconditioning (9) (z − W x m , v) = 0 for all v ∈K m Therefore, indefinitely preconditioned CG minimizes the error in the energy norm defined by the positive definite matrix W . In summary, our two-level scheme consists of two stages. First, we compute approximations to the negative eigenvalues and the corresponding invariant subspace (see Algorithm 1). This computation itself may be very expensive even if the invariant subspace has small dimension. However, in our motivating application many linear systems with the same coefficient matrix (or closely related coefficient matrices) need to be solved. Hence, this potentially expensive initial cost is quickly amortized. This is typical when solving eigenvalue problems with the shift-and-invert Arnoldi method as in [10]. The second stage in the iterative solution stage consists of nested MINRES and CG iterations (Algorithm 2). Note that while the outer MINRES iterations require matrix-vector multiplications with the original sparse coefficient matrix A, the inner CG iterations require matrix-vector multiplications of the form v = M mr u which are efficiently performed by using sparse matrix-vector multiplications and together with dense matrix-vector operations (BLAS Level 2) and vector-vector operations (BLAS Level 1) in the following procedure (11) M mr u = Au + 2(V (|Λ|(V * u))) The total cost of each such matrix multiplication operation is O(nnz + 2kn) where nnz,k and n are the number of nonzeros, negative eigenvalues and rows of A, respectively. We note that this extra operation is much more cache friendly than constructing an orthonormal basis in GMRES which rely on dot products (BLAS Level 1). Alternatively, to further speed up the convergence, the proposed scheme can be implemented using RMINRES [8,22], i.e. recycled and deflated MINRES, as the outer solver instead of plain MINRES. The trade-off would be increased storage and computation requirements due to the necessary orthogonalization against the recycled subspace and the updates of the recycled subspace [22]. Furthermore, finding subspaces that lead to improved convergence is considered to be a highly challenging task and application specific [8].
after applying the Sherman-Morrison-Woodbury formula and some algebraic manipulations, we obtain, Note that since we do not have the exact A −1 , but use an approximation of it, M −1 mr is not positive definite. Still, we can use it as the preconditioner for the CG iterations. In other words, we apply the preconditioner for CG, Solve Ax = b via MINRES using the preconditioner M mr = A + 2V |Λ|V T in each iteration of MINRES using the subroutines: • Compute matrix-vector products with A.
• Solve M mr z = y via preconditioned CG using as preconditioner M cg =Ã an approximation of A.
In each iteration of PCG: -Compute matrix-vector products: v = M mr u -Solve the system M cg t = g return x the action ofÃ −1 is approximated, such as by an incomplete factorization of A. In the improved scheme, application of the preconditioner involves an additional dense matrix-vector multiplication (BLAS Level 2) cost of 2kn but no additional storage requirement. Hereafter, we refer to this improved version as MINRES-CG * .
3. Application of the two-level method. In this section we describe the applications to which we apply the proposed two-level procedure.

Finite Element Models of Disk Brakes.
In the context of noise reduction in disk brakes, reduced order models are determined from the finite element model [10] by computing the eigenvalues in the right half plane and close to the imaginary axis of a parametric Quadratic Eigenvalue Problem (QEP) (14) ( where M and K E are symmetric positive definite, D G is skew-symmetric, D M , D R , K g are symmetric indefinite, and K R is general [10]. Here Ω denotes the angular velocity in the disk (2π < Ω < 4 × 2π) and Ω ref is the reference angular velocity. The QEP is solved by first rewriting it as a linear eigenvalue problem, using a companion linearization of (14) given by (17) 0 Audible brake squeal is associated with eigenvalues in the right half plane. For this reason we are interested in those eigenvalues that lie in a rectangular domain in the complex plane given by −50 < Re(λ) < 1, 000 and −1 < Im(λ) < 20, 000 corresponding to the audible range.
Solving the eigenvalue problem (17) via an eigensolver such as the shift-and-invert Arnoldi method [13], requires the solution of a shifted linear system of equations in each iteration, see [10] for details of the eigensolver. To apply our two-level linear system solver, we consider the solution of the following shifted linear system with complex shifts (γ inside the rectangular domain of interest), (18) C In [10] this complex linear system is solved with a sparse complex direct solver. To solve the problem iteratively, we follow [1] and map the complex system (18) to an equivalent double-size real system. Splitting into real and imaginary parts C =Â + iB and γ = γ r + iγ i with γ r = Re(γ) and γ i = Im(γ), we obtain for the real and complex parts of C, respectively. This leads to the real system which we then solve via a preconditioned Krylov subspace method with preconditioner andB =B. Note that both M and K E are symmetric and positive definite. The preconditioner can be block LU factorized as Hence, the major cost in solving systems involving the preconditioner M is the solution of two linear systems where the coefficient matrix is (i)B and (ii) S = (B +ÃB −1Ã ), namely the Schur complement. Since the solution of (i) is quite trivial, we only discuss how to solve systems involving the Schur complement matrix, which typically is dense see [2], but in our case it has the factorization Solving systems involving the Schur complement matrix, therefore, requires two steps: (i) scaling the right hand side vector with S −1 1 and (ii) solving systems where the coefficient matrix is S 2 .
Step (i) is again trivial, hence we now look into (ii) which we solve iteratively using a Krylov subspace method where the preconditioner is since in our case ||M|| F ||K E || F . Hence, the main cost in solving the block triangular systems lies in the solution of or after multiplying both sides of the system by −γ i we obtain where |γ| 2 = γ 2 i + γ 2 r . Even though M and K E are symmetric and positive definite, there is no guarantee that the symmetric coefficient matrix K E − |γ| 2 M is positive definite. However, system (28) is a perfectly suitable for the proposed MINRES-CG scheme, since in our application it only has few negative eigenvalues and they need to be computed only once. Furthermore, the preconditioner (22) is completely independent of the parameters Ω and Ω ref , and the coefficient matrix of inner systems that have to be solved (28) are the same for a given |γ|. This means that a factorization (incomplete or exact) or an approximation for the coefficient matrix K E − |γ| 2 M can be computed once and re-used for all values of γ of the same absolute value and for all corresponding Ω values.
Numerical experiments for this class of problems are presented in section 4.

Other applications.
As further applications we consider all symmetric indefinite problems in SuiteSparse Matrix Collection [5] of sizes between n = 1000 and n = 50, 000 and with at most 100 negative eigenvalues. Since this includes 7 matrices from the PARSEC group [3], we exclude the 3 smallest matrices from this group. Furthermore, since shifts around the so-called Fermi level are also of interest in the PARSEC group of matrices, we shift the largest matrix (SiO) by A−σI. For σ, we chose three values (0.25, 0.5 and 0.75) which approximately correspond to the gaps in the spectrum of A. The properties of these 11 matrices are given in Table 1. Note that we include two examples (*) that arise in finite element discretization of structural problems. These are not full eigenvalue problems but just mass matrices; solving these linear systems is useful if eigenvalues in inverse mass matrix inner product space are computed.

Numerical results.
In this section, we study the robustness of the proposed two-level scheme for indefinite linear systems described in the previous section. All experiments are performed using MATLAB R2018a.
In MINRES-CG, we use an indefinite preconditioner (M cg ) obtained either by an incomplete LDL T or LU factorization of the coefficient matrix for inner CG iterations. Former is the only suitable preconditioner available in MATLAB which we refer to as ILU . Hence, even though it does not exploit symmetry, we use it to show the robustness of the proposed scheme in Section 4.2.1. For the latter, on the other hand, we use symm-ildl which is an external package [11] that has an interface for MATLAB and is robust. Hereafter, we refer to this preconditioner as ILDL T . Therefore, we use it to show that the proposed scheme is competitive against other solvers in terms of number of iterations even when a much more robust preconditioner is used in Sections 4.1 and 4.2.2. For a fair comparison, exactly the same preconditioner is used for BiCGStab, GMRES(m) (m = 20, 40, 60 and 120) as well as another outerinner scheme with Flexible GMRES (FGMRES) as the outer solver and GMRES as the inner solver [16]. In Section 4.2.2, we also use a MINRES preconditioner with the modified ILDL T factorization. For FGMRES-GMRES we use a restart value of 120 for both inner and outer iterations. Iterative solvers, except FGMRES, are the implementations that are available in MATLAB. We note that MATLAB's BiCGStab implementation terminates early before completing a full iteration if the relative residual is already small enough. This counts as a half iteration. We modified GMRES to stop the iteration based on the true relative residual rather than the preconditioned relative residual. Storage requirements for MINRES, CG, BiCGStab, FGMRES and GMRES are given in [4,17,20,16,18], respectively. In Table 2, we illustrate the storage requirement of each of the iterative solvers via the number of vectors in addition to the coefficient matrix, the preconditioner (i.e. incomplete factors) and the right hand side vector which are common for all solvers.  4.1. Disk brake example. In the following we solve (28) for the small and large test problems of [10] of sizes n = 4, 669 and n = 842, 638, respectively, with Ω ref = 5. Note again that (28) is independent of Ω. For the first set of experiments we fix the shift γ to be the largest value in the range of values of interest, namely 1, 000 + 20, 000j. This also happens to be the most challenging case since the number of negative eigenvalues is also the largest, with k = 18 and k = 60, respectively.
For the proposed scheme, an ILDL T factorization of the coefficient matrix is used as the preconditioner (M cg ) of the inner CG iteration. We use the same preconditioner for BiCGStab, GMRES(m) and FGMRES-GMRES. For the smaller problem, we also use the ILU factorization with no fill-in (i.e. ILU (0)) preconditioner of MATLAB.
For all experiments a moderate outer stopping tolerance of relative residual norm less than or equal to 10 −3 is used. For MINRES-CG and FGMRES-GMRES schemes the inner stopping tolerance is 10 −2 . For all methods, the maximum (total) number of iterations are 2, 000 and 15, 000 for small and large problems, respectively. In all experiments, the right hand side vector is a random vector of size n.
The required number of iterations for the proposed scheme as well as for baseline algorithms are given in Table 3 for solving the small problem using ILU In Table 4, results are presented for solving the large problem using the seven iterative methods with the preconditioners LDL T (4, 10 −4 ), LDL T (5, 10 −5 ) as well as LDL T (5, 10 −6 ). Note that a much smaller dropping tolerance is required for the large problem. Incomplete factors contain 117.1 , 144.5, and 146.8 nonzeros per row which are relatively small considering complete LDL T factorization would produce 558.4 nonzeros per row. In fact, incomplete factorization may not be an efficient preconditioner for this problem. However, we still include these results here only to show the robustness of the proposed scheme in terms of number of iterations. While for all preconditioners GMRES(m), FGMRES-GMRES and BiCGStab reach the maximum number of iterations without converging, MINRES-CG still converges in 4 outer iterations albeit with a large number of inner iterations.
In Figure 1, the relative residual history is given when the ILU (0) preconditioner is used for three algorithms for the small test problem. Note that for MINRES-CG the relative residual is only available at each outer iteration. Hence, only those are presented in the figure.
As second application we fix the preconditioner to be ILU (0) and vary the shift γ in the complex domain of interest for the small test problem. Here γ is a parameter that we change in the context of the eigenvalue problem. It is of interest to see how the method behaves as γ is varied. In Figure 2, the total number of iterations is presented. Preconditioned BiCGStab fails to converge for some values of γ (shown as the white area in the figure) while MINRES-CG converges for all γ values. Figure 3 depicts the number of outer iterations and the average number of inner iterations for MINRES-CG.

4.2.
Test cases from the SuiteSparse matrix collection. In this subsection, the results are presented for systems that are obtained from the SuiteSparse Matrix Collection. In the first set we compare the proposed method against the classical general iterative schemes GMRES(m) and BiCGStab using an incomplete LU factorization based preconditioner. In the second set, we compare against the modified Table 3: Required number of iterations using various preconditioners and iterative methods for the small system ( †: maximum number of iterations is reached without convergence) Preconditioner Solver Outer its. Inner its. (Avg.) Total its.    incomplete Bunch-Parlett based preconditioned MINRES.

Comparison against ILU preconditioner.
We use ILU (0) for all cases except Meg4 where incomplete LU factorization fails due to a zero pivot. Therefore, we use the modified incomplete LU factorization in MATLAB with 10 −2 dropping tolerance (i.e. M ILU (10 −2 )) for this case only. Since in practice GMRES is always used with a value for the restart (m) we choose a restart value of m = 20, 40, 60 and 120. In FGMRES-GMRES, we use a restart of 120 for both inner and outer iterations. We stop the iterations when the relative residual norm is less than 10 −5 for all cases. The inner iteration stops when the relative residual norm is less than 10 −3 for CG and GMRES, in MINRES-CG and FGMRES-GMRES, respectively. Both in MINRES-CG and BiCGStab iterations stop when the true relative residual is less than the tolerance. For preconditioned GMRES the available residual is only the preconditioned residual. In order to have a fair comparison, we explicitly compute the true residual at each GMRES iteration and stop the iteration based on the true relative residual norm. For all methods, the maximum (total) number of iterations are 20, 000.
In Table 6, the detailed number of iterations for ILU preconditioned MINRES-CG, GMRES(m), FGMRES-GMRES and BiCGStab are given. GMRES (20) fails in 6 cases out of 11. For bcsstm10, GMRES(20) stagnates ( ‡), while for 5 other cases (namely bcsstm27, nasa1824, Si10H16, SiO(σ = 0.25) and Sio(σ = 0.75)), the maximum number of iterations is reached without convergence ( †). If the restart is increased to 40, 60 and 120, GMRES(m), fails in 4, 3 and 2 cases, respectively. BiCGStab fails for bcsstm27 and M eg4 due to the maximum of iterations being reached without convergence ( †) and a scalar quantity became too large or too small during the iteration ( * ), respectively. FGMRES-GMRES fails in 3 cases due to the maximum of iterations being reached without convergence ( †). The proposed MINRES-CG method does not fail in any of the test problems. Although the cost per iteration is different for each method, the total number of iterations are presented in Table 7. For the cases they do not fail, GMRES(120) and FGMRES-GMRES requires fewer number of iterations than MINRES-CG but they also require more storage. In 4 cases MINRES-CG requires fewer iterations than BiCGStab. It is possible to improve the total number of iterations of MINRES-CG via using the algorithm described in Section 2.1, Table 5 shows the improved number of iterations which is significant especially for the cases where the inner or outer number of iterations are high.
In order to study the effect of the inner stopping tolerance on the eigenvalues of the preconditioned matrix, we explicitly compute M −1 mr A using preconditioend CG iterations using stopping tolerances of 10 −2 , 10 −3 and 10 −4 . In Figure 5, a clear clustering of eigenvalues of the preconditioned matrix M −1 mr A is visible around +1 and −1 for bcsstm10 while the unpreconditioned coefficient matrix had no clustering of eigenvalues (see Figure 4). As expected, the clustering around −1 and +1 improves as the stopping tolerance for the inner CG iterations is decreased. In this subsection, we compare the proposed scheme against a robust incomplete LDL T (Bunch-Parlett) factorization [6].
We use the ILDL T implementation of [11] in MATLAB which computes the incomplete Bunch-Parlett factorization of the coefficient matrix. The default parameters are 3 and 10 −3 for the level of fill-in and the dropping tolerance, respectively. Furthermore, it uses the Approximate Minimum Degree reordering, Rook pivoting and Scaling are used to improve the numerical stability of the incomplete factors by default. Note that all of those enhancements that are implemented in ILDL T makes the preconditioner much more robust than the ILU (0) preconditioner. In the following experiments all methods are applied to the permuted and scaled linear systems.
After computing the ILDL T factorization the D matrix is modified as described in [9] in order to obtain a positive definite preconditioner to be used with MINRES. To have a fair comparison, the same ILDL T factorization (without the modification) is used as the preconditioner for GMRES(m), FGMRES-GMRES, BiCGStab and as the inner preconditioner for MINRES-CG. Stopping tolerances and the maximum number of iterations allowed are set exactly the same as in Section 4.2.
In Table 8 the total number of iterations for all methods are given. Even though it is a much more robust preconditioner, MINRES preconditioned with the modified ILDL T preconditioner stagnates ( ‡) for bcsstm27. For the same problem BiCGStab, FGMRES-GMRES and GMRES(m) (for all restart values m = 20, 40, 60, 120) reach the maximum number of iterations without converging ( †). On the other hand, MINRES-CG converges in all problems which confirms the robustness of the proposed scheme. In Table 9, the total number of iterations for all methods are given. GMRES(m) with larger restart values and BiCGStab require the fewest number of iterations. FGMRES-GMRES requires fewer number of iterations than MINRES-CG.  On the other hand, MINRES requires more iterations than MINRES-CG in 4 cases, and the required number of iterations are marginally better than that of MINRES-CG for 3 other cases.

Conclusions.
A two-level nested iterative scheme is proposed for solving sparse linear systems of equations where the coefficient matrix is symmetric indefinite with few negative eigenvalues. The first level is MINRES preconditioned via CG. The inner level CG is preconditioned via the original indefinite coefficient matrix. The robustness of the proposed scheme is presented for linear systems that arise in disk brake squeal as well as systems that arise in a variety of test cases from the SuiteSparse Matrix Collection.