Source author record

Zhongxiao Jia

Zhongxiao Jia appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

20works
2topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

20 published item(s)

preprint2021arXiv

A comparison of eigenvalue-based algorithms and the generalized Lanczos trust-region algorithm for Solving the trust-region subproblem

Solving the trust-region subproblem (TRS) plays a key role in numerical optimization and many other applications. Based on a fundamental result that the solution of TRS of size $n$ is mathematically equivalent to finding the rightmost eigenpair of a certain matrix pair of size $2n$, eigenvalue-based methods are promising due to their simplicity. For $n$ large, the implicitly restarted Arnoldi (IRA) and refined Arnoldi (IRRA) algorithms are well suited for this eigenproblem. For a reasonable comparison of overall efficiency of the algorithms for solving TRS directly and eigenvalue-based algorithms, a vital premise is that the two kinds of algorithms must compute the approximate solutions of TRS with (almost) the same accuracy, but such premise has been ignored in the literature. To this end, we establish close relationships between the two kinds of residual norms, so that, given a stopping tolerance for IRA and IRRA, we are able to determine a reliable one that GLTR should use so as to ensure that GLTR and IRA, IRRA deliver the converged approximate solutions with similar accuracy. We also make a convergence analysis on the residual norms by the Generalized Lanczos Trust-Region (GLTR) algorithm for solving TRS directly, the Arnoldi method and the refined Arnoldi method for the equivalent eigenproblem. A number of numerical experiments are reported to illustrate that IRA and IRRA are competitive with GLTR and IRRA outperforms IRA.

preprint2020arXiv

The Krylov Subspaces, Low Rank Approximations and Ritz Values of LSQR for Linear Discrete Ill-Posed Problems: the Multiple Singular Value Case

For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or $Ax=b$ with $b$ contaminated by white noise, the Golub-Kahan bidiagonalization based LSQR method and its mathematically equivalent CGLS, the Conjugate Gradient (CG) method applied to $A^TAx=A^Tb$, are most commonly used. They have intrinsic regularizing effects, where the iteration number $k$ plays the role of regularization parameter. The long-standing fundamental question is: {\em Can LSQR and CGLS find 2-norm filtering best possible regularized solutions}? The author has given definitive answers to this question for severely and moderately ill-posed problems when the singular values of $A$ are simple. This paper extends the results to the multiple singular value case, and studies the approximation accuracy of Krylov subspaces, the quality of low rank approximations generated by Golub-Kahan bidiagonalization and the convergence properties of Ritz values. For the two kinds of problems, we prove that LSQR finds 2-norm filtering best possible regularized solutions at semi-convergence. Particularly, we consider some important and untouched issues on best, near best and general rank $k$ approximations to $A$ for the ill-posed problems with the singular values $σ_k=\mathcal{O}(k^{-α})$ with $α>0$, and the relationships between them and their nonzero singular values. Numerical experiments confirm our theory. The results on general rank $k$ approximations and the properties of their nonzero singular values apply to several Krylov solvers, including LSQR, CGME, MINRES, MR-II, GMRES and RRGMRES.

preprint2019arXiv

Approximation Accuracy of the Krylov Subspaces for Linear Discrete Ill-Posed Problems

For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or $Ax=b$ with $b$ contaminated by Gaussian white noise, the Lanczos bidiagonalization based Krylov solver LSQR and its mathematically equivalent CGLS, the Conjugate Gradient (CG) method implicitly applied to $A^TAx=A^Tb$, are most commonly used, and CGME, the CG method applied to $\min\|AA^Ty-b\|$ or $AA^Ty=b$ with $x=A^Ty$, and LSMR, which is equivalent to the minimal residual (MINRES) method applied to $A^TAx=A^Tb$, have also been choices. These methods exhibit typical semi-convergence feature, and the iteration number $k$ plays the role of the regularization parameter. However, there has been no definitive answer to the long-standing fundamental question: {\em Can LSQR and CGLS find 2-norm filtering best possible regularized solutions}? The same question is for CGME and LSMR too. At iteration $k$, LSQR, CGME and LSMR compute {\em different} iterates from the {\em same} $k$ dimensional Krylov subspace. A first and fundamental step towards to answering the above question is to {\em accurately} estimate the accuracy of the underlying $k$ dimensional Krylov subspace approximating the $k$ dimensional dominant right singular subspace of $A$. Assuming that the singular values of $A$ are simple, we present a general $\sinΘ$ theorem for the 2-norm distances between these two subspaces and derive accurate estimates on them for severely, moderately and mildly ill-posed problems. We also establish some relationships between the smallest Ritz values and these distances. Numerical experiments justify the sharpness of our results.

preprint2019arXiv

Regularization Properties of the Krylov Iterative Solvers CGME and LSMR For Linear Discrete Ill-Posed Problems with an Application to Truncated Randomized SVDs

For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or $Ax=b$ with $b$ contaminated by Gaussian white noise, there are four commonly used Krylov solvers: LSQR and its mathematically equivalent CGLS, the Conjugate Gradient (CG) method applied to $A^TAx=A^Tb$, CGME, the CG method applied to $\min\|AA^Ty-b\|$ or $AA^Ty=b$ with $x=A^Ty$, and LSMR, the minimal residual (MINRES) method applied to $A^TAx=A^Tb$. These methods have intrinsic regularizing effects, where the number $k$ of iterations plays the role of the regularization parameter. In this paper, we establish a number of regularization properties of CGME and LSMR, including the filtered SVD expansion of CGME iterates, and prove that the 2-norm filtering best regularized solutions by CGME and LSMR are less accurate than and at least as accurate as those by LSQR, respectively. We also prove that the semi-convergence of CGME and LSMR always occurs no later and sooner than that of LSQR, respectively. As a byproduct, using the analysis approach for CGME, we improve a fundamental result on the accuracy of the truncated rank $k$ approximate SVD of $A$ generated by randomized algorithms, and reveal how the truncation step damages the accuracy. Numerical experiments justify our results on CGME and LSMR.

preprint2019arXiv

The Low Rank Approximations and Ritz Values in LSQR For Linear Discrete Ill-Posed Problems

LSQR and its mathematically equivalent CGLS have been popularly used over the decades for large-scale linear discrete ill-posed problems, where the iteration number $k$ plays the role of the regularization parameter. It has been long known that if the Ritz values in LSQR converge to the large singular values of $A$ in natural order until its semi-convergence then LSQR must have the same the regularization ability as the truncated singular value decomposition (TSVD) method and can compute a 2-norm filtering best possible regularized solution. However, hitherto there has been no definitive rigorous result on the approximation behavior of the Ritz values in the context of ill-posed problems. In this paper, for severely, moderately and mildly ill-posed problems, we give accurate solutions of the two closely related fundamental and highly challenging problems on the regularization of LSQR: (i) How accurate are the low rank approximations generated by Lanczos bidiagonalization? (ii) Whether or not the Ritz values involved in LSQR approximate the large singular values of $A$ in natural order? We also show how to judge the accuracy of low rank approximations reliably during computation without extra cost. Numerical experiments confirm our results.

preprint2018arXiv

A Joint Bidiagonalization Based Algorithm for Large Scale Linear Discrete Ill-posed Problems in General-Form Regularization

Based on the joint bidiagonalization process of a large matrix pair $\{A,L\}$, we propose and develop an iterative regularization algorithm for the large scale linear discrete ill-posed problems in general-form regularization: $\min\|Lx\| \ \mbox{\rm subject to} \ x\in\mathcal{S} = \{x|\ \|Ax-b\|\leq τ\|e\|\}$ with a Gaussian white noise $e$ and $τ>1$ slightly, where $L$ is a regularization matrix. Our algorithm is different from the hybrid one proposed by Kilmer {\em et al.}, which is based on the same process but solves the general-form Tikhonov regularization problem: $\min_x\left\{\|Ax-b\|^2+λ^2\|Lx\|^2\right\}$. We prove that the iterates take the form of attractive filtered generalized singular value decomposition (GSVD) expansions, where the filters are given explicitly. This result and the analysis on it show that the method must have the desired semi-convergence property and get insight into the regularizing effects of the method. We use the L-curve criterion or the discrepancy principle to determine $k^*$. The algorithm is simple and effective, and numerical experiments illustrate that it often computes more accurate regularized solutions than the hybrid one.

preprint2016arXiv

The Regularization Theory of the Krylov Iterative Solvers LSQR, CGLS, LSMR and CGME For Linear Discrete Ill-Posed Problems

For the large-scale linear discrete ill-posed problem $\min\|Ax-b\|$ or $Ax=b$ with $b$ contaminated by a white noise, Lanczos bidiagonalization based LSQR and its mathematically equivalent CGLS are most commonly used. They have intrinsic regularizing effects, where the number $k$ of iterations plays the role of regularization parameter. However, hitherto there has been no answer to the long-standing fundamental concern of Björck and Eldén in 1979: {\em for which kinds of problems LSQR and CGLS can find best possible regularized solutions}? Here a best possible regularized solution means that it is at least as accurate as the best regularized solution obtained by the truncated singular value decomposition (TSVD) method or by standard-form Tikhonov regularization and cannot be improved under certain conditions. In this paper we make a detailed analysis on the regularization of LSQR for severely, moderately and mildly ill-posed problems. For the first two kinds of problems, we prove that LSQR finds best possible solutions at semi-convergence and the following results hold until semi-convergence: (i) the $k$-step Lanczos bidiagonalization always generates a near best rank $k$ approximation to $A$; (ii) the $k$ Ritz values always approximate the first $k$ large singular values of $A$ in natural order; (iii) the $k$-step LSQR always captures the $k$ dominant SVD components of $A$; (iv) the diagonals and subdiagonals of the bidiagonal matrices generated by Lanczos bidiagonalization decay as fast as the singular values of $A$. However, for the third kind of problem, the above results do not hold generally. We also analyze the regularization of the other two Krylov solvers LSMR and CGME, proving that LSMR has similar regularizing effects to LSQR for each kind of problem and both are superior to CGME. Numerical experiments confirm our theory on LSQR.

preprint2015arXiv

An Approach to Making SPAI and PSAI Preconditioning Effective for Large Irregular Sparse Linear Systems

We investigate the SPAI and PSAI preconditioning procedures and shed light on two important features of them: (i) For the large linear system $Ax=b$ with $A$ irregular sparse, i.e., with $A$ having $s$ relatively dense columns, SPAI may be very costly to implement, and the resulting sparse approximate inverses may be ineffective for preconditioning. PSAI can be effective for preconditioning but may require excessive storage and be unacceptably time consuming; (ii) the situation is improved drastically when $A$ is regular sparse, that is, all of its columns are sparse. In this case, both SPAI and PSAI are efficient. Moreover, SPAI and, especially, PSAI are more likely to construct effective preconditioners. Motivated by these features, we propose an approach to making SPAI and PSAI more practical for $Ax=b$ with $A$ irregular sparse. We first split $A$ into a regular sparse $\tilde A$ and a matrix of low rank $s$. Then exploiting the Sherman--Morrison--Woodbury formula, we transform $Ax=b$ into $s+1$ new linear systems with the same coefficient matrix $\tilde A$, use SPAI and PSAI to compute sparse approximate inverses of $\tilde A$ efficiently and apply Krylov iterative methods to solve the preconditioned linear systems. Theoretically, we consider the non-singularity and conditioning of $\tilde A$ obtained from some important classes of matrices. We show how to recover an approximate solution of $Ax=b$ from those of the $s+1$ new systems and how to design reliable stopping criteria for the $s+1$ systems to guarantee that the approximate solution of $Ax=b$ satisfies a desired accuracy. Given the fact that irregular sparse linear systems are common in applications, this approach widely extends the practicability of SPAI and PSAI. Numerical results demonstrate the considerable superiority of our approach to the direct application of SPAI and PSAI to $Ax=b$.

preprint2015arXiv

Harmonic and Refined Harmonic Shift-Invert Residual Arnoldi and Jacobi--Davidson Methods for Interior Eigenvalue Problems

This paper concerns the harmonic shift-invert residual Arnoldi (HSIRA) and Jacobi--Davidson (HJD) methods as well as their refined variants RHSIRA and RHJD for the interior eigenvalue problem. Each method needs to solve an inner linear system to expand the subspace successively. When the linear systems are solved only approximately, we are led to the inexact methods. We prove that the inexact HSIRA, RHSIRA, HJD and RHJD methods mimic their exact counterparts well when the inner linear systems are solved with only low or modest accuracy. We show that (i) the exact HSIRA and HJD expand subspaces better than the exact SIRA and JD and (ii) the exact RHSIRA and RHJD expand subspaces better than the exact HSIRA and HJD. Based on the theory, we design stopping criteria for inner solves. To be practical, we present restarted HSIRA, HJD, RHSIRA and RHJD algorithms. Numerical results demonstrate that these algorithms are much more efficient than the restarted standard SIRA and JD algorithms and furthermore the refined harmonic algorithms outperform the harmonic ones very substantially.

preprint2015arXiv

Implicitly Restarted Generalized Second-order Arnoldi Type Algorithms for the Quadratic Eigenvalue Problem

We investigate the generalized second-order Arnoldi (GSOAR) method, a generalization of the SOAR method proposed by Bai and Su [{\em SIAM J. Matrix Anal. Appl.}, 26 (2005): 640--659.], and the Refined GSOAR (RGSOAR) method for the quadratic eigenvalue problem (QEP). The two methods use the GSOAR procedure to generate an orthonormal basis of a given generalized second-order Krylov subspace, and with such basis they project the QEP onto the subspace and compute the Ritz pairs and the refined Ritz pairs, respectively. We develop implicitly restarted GSOAR and RGSOAR algorithms, in which we propose certain exact and refined shifts for respective use within the two algorithms. Numerical experiments on real-world problems illustrate the efficiency of the restarted algorithms and the superiority of the restarted RGSOAR to the restarted GSOAR. The experiments also demonstrate that both IGSOAR and IRGSOAR generally perform much better than the implicitly restarted Arnoldi method applied to the corresponding linearization problems, in terms of the accuracy and the computational efficiency.

preprint2015arXiv

On Convergence of the Inexact Rayleigh Quotient Iteration with the Lanczos Method Used for Solving Linear Systems

For the Hermitian inexact Rayleigh quotient iteration (RQI), the author has established new local general convergence results, independent of iterative solvers for inner linear systems. The theory shows that the method locally converges quadratically under a new condition, called the uniform positiveness condition. In this paper we first consider the local convergence of the inexact RQI with the unpreconditioned Lanczos method for the linear systems. Some attractive properties are derived for the residuals, whose norms are $ξ_{k+1}$'s, of the linear systems obtained by the Lanczos method. Based on them and the new general convergence results, we make a refined analysis and establish new local convergence results. It is proved that the inexact RQI with Lanczos converges quadratically provided that $ξ_{k+1}\leqξ$ with a constant $ξ\geq 1$. The method is guaranteed to converge linearly provided that $ξ_{k+1}$ is bounded by a small multiple of the reciprocal of the residual norm $\|r_k\|$ of the current approximate eigenpair. The results are fundamentally different from the existing convergence results that always require $ξ_{k+1}<1$, and they have a strong impact on effective implementations of the method. We extend the new theory to the inexact RQI with a tuned preconditioned Lanczos for the linear systems. Based on the new theory, we can design practical criteria to control $ξ_{k+1}$ to achieve quadratic convergence and implement the method more effectively than ever before. Numerical experiments confirm our theory.

preprint2015arXiv

On Inner Iterations in the Shift-Invert Residual Arnoldi Method and the Jacobi--Davidson Method

Using a new analysis approach, we establish a general convergence theory of the Shift-Invert Residual Arnoldi (SIRA) method for computing a simple eigenvalue nearest to a given target $σ$ and the associated eigenvector. In SIRA, a subspace expansion vector at each step is obtained by solving a certain inner linear system. We prove that the inexact SIRA method mimics the exact SIRA well, that is, the former uses almost the same outer iterations to achieve the convergence as the latter does if all the inner linear systems are iteratively solved with {\em low} or {\em modest} accuracy during outer iterations. Based on the theory, we design practical stopping criteria for inner solves. Our analysis is on one step expansion of subspace and the approach applies to the Jacobi--Davidson (JD) method with the fixed target $σ$ as well, and a similar general convergence theory is obtained for it. Numerical experiments confirm our theory and demonstrate that the inexact SIRA and JD are similarly effective and are considerably superior to the inexact SIA.

preprint2015arXiv

On the Condition Number of the Total Least Squares Problem

This paper concerns singular value decomposition (SVD)-based computable formulas and bounds for the condition number of the Total Least Squares (TLS) problem. For the TLS problem with the coefficient matrix $A$ and the right-hand side $b$, a new closed formula is presented for the condition number. Unlike an important result in the literature that uses the SVDs of both $A$ and $[A,\ b]$, our formula only requires the SVD of $[A,\ b]$. Based on the closed formula, both lower and upper bounds for the condition number are derived. It is proved that they are always sharp and estimate the condition number accurately. A few lower and upper bounds are further established that involve at most the smallest two singular values of $A$ and of $[A,\ b]$. Tightness of these bounds is discussed, and numerical experiments are presented to confirm our theory and to demonstrate the improvement of our upper bounds over the two upper bounds due to Golub and Van Loan as well as Baboulin and Gratton. Such lower and upper bounds are particularly useful for large scale TLS problems since they require the computation of only a few singular values of $A$ and $[A, \ b]$ other than all the singular values of them.

preprint2015arXiv

On the Convergence of Ritz Pairs and Refined Ritz Vectors for Quadratic Eigenvalue Problems

For a given subspace, the Rayleigh-Ritz method projects the large quadratic eigenvalue problem (QEP) onto it and produces a small sized dense QEP. Similar to the Rayleigh-Ritz method for the linear eigenvalue problem, the Rayleigh-Ritz method defines the Ritz values and the Ritz vectors of the QEP with respect to the projection subspace. We analyze the convergence of the method when the angle between the subspace and the desired eigenvector converges to zero. We prove that there is a Ritz value that converges to the desired eigenvalue unconditionally but the Ritz vector converges conditionally and may fail to converge. To remedy the drawback of possible non-convergence of the Ritz vector, we propose a refined Ritz vector that is mathematically different from the Ritz vector and is proved to converge unconditionally. We construct examples to illustrate our theory.

preprint2015arXiv

Robust Dropping Criteria for F-norm Minimization Based Sparse Approximate Inverse Preconditioning

Dropping tolerance criteria play a central role in Sparse Approximate Inverse preconditioning. Such criteria have received, however, little attention and have been treated heuristically in the following manner: If the size of an entry is below some empirically small positive quantity, then it is set to zero. The meaning of "small" is vague and has not been considered rigorously. It has not been clear how dropping tolerances affect the quality and effectiveness of a preconditioner $M$. In this paper, we focus on the adaptive Power Sparse Approximate Inverse algorithm and establish a mathematical theory on robust selection criteria for dropping tolerances. Using the theory, we derive an adaptive dropping criterion that is used to drop entries of small magnitude dynamically during the setup process of $M$. The proposed criterion enables us to make $M$ both as sparse as possible as well as to be of comparable quality to the potentially denser matrix which is obtained without dropping. As a byproduct, the theory applies to static F-norm minimization based preconditioning procedures, and a similar dropping criterion is given that can be used to sparsify a matrix after it has been computed by a static sparse approximate inverse procedure. In contrast to the adaptive procedure, dropping in the static procedure does not reduce the setup time of the matrix but makes the application of the sparser $M$ for Krylov iterations cheaper. Numerical experiments reported confirm the theory and illustrate the robustness and effectiveness of the dropping criteria.

preprint2015arXiv

Some Results on Regularization of LSQR and CGLS for Large-Scale Discrete Ill-Posed Problems

For large-scale discrete ill-posed problems, LSQR, a Lanczos bidiagonalization process based Krylov method, is most often used. It is well known that LSQR has natural regularizing properties, where the number of iterations plays the role of the regularization parameter. In this paper, for severely and moderately ill-posed problems, we establish quantitative bounds for the distance between the $k$-dimensional Krylov subspace and the subspace spanned by $k$ dominant right singular vectors. They show that the $k$-dimensional Krylov subspace may capture the $k$ dominant right singular vectors for severely and moderately ill-posed problems, but it seems not the case for mildly ill-posed problems. These results should be the first step towards to estimating the accuracy of the rank-$k$ approximation generated by Lanczos bidiagonalization. We also derive some other results, which help further understand the regularization effects of LSQR. We draw to a conclusion that a hybrid LSQR should generally be used for mildly ill-posed problems. We report numerical experiments to confirm our theory. We present more definitive and general observed phenomena, which will derive more research.

preprint2014arXiv

A Posteriori Error Estimates of Krylov Subspace Approximations to Matrix Functions

Krylov subspace methods for approximating a matrix function $f(A)$ times a vector $v$ are analyzed in this paper. For the Arnoldi approximation to $e^{-τA}v$, two reliable a posteriori error estimates are derived from the new bounds and generalized error expansion we establish. One of them is similar to the residual norm of an approximate solution of the linear system, and the other one is determined critically by the first term of the error expansion of the Arnoldi approximation to $e^{-τA}v$ due to Saad. We prove that each of the two estimates is reliable to measure the true error norm, and the second one theoretically justifies an empirical claim by Saad. In the paper, by introducing certain functions $ϕ_k(z)$ defined recursively by the given function $f(z)$ for certain nodes, we obtain the error expansion of the Krylov-like approximation for $f(z)$ sufficiently smooth, which generalizes Saad's result on the Arnoldi approximation to $e^{-τA}v$. Similarly, it is shown that the first term of the generalized error expansion can be used as a reliable a posteriori estimate for the Krylov-like approximation to some other matrix functions times $v$. Numerical examples are reported to demonstrate the effectiveness of the a posteriori error estimates for the Krylov-like approximations to $e^{-τA}v$, $\cos(A)v$ and $\sin(A)v$.

preprint2011arXiv

A contribution to the condition number of the total least squares problem

This paper concerns cheaply computable formulas and bounds for the condition number of the TLS problem. For a TLS problem with data $A$, $b$, two formulas are derived that are simpler and more compact than the known results in the literature. One is derived by exploiting the properties of Kronecker products of matrices. The other is obtained by making use of the singular value decomposition (SVD) of $[A \,\,b]$, which allows us to compute the condition number cheaply and accurately. We present lower and upper bounds for the condition number that involve the singular values of $[A \,\, b]$ and the last entries of the right singular vectors of $[A \,\, b]$. We prove that they are always sharp and can estimate the condition number accurately by no more than four times. Furthermore, we establish a few other lower and upper bounds that involve only a few singular values of $A$ and $[A \,\, b]$. We discuss how tight the bounds are. These bounds are particularly useful for large scale TLS problems since for them any formulas and bounds for the condition number involving all the singular values of $A$ and/or $[A \ b]$ are too costly to be computed. Numerical experiments illustrate that our bounds are sharper than a known approximate condition number in the literature.

preprint2010arXiv

On Convergence of the Inexact Rayleigh Quotient Iteration with MINRES

For the Hermitian inexact Rayleigh quotient iteration (RQI), we present a new general theory, independent of iterative solvers for shifted inner linear systems. The theory shows that the method converges at least quadratically under a new condition, called the uniform positiveness condition, that may allow inner tolerance $ξ_k\geq 1$ at outer iteration $k$ and can be considerably weaker than the condition $ξ_k\leqξ<1$ with $ξ$ a constant not near one commonly used in literature. We consider the convergence of the inexact RQI with the unpreconditioned and tuned preconditioned MINRES method for the linear systems. Some attractive properties are derived for the residuals obtained by MINRES. Based on them and the new general theory, we make a more refined analysis and establish a number of new convergence results. Let $\|r_k\|$ be the residual norm of approximating eigenpair at outer iteration $k$. Then all the available cubic and quadratic convergence results require $ξ_k=O(\|r_k\|)$ and $ξ_k\leqξ$ with a fixed $ξ$ not near one, respectively. Fundamentally different from these, we prove that the inexact RQI with MINRES generally converges cubically, quadratically and linearly provided that $ξ_k\leqξ$ with a constant $ξ<1$ not near one, $ξ_k=1-O(\|r_k\|)$ and $ξ_k=1-O(\|r_k\|^2)$, respectively. Therefore, the new convergence conditions are much more relaxed than ever before. The theory can be used to design practical stopping criteria to implement the method more effectively. Numerical experiments confirm our results.

preprint2009arXiv

A Refined Harmonic Lanczos Bidiagonalization Method and an Implicitly Restarted Algorithm for Computing the Smallest Singular Triplets of Large Matrices

The harmonic Lanczos bidiagonalization method can be used to compute the smallest singular triplets of a large matrix $A$. We prove that for good enough projection subspaces harmonic Ritz values converge if the columns of $A$ are strongly linearly independent. On the other hand, harmonic Ritz values may miss some desired singular values when the columns of $A$ almost linearly dependent. Furthermore, harmonic Ritz vectors may converge irregularly and even may fail to converge. Based on the refined projection principle for large matrix eigenproblems due to the first author, we propose a refined harmonic Lanczos bidiagonalization method that takes the Rayleigh quotients of the harmonic Ritz vectors as approximate singular values and extracts the best approximate singular vectors, called the refined harmonic Ritz approximations, from the given subspaces in the sense of residual minimizations. The refined approximations are shown to converge to the desired singular vectors once the subspaces are sufficiently good and the Rayleigh quotients converge. An implicitly restarted refined harmonic Lanczos bidiagonalization algorithm (IRRHLB) is developed. We study how to select the best possible shifts, and suggest refined harmonic shifts that are theoretically better than the harmonic shifts used within the implicitly restarted Lanczos bidiagonalization algorithm (IRHLB). We propose a novel procedure that can numerically compute the refined harmonic shifts efficiently and accurately. Numerical experiments are reported that compare IRRHLB with five other algorithms based on the Lanczos bidiagonalization process. It appears that IRRHLB is at least competitive with them and can be considerably more efficient when computing the smallest singular triplets.