A New Iterative Solution Method for Solving Multiple Linear Systems

In this paper, we consider using conjugate gradient (CG) methods for solving multiple linear systems $A^{(i)} x^{(i)} = b^{(i)},$ for $1 \le i \le s,$ where the coefficient matrices $A^{(i)}$ and the right-hand sides $b^{(i)}$ are different in general. In particular, we focus on the seed projection method which generates a Krylov subspace from a set of direction vectors obtained by solving one of the systems, called the seed system, by the CG method and then projects the residuals of other systems onto the generated Krylov subspace to get the approximate solutions. The whole process is repeated until all the systems are solved. Most papers in the literature [T. F. Chan and W. L. Wan, {\it SIAM J. Sci. Comput.}, 18 (1997), 1698--1721; B. Parlett {\it Linear Algebra Appl.}, 29 (1980), 323--346; Y. Saad, {\it Math. Comp.}, 48 (1987), 651--662; V. Simoncini and E. Gallopoulos, {\it SIAM J. Sci. Comput.}, 16 (1995), 917--933; C. Smith, A. Peterson, and R. Mittra, {\it IEEE Trans. Antennas and Propagation}, 37 (1989), pp. 1490--1493] considered only the case where the coefficient matrices $A^{(i)}$ are the same but the right-hand sides are different. We extend and analyze the method to solve multiple linear systems with varying coefficient matrices and right-hand sides. A theoretical error bound is given for the approximation obtained from a projection process onto a Krylov subspace generated from solving a previous linear system. Finally, numerical results for multiple linear systems arising from image restorations and recursive least squares computations are reported to illustrate the effectiveness of the method.

We study the role of preconditioning strategies recently developed for coercive problems in connection with a two-step iterative method, based on the Hermitian skew-Hermitian splitting (HSS) of the coefficient matrix, proposed by Bai, Golub and Ng for the solution of nonsymmetric linear systems whose real part is coercive. As a model problem we consider Finite Differences (FD) matrix sequences {A n (a,p)} n discretizing the elliptic (convection-diffusion) problem with Ω being a plurirectangle of R d with a(x) being a uniformly positive function and p(x) denoting the Reynolds function: here for plurirectangle we mean a connected union of rectangles in d dimensions with edges parallel to the axes. More precisely, in connection with preconditioned HSS/GMRES like methods, we consider the preconditioning sequence {P n (a)} n , P n (a):= D n 1/2(a)A n (1,0) D n 1/2(a) where D n (a) is the suitably scaled main diagonal of A n (a,0). If a(x) is positive and regular enough, then the preconditioned sequence shows a strong clustering at unity so that the sequence {P n (a)} n turns out to be a superlinear preconditioning sequence for {A n (a,0)} n where A n (a,0) represents a good approximation of Re(A n (a,p)) namely the real part of A n (a,p). The computational interest is due to the fact that the preconditioned HSS method has a convergence behavior depending on the spectral properties of {P n -1(a)Re(A n (a,p))} n ≈ {P n -1(a)A n (a,0)} n : therefore the solution of a linear system with coefficient matrix A n (a,p) is reduced to computations involving diagonals and to the use of fast Poisson solvers for {A n (1,0)} n . Some numerical experimentations confirm the optimality of the discussed proposal and its superiority with respect to existing techniques.

This chapter discusses iterative methods for the solution of very large severely ill-conditioned linear systems of equations that arise from the discretization of linear ill-posed problems. The right-hand side vector represents the given data and is assumed to be contaminated by errors. Solution methods proposed in the literature employ some form of filtering to reduce the influence of the error in the right-hand side on the computed approximate solution. The amount of filtering is determined by a parameter often referred to as the regularization parameter. We discuss how the filtering affects the computed approximate solution and consider the selection of regularization parameter. Methods in which a suitable value of the regularization parameter is determined during the computation without user intervention are emphasized. New iterative solution methods based on expanding explicitly chosen filter functions in terms of Chebyshev polynomials are presented. The properties of these methods are illustrated with applications to image restoration.

Solution of a large sparse or dense linear systems often constitutes the most time consuming step when solving large multidimensional industrial problems on modern high performance supercomputers. In most practical cases these problems produce linear systems with ill-conditioned and off-diagonal dominant coefficient matrices. As the problem size grows, especially in the 3-dimensional case, direct methods for solving these systems (such as triangular decomposition) begin to lose their effectiveness due to the large arithmetic costs and memory requirements. These difficulties motivate resorting to iterative solution methods which are, at least theoretically, more efficient with respect to the arithmetic cost and memory requirements. On March 12, 1993 Qasim Sheikh from Cray Research presented numerical experiments with iterative solvers that are capable of solving large linear systems which arise in industrial applications. He then compared the effectiveness and some performance aspects of these iterative methods with other state of the art, highly optimized, iterative as well as direct solvers. The main results of his talk are based on a recent paper [1] written jointly with Kharchenko, Nikishin, Yeremin and Heroux.KeywordsIterative SolverPreconditioned MatrixLarge Linear SystemNonsymmetric Linear SystemTurbine PumpThese keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

Greedy low-rank approximation in Tucker format of solutions of tensor linear systems

In this paper, an iterative method is constructed to solve the following constrained linear matrix equations system: [A1(X),A2(X),... ,Ar(X)]=[E1,E2, ... ,Er ], X ∈ I={X |X= U(X)}, where Ai is a linear operator from Cmxn onto Cpixqi, Ei ∈ Cpixqi, i=1 , 2,..., r , and U is a linear self-conjugate involution operator. When the above constrained matrix equations system is consistent, for any initial matrix X0 ∈ I, a solution can be obtained by the proposed iterative method in finite iteration steps in the absence of roundoff errors, and the least Frobenius norm solution can be derived when a special kind of initial matrix is chosen. Furthermore, the optimal approximation solution to a given matrix can be derived. Several numerical examples are given to show the efficiency of the presented iterative method. Mathematical subject classification: 15A24, 65D99,65F30.

Dual-energy computed tomography (DECT) scans objects using two different X-ray spectra to acquire more information, which is also called dual spectral CT (DSCT) in some articles. Compared to traditional CT, DECT exhibits superior material distinguishability. Therefore, DECT can be widely used in the medical and industrial domains. However, owing to the nonlinearity and ill condition of DECT, studies are underway on DECT reconstruction to obtain high quality images and achieve fast convergence speed. Therefore, in this study, we propose an iterative reconstruction method based on monochromatic images (IRM-MI) to rapidly obtain high-quality images in DECT reconstruction. An IRM-MI is proposed for DECT. The proposed method converts DECT reconstruction problem from the basis material images decomposition to monochromatic images decomposition to significantly improve the convergence speed of DECT reconstruction by changing the coefficient matrix of the original equations to increase the angle of the high- and low-energy projection curves or reduce the condition number of the coefficient matrix. The monochromatic images were then decomposed into basis material images. Furthermore, we conducted numerical experiments to evaluate the performance of the proposed method. The decomposition results of the simulated data and real data experiments confirmed the effectiveness of the proposed method. Compared to the extended algebraic reconstruction technique (E-ART) method, the proposed method exhibited a significant increase in the convergence speed by increasing the angle of polychromatic projection curves or decreasing the condition number of the coefficient matrix, when choosing the appropriate monochromatic images. Therefore, the proposed method is also advantageous in acquiring high quality and rapidly converged images. We developed an iterative reconstruction method based on monochromatic images for the material decomposition for DECT. The numerical experiments using the proposed method validated its capability of decomposing the basis material images. Furthermore, the proposed method achieved faster convergence speed compared to the E-ART method.

ABSTRACTTriCG and TriMR are two short‐recurrence iterative methods proposed in recent years for solving symmetric quasi‐definite (SQD) linear systems. The basic mechanism underlying TriCG and TriMR is the generalized Saunders–Simon–Yip tridiagonalization process (gSSY). In this work, we show that TriCG and TriMR may fail to find the solution of SQD linear systems when gSSY terminates. To avoid this issue, we propose an improved generalized Saunders–Simon–Yip tridiagonalization process (igSSY). Based on igSSY, we introduce two new iterative methods named iTriCG and iTriMR for solving SQD linear systems in the same fashion as TriCG and TriMR. We show that iTriCG and iTriMR always terminate with the solution of SQD linear systems when igSSY terminates. In addition, we prove that the upper triangular factor of the QR factorization used in TriMR has only three nonzero diagonals, and based on this fact, we provide simplified short recurrences for TriMR, which reduce the work per iteration. These simplified short recurrences also apply to iTriMR. Numerical experiments are given to illustrate the merits of iTriCG and iTriMR.MSC2020 Classification: 15A06, 65F10, 65F25, 65F50

Efficient approximate solution of sparse linear systems

Developments and trends in the parallel solution of linear systems

This study presents an efficacy comparison of iterative solution methods for solving the compressible-fluid Reynolds equation in modeling air- or gas-lubricated bearings. A direct fixed-point iterative (DFI) method and Newton’s method are employed to transform the Reynolds equation in a form that can be solved iteratively. The iterative solution methods examined are the Gauss–Seidel method, the successive over-relaxation (SOR) method, the preconditioned conjugate gradient (PCG) method, and the multigrid method. The overall solution time is affected by both the transformation method and the iterative method applied. In this study, Newton’s method shows its effectiveness over the straightforward DFI method when the same iterative method is used. It is demonstrated that the use of an optimal relaxation factor is of vital importance for the efficiency of the SOR method. The multigrid method is an order faster than the PCG and optimal SOR methods. Also, the multigrid and PCG methods involve an extended coding work and are less flexible in dealing with gridwork and boundary conditions. Consequently, a compromise has to be made in terms of ease of use as well as programming effort for the solution of the compressible-fluid Reynolds equation.

The centered difference discretization of the spatial fractional coupled nonlinear Schrödinger equations obtains a discretized linear system whose coefficient matrix is the sum of a real diagonal matrix D and a complex symmetric Toeplitz matrix T̃ which is just the symmetric real Toeplitz T plus an imaginary identity matrix iI. In this study, we present a medium-shifted splitting iteration method to solve the discretized linear system, in which the fast algorithm can be utilized to solve the Toeplitz linear system. Theoretical analysis shows that the new iteration method is convergent. Moreover, the new splitting iteration method naturally leads to a preconditioner. Analysis shows that the eigenvalues of the corresponding preconditioned matrix are tighter than those of the original coefficient matrix A. Finally, compared with the other algorithms by numerical experiments, the new method is more effective.

In this dissertation we develop numerical algorithms and software tools to facilitate parallel solutions of nonlinear programming (NLP) problems. In particular, we address large-scale, block-structured problems with an intrinsic decomposable configuration. These problems arise in a great number of engineering applications, including parameter estimation, optimal control, network optimization, and stochastic programming. The structure of these problems can be leveraged by optimization solvers to accelerate solutions and overcome memory limitations, and we propose variants to two classes of optimization algorithms: augmented Lagrangian (AL) schemes and Schur-complement interior-point methods. The convergence properties of augmented Lagrangian decomposition schemes like the alternating direction method of multipliers (ADMM) and progressive hedging (PH) are well established for convex optimization but convergence guarantees in non-convex settings are still poorly understood. In practice, however, ADMM and PH often perform satisfactorily in complex non-convex NLPs. In this work, we study connections between the method of multipliers (MM), ADMM, and PH to derive benchmarking metrics that explain why PH and ADMM work in practice. We illustrate the concepts using challenging dynamic optimization problems. Our exposition seeks to establish more formalism in benchmarking ADMM, PH, and AL schemes and to motivate algorithmic improvements. The effectiveness of nonlinear interior-point solvers for solving large-scale problems relies quite heavily on the solution of the underlying linear algebra systems. The schur-complement decomposition is very effective for parallelizing the solution of linear systems with modest coupling. However, for systems with large number of coupling variables the schur-complement method does not scale favorably. We implement an approach that uses a Krylov solver (GMRES) preconditioned with ADMM to solve block-structured linear systems that arise in the interior-point method. We show that this ADMM-GMRES approach overcomes the well-known scalability issues of Schur decomposition. One important drawback of using decomposition approaches like ADMM and PH is their convergence rate. Unlike Schur-complement interior-point algorithms that have super-linear convergence, augmented Lagrangian approaches typically exhibit linear and sublinear rates. We exploit connections between ADMM and the Schur-complement decomposition to derive an accelerated version of ADMM. Specifically, we study the effectiveness of performing a Newton-Raphson algorithm to compute multiplier estimates for augmented Lagrangian methods. We demonstrate using two-stage stochastic programming problems that our multiplier update achieves convergence in fewer iterations for MM on general nonlinear problems. In the case of ADMM, the newton update significantly reduces the number of subproblem solves for convex quadratic programs (QPs). Moreover, we show that using newton multiplier updates makes the method robust to the selection of the penalty parameter. Traditionally, state-of-the-art optimization solvers are implemented in low-level programming languages. In our experience, the development of decomposition algorithms in these frameworks is challenging. They present a steep learning curve and can slow the development and testing of new numerical algorithms. To mitigate these challenges, we developed PyNumero, a new open source framework implemented in Python and C++. The package seeks to facilitate development of optimization algorithms for large-scale optimization within a high-level programming environment while at the same time minimizing the computational burden of using Python. The efficiency of PyNumero is illustrated by implementing algorithms for problems arising in stochastic programming and optimal control. Timing results are presented for both serial and parallel implementations. Our computational studies demonstrate that with the appropriate balance between compiled code and Python, efficient implementations of optimization algorithms are achievable in these high-level languages.

In literature, it has been reported that the convergence of some preconditioned stationary iterative methods using certain type upper triangular matrices as preconditioners are faster than the basic iterative methods. In this paper, a new preconditioned iterative method for the numerical solution of linear systems has been introduced, and the convergence analysis of the proposed method and an existing one have been done. Some numerical examples have also been given, which show the effectiveness of both of the methods.

This paper is a matrix iterative method presented to compute the solutions of the matrix equation, AXB=C, with unknown matrix X∈𝒮, where 𝒮 is the constrained matrices set like symmetric, symmetric-R-symmetric and (R, S)-symmetric. By this iterative method, for any initial matrix X 0∈𝒮, a solution X* can be obtained within finite iteration steps if exact arithmetics were used, and the solution X* with the minimum Frobenius norm can be obtained by choosing a special kind of initial matrix. The solution , which is nearest to a given matrix ˜X in Frobenius norm, can be obtained by first finding the minimum Frobenius norm solution of a new compatible matrix equation. The numerical examples given here show that the iterative method proposed in this paper has faster convergence and higher accuracy than the iterative methods proposed in [G.-X. Huang, F. Yin, and K. Guo, An iterative method for the skew-symmetric solution and the optimal approximate solution of the matrix equation AXB=C , J. Comput. Appl. Math. 212 (2008), pp. 231–244; Y. Lei and A.-P. Liao, A minimal residual algorithm for the inconsistent matrix equation AXB=C over symmetric matrices, Appl. Math. Comput. 188 (2007), pp. 499–513; Z.-Y. Peng, An iterative method for the least squares symmetric solution of the linear matrix equation AXB=C , Appl. Math. Comput. 170 (2005), pp. 711–723; Y.-X. Peng, X.-Y. Hu, and L. Zhang, An iteration method for the symmetric solutions and the optimal approximation solution of the matrix equation AXB=C , Appl. Math. Comput. 160 (2005), pp. 763–777].