Reorthogonalized randomized power subspace iteration

Input $F\in \{\mathbb{R},\mathbb{C}\}$, a fixed matrix $A\in F^{m\times n}$, integers $\mathrm{\ell }\ge 1$ and $q\ge 0$, and draw $\mathrm{\Omega }\in F^{n\times \mathrm{\ell }}$ with iid standard real or normalized complex Gaussian entries as in Standard real and normalized complex Gaussian matrix laws. Put $Q_0$ equal to the orthonormal basis selected by the finite scan in Exact Gaussian randomized range finder for $\mathrm{Ran}(A\mathrm{\Omega })$. For $j=1,\dots ,q$, form $A^{\ast }{Q}_{j-1}$ and select ${\tilde{Q}}_j$ from its range by the same scan, then form $A{\tilde{Q}}_j$ and select $Q_j$ from its range by that scan. Return $Q_q$. Thus exact QR/range orthogonalization follows each application of $A^{\ast }$ and $A$; empty or rank-deficient ranges return bases of their actual ranges. The procedure is defined in exact arithmetic and makes no finite-precision error guarantee.