4
$\begingroup$

After having googled for several days without locating a definitive answer, I will try my luck here!

I have implemented a version of the QR algorithm to calculate Eigenvalues and hopefully Eigenvectors of a matrix $A$ of dimension $n\times n$.

In order to speed up the convergence rate i have applied a version of the algorithm that utilizes several improvenments, mainly inspired by https://www.math.kth.se/na/SF2524/matber15/qrmethod.pdf.

Firstly, we calculate $H=Hessenberg(A)$, which transforms $A$ to upper Hesenberg form $H_{n\times n}$.

Now, using $H$ we estimate the Eigenvalues using a QR-algorithm applying a Wilkinson shift and deflation. This can loosely be described as the following pseudo procedure, where $\lambda_i$ denotes the i'th Eigenvalue:

$\text{set}\ H_0:=H\\ \text{for}\ m=n,\ldots,2 \ \text{do} \\ \quad k=0\\ \qquad \text{repeat}\\ \qquad \quad k=k+1\\ \qquad \quad \sigma_k=Wilkinson(H_{k-1})\\ \qquad \quad H_{k-1}-\sigma_kI=:Q_kR_k\quad (*)\\ \qquad \quad H_k=R_kQ_k+\sigma_kI\\ \qquad \text{until}\ \vert h^{(k)}_{m,m-1}\vert<\epsilon\\ \qquad \lambda_m=h^{(k)}_{m,m}\\ \qquad H^{(0)}=H^{(k)}_{1:(m-1),1:(m-1)}\\ \text{end for}$

The last step is simply a deflation step that drops the last row and column of $H$ upon satisfactory convergence towards an eigenvalue.

The function $Wilkinson$ calculates the shift, and the QR factorization $(*)$ is done using Givens rotations, which is a standard procedure.

My question is, how do i determine the corresponding eigenvectors? In the standard QR-algorithm this can be computed as $\Pi _iQ_i$, however due to the deflation step I don't know how to proceed?

I have verified that eigenvalues are calculated correctly.

In advance, thank you very much for any help!

EDIT: Implementing Federico's solution below works as intended. Make sure you use your initial similarity transform $H=UAU^*$, where $H$ is upper Hessenberg, i.e. let $\bar{Q}_H$ be the eigenmatrix of $H$ then $\bar{Q}_A=U^*\bar{Q}_H$ yields the eigenmatrix of $A$.

$\endgroup$
  • $\begingroup$ What is $\Pi_i Q_i$? $\endgroup$ – Federico Poloni Jan 5 '17 at 17:38
1
$\begingroup$

Instead of dropping one row and one column, compute at each step a $(n-1)\times(n-1)$ orthogonal transformation (or $(n-k)\times(n-k)$, after $k$ deflation steps) $Q$ by working to the reduced matrix, and then apply it to the full matrix as $$ \begin{bmatrix} Q^* \\& I \end{bmatrix} \begin{bmatrix} H_{11} & H_{12}\\ 0 & H_{22} \end{bmatrix} \begin{bmatrix} Q \\& I \end{bmatrix} = \begin{bmatrix} Q^*H_{11}Q & Q^*H_{12}\\ 0 & H_{22} \end{bmatrix}. $$ In practice all you have to do is operating on the leading $(n-k)\times(n-k)$ block as you were doing before, and then multiplying $H_{12}$ by the orthogonal transformation $Q$ that you have generated.

In this way, your algorithm computes explicitly a sequence of $n\times n$ orthogonal transformations $Q_1, Q_2, \dots, Q_m$ that turns $A$ into a triangular matrix (Schur form). You can accumulate the product $Q_1Q_2\dotsm Q_m$ with $O(n^2)$ additional operations per step (so $O(n^3)$ in total during the algorithm, under the usual assumptions that $O(1)$ iterations per eigenvalue are sufficient). After that, all you have to do is recover the eigenvectors from the Schur form.

I hope this is sufficiently clear!

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you, I will try to alter the algorithm tomorrow - if I can get it to work I will accept your answer. If I understand your answer correctly I (in theory) already have all available components at disposal and it's just a matter of keeping track of the right hand side of the equation? $\endgroup$ – Rasmus Jan 5 '17 at 21:15
  • $\begingroup$ @Rasmus There is no right-hand side here. To get eigenvectors, you still have to implement three things, all reasonably easy: (1) update at every step also the top-right block of the $n\times n$ working matrix instead of focusing only the top-left one, (2) compute the product of all orthogonal transformations one by one while you apply them, and (3) compute the eigenvector matrix of the triangular matrix in the Schur form (and multiply it by the matrix in (2)). $\endgroup$ – Federico Poloni Jan 5 '17 at 21:31
0
$\begingroup$

@Federico Poloni I apologize for being a bit slow here, how do I calculate the product of the orthogoanl matrices $Q_1,Q_2,\ldots,Q_m$, when at some point during the algorithm a deflation is performed which also reduces the dimension of the applied $Q$'s accordingly.

To clear things for me, I will introduce a double iteration index $m,k$, where $k$ corresponds to th repeatedly applied QR-factorization before satisfactory convergence towards eigenvalue $\lambda_m$, $m=1,\ldots,n$, which is associated with the outer for loop. Hence, denote for each iteration of the QR-factorization $Q_{m,k}$ which corresponds to iteration $k$ looking for eigenvalue $m$.

Now, whenever convergence is obtained for some $m$ we focus on the leading $(m-1)\times(m-1)$-block matrix and hence $\dim Q_{m,k}=m >\dim Q_{m-1,k}, \ k=1,\ldots$. However, as I understand I ought to calculate the product $\prod_{m=n}^1\prod_kQ_{m,k}$?

How can I alter my orthogonal matrices $Q_{m-1,k}$ so I can multiply them with $Q_{m,k}$ in the correct way?

I hope you understand my somewhat clumsy question - and I should certainly have paid more attention to linear algebra.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Pad them with identities, as in the displayed equation in my answer. Don't think of it as "deflating", think "working on the top-left block". If you think to them in terms of Givens rotations, it's more natural (and in any case you have to compute this product one Givens rotation at a time anyway, for complexity reasons). $\endgroup$ – Federico Poloni Jan 6 '17 at 9:21
  • 2
    $\begingroup$ Also, please try not to use answers for something that is not an answer, on this site. Yours should have been a comment to my answer, or, if you needed more characters, a p.s. added to your question with an edit. $\endgroup$ – Federico Poloni Jan 6 '17 at 9:23
  • $\begingroup$ I will keep that in mind :-) Okay, so if i define $\bar{Q}^{(j)}$ as the product of all the preceeding orthogonal matrices $Q$ i can "update" this for some $Q_{m,k}$, $m<n$ similarly to your answer where the center matrix consists of the submatrices that form $\bar{Q}^{(j)}$, namley $\bar{Q}_{11},\bar{Q}_{12},\bar{Q}_{21},\bar{Q}_{22}$ and $\dim \bar{Q}_{11} = \dim Q_{m,k}$? $\endgroup$ – Rasmus Jan 6 '17 at 11:11
  • $\begingroup$ No, this "update" is from one side only. If you make similarity transformations $A_1=Q_1^*AQ_1$, $A_2=Q_2^*A_1Q_2$, $A_3=Q_3^*A_2Q_3$ they compose as $A_3 = (Q_1Q_2Q_3)^*A(Q_1Q_2Q_3)$, so you have to keep track of the matrix $\bar{Q} = Q_1Q_2Q_3$, updating it at each step with a one-sided product only. (Each $Q_i$ may be of the form $\begin{bmatrix}* & 0 \\ 0 & I\end{bmatrix}$ here.) $\endgroup$ – Federico Poloni Jan 6 '17 at 13:46
  • $\begingroup$ First of all thank you very much for taking your time to help me! I have implemented this approach, and at the end of my algorithm the (columns) of $\bar{Q}$ almost converges to the correct eigenvectors of the upper Hessenberg matrix $H$. I say almost as the column vectors somehow have incorrect signs for only some of the entries - however the numerical entry-values are correct. Yet, I have no idea what causes this "arbitrary" sign-flipping. Finally, if $H=hessenberg(A)$ then the eigenmatrices of $H$ and $A$ do not coincide. What is the relationship between the two? $\endgroup$ – Rasmus Jan 6 '17 at 14:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.