<s>[INST] Find the complex roots of $4x^2+16$ using QR method $$p(x)=4x^2+16$$
$$m(x)=\frac{p(x)}{4}=x^2+4$$
Now, I can determine the Companion matrix of $x^2+4$:
$$A=A_1=\begin{pmatrix} 0 \ -4 \\ 1 \ \ \ \ \ \ 0 \end{pmatrix}$$

QR decomposition (using Octave):
[Q,R]=qr(A)

$$Q_1\begin{pmatrix} 0 \ -1 \\ -1  \ \  \ \ 0 \end{pmatrix}$$ 
$$R_1=\begin{pmatrix} -1 \ \ \  0 \\ 0  \ \  \ \ 4 \end{pmatrix}$$

$$A_2=R_1 Q_1=\begin{pmatrix} 0 \ 1 \\ -4     \ 0 \end{pmatrix}$$ 
$$R_2=\begin{pmatrix} -4 \ \ \  0 \\ 0  \ \  \ \ 1 \end{pmatrix}$$
$$Q_2=\begin{pmatrix} 0 \ \ \  1 \\ 1 \ \  \ \ 0 \end{pmatrix}$$

$$A_3=R_2 Q_2=\begin{pmatrix} 0 \ -4 \\ 1     \ 0 \end{pmatrix}$$ 
$$R_3=\begin{pmatrix} -1 \ \ \  0 \\ 0  \ \  \ \ 4 \end{pmatrix}$$
$$Q_3=\begin{pmatrix} 0 \ \ \  -1 \\ -1  \ \  \ \ 0 \end{pmatrix}$$

How can I find the complex roots $2i$,$-2i$ (in general, complex roots of polynomials) using this method?
Thanks!
 [/INST] With the real QR algorithm, you only get a upper triangular matrix as limit if all eigenvalues are real. With complex eigenvalues, you get $2×2$ matrix blocks on the diagonal for any complex conjugate pair.
In the complex case, you get a unitarian matrix for $Q$ by employing complex shifts. A first strategy is to compute the shift as eigenvalue of the bottom-most $2×2$ diagonal block. Which again makes the matrix formulation for degree $2$ one big circular computation.
</s>