For a non-Hermitian positive definite linear system Ax=bAx=b with Hermitian part HH and skew-Hermitian part SS, the Hermitian/skew-Hermitian splitting (HSS) iteration converges for every parameter α>0\alpha>0, and the same splitting matrix M(α)=12α(αI+H)(αI+S)M(\alpha)=\frac{1}{2\alpha}(\alpha I+H)(\alpha I+S) serves as the HSS preconditioner for Krylov subspace methods such as GMRES and BiCGSTAB. Chen (2015) proposed estimating α\alpha for the HSS iteration by equalizing the Euclidean condition numbers κ(αI+H)\kappa(\alpha I+H) and κ(αI+S)\kappa(\alpha I+S), which reduces to a cubic equation whose coefficients depend only on the extreme eigenvalues of HH and extreme singular values of SS. The source's closing remarks left open whether this estimate, tuned for the iteration, is effective when M(α)M(\alpha) is instead used as a preconditioner accelerating Krylov methods, and how to select the parameter when the extreme spectral quantities are too costly to compute and only matrix-vector products are available.

Background

Large sparse linear systems Ax=bAx=b with a non-Hermitian positive definite coefficient matrix (a matrix A=H+SA=H+S whose Hermitian part H=(A+A∗)/2H=(A+A^{*})/2 is positive definite and whose skew-Hermitian part S=(A−A∗)/2S=(A-A^{*})/2 is skew-Hermitian; every eigenvalue of such a matrix has positive real part, though the converse need not hold) arise from finite-difference and finite-element discretizations of convection-diffusion and related PDEs. The Hermitian and skew-Hermitian splitting (HSS) iteration of Bai, Golub and Ng (2003) writes A=H+SA=H+S and alternates two shifted half-steps

(αI+H)xk+1/2=(αI−S)xk+b,(αI+S)xk+1=(αI−H)xk+1/2+b,(\alpha I+H)x_{k+1/2}=(\alpha I-S)x_k+b,\qquad (\alpha I+S)x_{k+1}=(\alpha I-H)x_{k+1/2}+b,

where α>0\alpha>0 is a free parameter. The iteration converges for every α>0\alpha>0, but its speed is governed by the spectral radius ρ\rho (the largest modulus of the eigenvalues) of the iteration matrix, which depends critically on α\alpha. The classical convergence bound of Bai, Golub and Ng, δ(α)=max⁡λ∈λ(H)∣α−λα+λ∣\delta(\alpha)=\max_{\lambda\in\lambda(H)}|\frac{\alpha-\lambda}{\alpha+\lambda}|, depends only on the spectrum of HH and is minimized at α=λmin⁡λmax⁡\alpha=\sqrt{\lambda_{\min}\lambda_{\max}}, but it is often loose. Exact or quasi-optimal parameters are known only for special structures: Bai, Golub and Li (2006) solved certain two-by-two block matrices, Benzi, Gander and Golub (2003) optimized the HSS iteration for saddle-point problems, and Huang (2014) derived a practical formula that nevertheless performed poorly on strongly convective tests. Chen (2015) proposed instead to balance the conditioning of the two subproblems solved at each half-step: minimize φ(α)=∣κ(αI+H)−κ(αI+S)∣\varphi(\alpha)=|\kappa(\alpha I+H)-\kappa(\alpha I+S)|, where κ\kappa is the Euclidean (spectral) condition number. Setting φ(α)=0\varphi(\alpha)=0 reduces to the cubic

f(α)=2α3+(λmax⁡+λmin⁡−v)α2+2(σmin⁡2−λmin⁡v)α+σmin⁡2(λmax⁡+λmin⁡)−λmin⁡2v,f(\alpha)=2\alpha^3+(\lambda_{\max}+\lambda_{\min}-v)\alpha^2+2(\sigma_{\min}^2-\lambda_{\min}v)\alpha+\sigma_{\min}^2(\lambda_{\max}+\lambda_{\min})-\lambda_{\min}^2v,

with v=(σmax⁡2−σmin⁡2)/(λmax⁡−λmin⁡)v=(\sigma_{\max}^2-\sigma_{\min}^2)/(\lambda_{\max}-\lambda_{\min}), where λmin⁡,λmax⁡\lambda_{\min},\lambda_{\max} are the extreme eigenvalues of HH and σmin⁡,σmax⁡\sigma_{\min},\sigma_{\max} the extreme singular values of SS; a positive root gives the suggested parameter αsug\alpha_{sug}. On the standard two-dimensional convection-diffusion benchmark this parameter was competitive with the experimentally tuned optimum and far better than Huang's estimate. Notably, the source itself establishes no theoretical link between condition-number balancing and the spectral radius of the iteration matrix.

The same splitting matrix M(α)=12α(αI+H)(αI+S)M(\alpha)=\frac{1}{2\alpha}(\alpha I+H)(\alpha I+S) can instead be used as a preconditioner (an approximate inverse applied to the system so that the iterative method converges faster) for a Krylov subspace method — an iterative solver whose approximations are drawn from the subspaces spanned by successive applications of the coefficient matrix to the initial residual — such as GMRES (the generalized minimal residual method, which minimizes the residual norm over these growing subspaces; restarted variants cap the subspace size to limit memory) or BiCGSTAB (a stabilized biconjugate-gradient variant) — applied to the preconditioned system M(α)−1Ax=M(α)−1bM(\alpha)^{-1}Ax=M(\alpha)^{-1}b, and this is frequently more effective than the plain iteration; Bertaccini, Golub, Serra-Capizzano and Tablino-Possio (2005) studied such preconditioned HSS-type methods for discrete convection-diffusion equations. A structural identity links the two uses: writing the splitting as A=M(α)−N(α)A=M(\alpha)-N(\alpha), the eigenvalues of the preconditioned matrix satisfy λ(M(α)−1A)=1−λ(M(α)−1N(α))\lambda(M(\alpha)^{-1}A)=1-\lambda(M(\alpha)^{-1}N(\alpha)), so the parameter does steer the preconditioned spectrum. However, convergence of GMRES and BiCGSTAB depends on the eigenvectors and the non-normality (departure from orthogonality) of the preconditioned operator, not on eigenvalues alone, so a parameter that is good for the iteration — or for an eigenvalue enclosure — need not minimize the Krylov iteration count. Spectral characterizations of HSS preconditioners exist only for special classes, chiefly saddle-point matrices (Simoncini and Benzi 2004; Bai 2018). For the general non-Hermitian positive definite class no theory identifies the parameter of the HSS preconditioner that minimizes, or provably bounds, the cost of preconditioned GMRES or BiCGSTAB, and practitioners tune α\alpha experimentally. Moreover, all a priori estimates, including the cubic above, require the extreme eigenvalues of HH and the extreme singular values of SS, which are themselves expensive to compute for very large systems; in matrix-free settings, where the solver can only apply AA and A∗A^{*} to vectors, a parameter-estimation procedure with characterized accuracy is missing.

Problem Statement

For a non-Hermitian positive definite matrix A=H+SA=H+S (HH Hermitian positive definite, SS skew-Hermitian) and the HSS splitting matrix M(α)=12α(αI+H)(αI+S)M(\alpha)=\frac{1}{2\alpha}(\alpha I+H)(\alpha I+S), two coupled questions arising from Chen (2015) remain unanswered. Both parts are required, and each is resolved on its own terms: an affirmative or negative answer to part (1) fully resolves that part, while part (2) is resolved by a matrix-free estimation procedure together with a characterization of its accuracy.

(1) Transfer to preconditioning. Does the condition-number-balancing estimate — the positive root of the cubic obtained from κ(αI+H)=κ(αI+S)\kappa(\alpha I+H)=\kappa(\alpha I+S) with coefficients in λmin⁡(H)\lambda_{\min}(H), λmax⁡(H)\lambda_{\max}(H), σmin⁡(S)\sigma_{\min}(S), σmax⁡(S)\sigma_{\max}(S) — also select a good parameter when M(α)M(\alpha) is used as a preconditioner for GMRES or BiCGSTAB applied to M(α)−1Ax=M(α)−1bM(\alpha)^{-1}Ax=M(\alpha)^{-1}b? Here "good" is measured by the preconditioned iteration count of GMRES or BiCGSTAB (or a provable bound on it) over a stated matrix class, compared with the best empirically tuned parameter for that class.

(2) Matrix-free parameter estimation. How can the parameter of the HSS preconditioner (and of the HSS iteration) be estimated using only matrix-vector products with AA and A∗A^{*}, i.e., without computing the extreme eigenvalues of HH or the extreme singular values of SS explicitly?

Both parts serve one scientific objective: principled, practically computable parameter selection that makes HSS-based iteration and HSS-preconditioned Krylov solution of non-Hermitian positive definite systems reliable without case-by-case experimental tuning.

References

  1. Fang Chen (2015), On choices of iteration parameter in HSS method, Applied Mathematics and Computation 271:832–837, DOI 10.1016/j.amc.2015.09.003, https://doi.org/10.1016/j.amc.2015.09.003 (source paper)
  2. Zhong-Zhi Bai, Gene H. Golub, Michael K. Ng (2003), Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems, SIAM Journal on Matrix Analysis and Applications 24:603–626, DOI 10.1137/S0895479801395458, https://doi.org/10.1137/S0895479801395458
  3. Zhong-Zhi Bai, Gene H. Golub, Chi-Kwong Li (2006), Optimal parameter in Hermitian and skew-Hermitian splitting method for certain two-by-two block matrices, SIAM Journal on Scientific Computing 28:583–603, DOI 10.1137/050623644, https://doi.org/10.1137/050623644
  4. Yu-Mei Huang (2014), A practical formula for computing optimal parameters in the HSS iteration methods, Journal of Computational and Applied Mathematics 255:142–149, DOI 10.1016/j.cam.2013.01.023, https://doi.org/10.1016/j.cam.2013.01.023
  5. Ai-Li Yang (2018), Scaled norm minimization method for computing the parameters of the HSS and the two-parameter HSS preconditioners, Numerical Linear Algebra with Applications 25:e2169, DOI 10.1002/nla.2169, https://doi.org/10.1002/nla.2169
  6. Yuebin Hao, Ai-Li Yang, Yu-Jiang Wu (2018), How to compute the optimal parameters of the parameterized PMHSS preconditioner?, Applied Mathematics Letters 86:8–13, DOI 10.1016/j.aml.2018.06.009, https://doi.org/10.1016/j.aml.2018.06.009
  7. Qinmeng Zou, Frédéric Magoulès (2020), Parameter estimation in the Hermitian and skew-Hermitian splitting method using gradient iterations, Numerical Linear Algebra with Applications 27:e2304, DOI 10.1002/nla.2304, https://doi.org/10.1002/nla.2304
  8. Na Huang (2022), Variable-parameter HSS methods for non-Hermitian positive definite linear systems, Linear and Multilinear Algebra 70:6664–6681, DOI 10.1080/03081087.2021.1968328, https://doi.org/10.1080/03081087.2021.1968328
  9. Mehdi Dehghan, Akbar Shirilord (2022), Approximating optimal parameters for generalized preconditioned Hermitian and skew-Hermitian splitting (GPHSS) method, Computational and Applied Mathematics 41, DOI 10.1007/s40314-021-01747-5, https://doi.org/10.1007/s40314-021-01747-5
  10. Fang Chen, Tian-Yi Li, Galina V. Muratova (2021), Lopsided scaled HSS preconditioner for steady-state space-fractional diffusion equations, Calcolo 58, DOI 10.1007/s10092-021-00419-4, https://doi.org/10.1007/s10092-021-00419-4
  11. Tian-Yi Li, Fang Chen, Zhi-Wei Fang, Hai-Wei Sun, Zhi Wang (2024), Two-parameter modified matrix splitting iteration method for Helmholtz equation, International Journal of Computer Mathematics 101:1205–1218, DOI 10.1080/00207160.2023.2301570, https://doi.org/10.1080/00207160.2023.2301570
  12. Valeria Simoncini, Michele Benzi (2004), Spectral properties of the Hermitian and skew-Hermitian splitting preconditioner for saddle point problems, SIAM Journal on Matrix Analysis and Applications 26:377–389, DOI 10.1137/S0895479803434926, https://doi.org/10.1137/S0895479803434926
  13. Daniele Bertaccini, Gene H. Golub, Stefano Serra-Capizzano, Cristina Tablino-Possio (2005), Preconditioned HSS methods for the solution of non-Hermitian positive definite linear systems and applications to the discrete convection-diffusion equation, Numerische Mathematik 99:441–484, DOI 10.1007/s00211-004-0574-1, https://doi.org/10.1007/s00211-004-0574-1
  14. Zhong-Zhi Bai (2018), On spectral clustering of HSS preconditioner for generalized saddle-point matrices, Linear Algebra and its Applications 555:285–300, DOI 10.1016/j.laa.2018.06.022, https://doi.org/10.1016/j.laa.2018.06.022
  15. Michele Benzi (2009), A generalization of the Hermitian and skew-Hermitian splitting iteration, SIAM Journal on Matrix Analysis and Applications 31:360–374, DOI 10.1137/080723181, https://doi.org/10.1137/080723181
  16. Michele Benzi, Martin J. Gander, Gene H. Golub (2003), Optimization of the Hermitian and skew-Hermitian splitting iteration for saddle-point problems, BIT Numerical Mathematics 43:881–900, DOI 10.1023/B:BITN.0000014548.26616.65, https://doi.org/10.1023/B:BITN.0000014548.26616.65