Block Newton Method and Block Rayleigh Quotient Iteration for computing invariant subspaces of general complex matrices

Abstract

We consider the Block Newton Method and a modification of it, the Block Rayleigh Quotient Iteration, for approximating a simple p-dimensional invariant subspace X=im(X) and the corresponding eigenvalues collected in the projection L=(XHX)−1XHAX of an arbitrary complex matrix A∈Cn×n. Both methods generate a sequence of bases {Uk} such that the subspaces im(Uk) approximate the target subspace X. The Block Newton Method is Newton's method applied to the extended system AX−XL=0, WHX=I with respect to (X,L) where W is a normalization matrix. We give a direct convergence analysis which, in contrast to the approach by Beyn/Kless/Thümmler [2001], explicitly exploits the product form of the second derivative and, therefore, leads to new and essentially sharper error recursions for each of the sequences {Uk} and {Θk} where Uk and Θk are the k-th X and L iterates, resp. Then, following Lösche/Schwetlick/Timmermann [1998] where the case of real symmetric A is considered, we modify the Block Newton Method as follows: In step k we set W=Uk and use only the X-update from the Newton method to improve Uk and then, with the improved basis, perform a Rayleigh–Schur process and take the orthonormal basis and the generalized Rayleigh quotient matrix obtained from it as new iterates. So the method can be considered as a Block Rayleigh Quotient Iteration. We show that for a sufficiently good initial approximation U0 the method converges in that the sequence {sin⁡ξk} with ξk=∡(im(Uk),X) converges Q-quadratically to zero provided that the invariant subspace X is simple, i.e., that the eigenvalues corresponding to X which are the eigenvalues of the projection L=(XHX)−1(XHAX) of A onto X are separated from those corresponding to X⊥ which are given by λ(A)−λ(L). The convergence results are directly proven – without exploiting standard Newton convergence results – by using the angles between the corresponding subspaces instead of working with the norm differences of the bases as done by Beyn/Kless/Thümmler [2001]. As a byproduct the method is shown to converge cubically if A=AH as in the case p=1.