An efficient approach for solving saddle point problems using block structure

Saddle point problems appear in many scientific and engineering applications. Hence, it is important to solve them fast with high accuracy. We have proposed the method for saddle point problems using block structure in order to solve them fast. In our proposed method, first, a linear system with multiple right-hand sides is solved instead of solving the saddle point problems directly. After that, a small linear system with a dense matrix is solved. In our previous work, it has been observed that our method is faster than the conventional approach, but the accuracy of the obtained approximate solutions is worse than the conventional one. In this paper, we propose to improve the accuracy of the approximate solutions of the saddle point problems by improving the accuracy of the solution of the small linear system using the mixed precision iterative refinement technique. Numerical experiments illustrate that the proposed approach improves the accuracy of the approximate solutions to the same extent as the conventional approach.

Block Krylov subspace methods are efficient iterative methods for solving linearsystems with multiple right-hand sides. The Block BiCGGR method is one of the Block Krylov subspace methods, and can generate highly accurate approximatesolutions. However, when the number of right-hand sides is large, this method becomes numerically unstable. In this paper, the Block BiCGGR2 method is proposed. This method is a numerically stable and highly accurate method than the Block BiCGGR method and the Block BiCGSTAB method. Numerical results illustrate the robustness and the accuracy of the approximate solutions of the Block BiCGGR2 method.

For Hermitian positive definite linear systems and eigenvalue problems, the eigCG algorithm is a memory efficient algorithm that solves the linear system and simultaneously computes some of its eigenvalues. The algorithm is based on the Conjugate-Gradient (CG) algorithm, however, it uses only a window of the vectors generated by the CG algorithm to compute approximate eigenvalues. The number and accuracy of the eigenvectors can be increased by solving more right-hand sides. For Hermitian systems with multiple right-hand sides, the computed eigenvectors can be used to speed up the solution of subsequent systems. The algorithm was tested on Lattice QCD problems by solving the normal equations and was shown to give large speed up factors and to remove the critical slowing down as we approach light quark masses. Here, an extension to the non-symmetric case based on the two-sided Lanczos algorithm is given. The new algorithm is tested on Lattice QCD problems and is shown to give promising results. We also study the removal of the critical slowing down and compare results with those of the eigCG algorithm. We also discuss the case when the system is gamma5-Hermitian.

This paper focuses on total least squares (TLS) problems $AX\approx B$ with multiple right-hand sides. Existence and uniqueness of a TLS solution for such problems was analyzed in the paper [I. Hnětynková et al., SIAM J. Matrix Anal. Appl., 32, 2011, pp. 748--770]. For TLS problems with single right-hand sides the paper [C. C. Paige and Z. Strakoš, SIAM J. Matrix Anal. Appl., 27, 2006, pp. 861--875] showed how necessary and sufficient information for solving $Ax\approx b$ can be revealed from the original data through the so-called core problem concept. In this paper we present a theoretical study extending this concept to problems with multiple right-hand sides. The data reduction we present here is based on the singular value decomposition of the system matrix $A$. We show minimality of the reduced problem; in this sense the situation is analogous to the single right-hand side case. Some other properties of the core problem, however, cannot be extended to the case of multiple right-hand sides.

A deflated and restarted Lanczos algorithm to solve hermitian linear systems, and at the same time compute eigenvalues and eigenvectors for application to multiple right-hand sides, is described.For the first right-hand side, eigenvectors with small eigenvalues are computed while simultaneously solving the linear system.Two versions of this algorithm are given.The first is called Lan-DR and is based on conjugate gradient (CG) implementation of the Lanczos algorithm.This version will be optimal for the hermitian positive definite case.The second version is called MinRes-DR and is based on the minimum residual (MinRes) implementation of Lanczos algorithm.This version is optimal for indefinite hermitian systems where the CG algorithm is subject to instabilities.For additional right-hand sides, we project over the calculated eigenvectors to speed up convergence.The algorithms used for subsequent right-hand sides are called D-CG and D-MinRes respectively.After some introductory examples are given, we show tests for the case of Wilson fermions at kappa critical.A considerable speed up in the convergence is observed compared to unmodified CG and MinRes.

Contemporary large-scale Partial Differential Equation (PDE) simulations usually require the solution of large and sparse linear systems. Moreover, it is often needed to solve these linear systems with different or multiple Right-Hand Sides (RHSs). In this paper, various strategies will be presented to extend the scalability of existing multigrid or domain decomposition linear solvers using appropriate recycling strategies or block methods—i.e., by treating multiple right-hand sides simultaneously. The scalability of this work is assessed by performing simulations on up to 8,192 cores for solving linear systems arising from various physical phenomena modeled by Poisson's equation, the system of linear elasticity, or Maxwell's equation. This work is shipped as part of on open-source software, readily available and usable in any C/C++, Python, or Fortran code. In particular, some simulations are performed on top of a well-established library, PETSc, and it is shown how our approaches can be used to decrease time to solution down by 30%.

The controlled Cholesky factorization (CCF) has been shown to be a robust preconditioner for the conjugate gradient (CG) method. In this scheme the amount of fill-in is defined in terms of a parameter $\eta$, the number of extra nonzero elements allowed per column. This parameter $\eta$ can be chosen in such a way that the preconditioning factor will require more or less storage than the coefficient matrix. The selection of nonzero elements kept for the preconditioner is made by value instead of by position, so the sparsity pattern of the original matrix is never considered by the proposed factorization. This preconditioned method is termed the CCCG($\eta$) method and is a generalization of a method developed by Jones and Plassmann [\textit{ACM Trans.\ Math.\ Software}, 21 (1995), pp.\ 5--17]. It is demonstrated how an optimum value of $\eta$ can be automatically determined when solving systems with multiple right-hand sides. Typically such systems occur when using time-stepping procedures. Numerical examples are presented for some time-dependent parabolic partial differential equations discretized using linear triangular finite elements on unstructured grids. Comparisons made between CCCG($\eta$), standard ICCG, and the Jones and Plassmann method show CCCG($\eta$) to be an efficient solver.

Recently, by combining the preconditioned accelerated Hermitian and skew-Hermitian splitting (PAHSS) and the parameterized Uzawa (PU) methods, Zheng and Ma (Appl Math Comput 273, 217–225, 2016b) presented the PAHSS–PU method for saddle point problems. By adding a block lower triangular matrix to the coefficient matrix on two sides of the first equation of the PAHSS–PU iterative scheme, the modified PAHSS–PU (MPAHSS–PU) iteration method is proposed in this paper, which has a faster convergence rate than the PAHSS–PU one. Furthermore, changing the position of the parameters in the MPAHSS–PU method, we develop another new method referred to as the modified PPHSS-SOR (MPPHSS-SOR) iteration method for solving saddle point problems. We provide the convergence properties of the MPAHSS–PU and the MPPHSS-SOR iteration methods, which show that the new methods are convergent if the related parameters satisfy suitable conditions. Meanwhile, practical ways to choose iteration parameters for the proposed methods are developed. Finally, numerical experiments demonstrate that the MPAHSS–PU and the MPPHSS-SOR methods outperform some existing ones both on the number of iterations and the computational times.

We present a seed method for solving large nonsymmetric linear systems with multiple right-hand sides. The method uses a single augmented Krylov subspace corresponding to a seed system as a generator of approximations to the nonseed systems. The residual evaluate of the method is shown, and a new strategy to form a seed system which could supply information shareable among the right-hand sides is given. Numerical experiments indicate that our seed selection strategy is more efficient than two existing strategies and our method has significant time saving compared with the block GMRES method and the GMRES method with a projection process.

Linear algebra operations play an important role in scientific computing and data analysis. With increasing data volume and complexity in the "Big Data" era, linear algebra operations are important tools to process massive datasets. On one hand, the advent of modern high-performance computing architectures with increasing computing power has greatly enhanced our capability to deal with a large volume of data. One the other hand, many classical, deterministic numerical linear algebra algorithms have difficulty to scale to handle large data sets. Monte Carlo methods, which are based on statistical sampling, exhibit many attractive properties in dealing with large volume of datasets, including fast approximated results, memory efficiency, reduced data accesses, natural parallelism, and inherent fault tolerance. In this dissertation, we present new Monte Carlo methods to accommodate a set of fundamental and ubiquitous large-scale linear algebra operations, including solving large-scale linear systems, constructing low-rank matrix approximation, and approximating the extreme eigenvalues/ eigenvectors, across modern distributed and parallel computing architectures. First of all, we revisit the classical Ulam-von Neumann Monte Carlo algorithm and derive the necessary and sufficient condition for its convergence. To support a broad family of linear systems, we develop Krylov subspace Monte Carlo solvers that go beyond the use of Neumann series. New algorithms used in the Krylov subspace Monte Carlo solvers include (1) a Breakdown-Free Block Conjugate Gradient algorithm to address the potential rank deficiency problem occurred in block Krylov subspace methods; (2) a Block Conjugate Gradient for Least Squares algorithm to stably approximate the least squares solutions of general linear systems; (3) a BCGLS algorithm with deflation to gain convergence acceleration; and (4) a Monte Carlo Generalized Minimal Residual algorithm based on sampling matrix-vector products to provide fast approximation of solutions. Secondly, we design a rank-revealing randomized Singular Value Decomposition (R3SVD) algorithm for adaptively constructing low-rank matrix approximations to satisfy application-specific accuracy. Thirdly, we study the block power method on Markov Chain Monte Carlo transition matrices and find that the convergence is actually depending on the number of independent vectors in the block. Correspondingly, we develop a sliding window power method to find stationary distribution, which has demonstrated success in modeling stochastic luminal Calcium release site. Fourthly, we take advantage of hybrid CPU-GPU computing platforms to accelerate the performance of the Breakdown-Free Block Conjugate Gradient algorithm and the randomized Singular Value Decomposition algorithm. Finally, we design a Gaussian variant of Freivalds’ algorithm to efficiently verify the correctness of matrix-matrix multiplication while avoiding undetectable fault patterns encountered in deterministic algorithms.

A new resolution method for solving linear systems resulting from the discrete adjoint approach in aerodynamic shape optimization with a single objective and multiple constraint functions is presented. The steady flow is governed by the three-dimensional compressible Favre-averaged Navier–Stokes equations combined with a two-equation turbulence model. Block flexible variants of the generalized minimal residual method are Krylov methods designed for the solution of linear systems with multiple right-hand sides allowing variable preconditioning. Then, they can be applied to the discrete adjoint equations of Favre-averaged Navier–Stokes equations for all objective and constraint functions as a right-hand side. The paper focuses on the development of a new algorithm combining the three following features: flexibility, multiple right-hand-side concurrent processing, and spectral deflation. This latter property enables to recycle approximate spectral information from one restart to the next one in order to avoid generalized minimal residual stagnation when short restart recurrences are performed. This novel block method is named block flexible generalized minimal residual with deflated restarting. It is evaluated on two aerodynamic configurations showing its capabilities to reduce computational times over standard Krylov methods.

Evaluating mixed-integer programming models over multiple right-hand sides

We present a fast direct solver for structured linear systems based on multilevel matrix compression. Using the recently developed interpolative decomposition of a low-rank matrix in a recursive manner, we embed an approximation of the original matrix into a larger, but highly structured sparse one that allows fast factorization and application of the inverse. The algorithm extends the Martinsson/Rokhlin method developed for 2D boundary integral equations and proceeds in two phases: a precomputation phase, consisting of matrix compression and factorization, followed by a solution phase to apply the matrix inverse. For boundary integral equations which are not too oscillatory, e.g., based on the Green's functions for the Laplace or low-frequency Helmholtz equations, both phases typically have complexity O(N) in two dimensions, where $N$ is the number of discretization points. In our current implementation, the corresponding costs in three dimensions are $O(N^{3/2})$ and $O(N \log N)$ for precomputation and solution, respectively. Extensive numerical experiments show a speedup of $\sim 100$ for the solution phase over modern fast multipole methods; however, the cost of precomputation remains high. Thus, the solver is particularly suited to problems where large numbers of iterations would be required. Such is the case with ill-conditioned linear systems or when the same system is to be solved with multiple right-hand sides. Our algorithm is implemented in Fortran and freely available.

Most calculations in lattice Quantum Chromodynamics (QCD) involve the solution of a series of linear systems of equations with exceedingly large matrices and a large number of right hand sides. Iterative methods for these problems can be sped up significantly if we deflate approximations of appropriate invariant spaces from the initial guesses. Recently we have developed eigCG, a modification of the Conjugate Gradient (CG) method, which while solving a linear system can reuse a window of the CG vectors to compute eigenvectors almost as accurately as the Lanczos method. The number of approximate eigenvectors can increase as more systems are solved. In this paper we review some of the characteristics of eigCG and show how it helps remove the critical slowdown in QCD calculations. Moreover, we study scaling with lattice volume and an extension of the technique to nonsymmetric problems.

In this paper, we introduce two new methods for solving large sparse nonsymmetric linear systems with several right-hand sides. These methods are the global Hessenberg and global CMRH methods. Using the global Hessenberg process, these methods are less expensive than the global FOM and global GMRES methods [9]. Theoretical results about the new methods are given, and experimental results that show good performances of these new methods are presented.