Convergence of the complex QR algorithm to Schur decomposition

I study the complex Schur decomposition of a complex matrix $$A \in \mathbb{C}^{n \times n}$$, that is: $$A = U T U^H$$ where $$T$$ is upper-triangular (the eigenvalues of $$A$$ appear on its diagonal, and may be complex) and $$U$$ is unitary $$U^H = U^{-1}$$ (recall that $$U^H = U^* = \bar{U}^T$$). Such a decomposition always exists.

The standard way to compute the Schur decomposition (at least, the one used in all the sources I found, including Golub & Van Loan 2013, "Matrix Calculations" and Wikipedia) is using the QR algorithm (see Golub & Van Loan, chapter 7, or Wikipedia). Basically, one repeats the following steps:

1. Set $$T_0 = A$$ and $$U = I$$ (in practice, one starts with reduction of $$A$$ to an Hessenberg form).
2. Compute the QR factorization $$T_k = Q_k R_k$$.
3. Form $$T_{k+1} = R_k Q_K = Q_k^H Q_k R_k Q_k = Q_k^H T_k Q_k$$.
4. Accumulate $$U = U Q_k$$.

Then, in most cases (and this is an issue I need clarifications about) $$T_{k+1} \to T$$ converges to the upper-triangular matrix $$T$$ of the Schur decomposition and $$U = \lim_{k \to \infty} Q_0 Q_1 \cdots Q_k$$.

In this question the poster quotes a theorem that ensure convergence: let $$|\lambda_1| > ... > |\lambda_n | \ge 0$$ be the eigenvalues of $$A$$. Then $$(T_{k})_{i,j} = O\left(\left(\dfrac{\vert \lambda_{i} \vert}{\vert \lambda_{j} \vert}\right)^k\right)$$ for $$i > j$$.

This means that if there are $$\lambda_{i_1}$$ and $$\lambda_{i_2}$$ such that $$|\lambda_{i_1}| = |\lambda_{i_2}|$$, or they are very close, $$\left| \frac{|\lambda_{i_1}|}{|\lambda_{i_2}|} - 1 \right| < \mathrm{tol}$$ (here tol is the computer's numerical tolerance) the matrices $$T_k$$ will fail to converge to $$0$$ in the element $$(T_k)_{i_1,i_2}$$ and thus won't converge to an upper-triangular matrix.

I have written a code to implement the QR algorithm for complex matrices, and encountered this phenomena: a matrix that failed to convergence to upper-triangular with eigenvalues almost the same in their absolute value.

Python's SciPy package has a special function $$\texttt{scipy.linalg.schur}$$ that computes the Schur decomposition even for such bad matrices. I have three questions:

1. How to overcome this problem in the framework of the QR algorithm? Or is it impossible and I am forced to accept that for some matrices the QA algorithm just won't work?
2. What are the necessary and/or sufficient condition for the convergence of the QR algorithm to the Schur decomposition of $$A \in \mathbb{C}^{n \times n}$$?
3. How Python's SciPy schur function (which uses LAPACK function at the back stage) overcome this problem and compute Schur decomposition even for such bad matrices?

Related Questions:

• This problem appears also for $\lambda(A) = \{0.62534345+0.61656135i, -0.25075095+0.07059212i, 0.07124168+0.01315155i, 0.2578328 -0.04137332i\}$ which satisfy $|\lambda|(A) = \{ 0.87818126, 0.26049815, 0.07244543, 0.2611312 \}$. We see $A$ has 2 eigenvalues close to each other $|\lambda| \approx 0.26$ but still not as close as $10^{-15}$. Commented May 31, 2020 at 8:21