where $\boldsymbol{Q}_{(i,j)}$ is a Householder transformation matrix; $\boldsymbol{I}$ is the identity matrix; the scalar $\tau_{(i,j)}$ and vector $\boldsymbol{v}_{(i,j)}$ are computed following the standard Householder method [@numerical_press_2007; @householder_householder_1958]. Most Householder vectors $\boldsymbol{v}_{(i,j)}$ have a length equal to $b$, except that $\boldsymbol{v}_{(i,j)}$ generated in the last sweep of each stage may be shorter. Matrix-vector operations are still needed in order to apply these Householder transformations, but the computational cost is much less than in the one-stage tridiagonalization algorithm, as the vectors are much shorter. The left panel of Fig. \[fig:step4\] shows where vectors $\boldsymbol{v}_{(i,j)}$ are generated with an example of $N=17$, $b=4$, and $N-2=15$ stages of bulge chasing. In the actual code, the Householder matrices $\boldsymbol{Q}_{(i,j)}$ are never explicitly constructed. Instead, $\tau_{(i,j)}$ and $\boldsymbol{v}_{(i,j)}$ are stored and used for the transformations in Eqs. \[eq:2stage\_fwd2\] and \[eq:2stage\_bkwd1\], where $\boldsymbol{Q}$ is the product of all $\boldsymbol{Q}_{(i,j)}$.

![image](fig/step4.png){width="99.00000%"}

After the tridiagonal eigenproblem (Eq. \[eq:2stage\_solve\]) is solved, we get the eigenvalues $\boldsymbol{\Sigma}$ and eigenvectors $\boldsymbol{X}$ of the tridiagonal matrix $\boldsymbol{T}$. The main task of the tridiagonal-to-banded back-transformation is to apply all $\boldsymbol{Q}^*_{(i,j)}$ to $\boldsymbol{X}$ (Eq. \[eq:2stage\_bkwd1\]). It is obvious that the eigenvectors can be back-transformed independently, leading to a trivial parallelism over eigenvectors. Moreover, as the 