Strong rank-revealing QR with finite swaps and explicit statuses

Input a real or complex matrix $A\in {\mathbb{F}}^{m\times n}$ with $m\ge n$, a target $k$ with $0\le k<n$, a factor $f\ge 1$, and a finite swap cap $S\in \mathbb{N}\cup \{0\}$. The exact-arithmetic version is specified below. In finite precision, $A$ denotes the represented matrix: each real and imaginary component, and $f$, is a finite binary floating-point number, interpreted as its exact dyadic rational value. The arithmetic format and rounding rule are fixed for a run. A nonfinite input or required rounded intermediate returns range-failure. This mode makes no rank claim about an unspecified matrix before input rounding.

For the finite-precision rank test, copy the represented $A$ into exact rational complex arithmetic. Scan columns from left to right with a current pivot row, select the first nonzero entry at or below that row, interchange rows, and eliminate entries below it by exact rational row operations. If no such entry exists, skip that column. Increment the pivot row only after a nonzero pivot. The number $r$ of pivots is the reported exact represented-input rank. Return rank-deficient-input if $r<n$; otherwise continue with rank status full-column-rank. This test uses unbounded-size rational integers, not the rounded QR diagonal or a numerical-rank threshold. In exact mode, use the same finite row-elimination instructions over the exact input field to decide the rank condition in step 1.

The finite-precision mode additionally takes a finite real binary format $\mathcal{F}$ containing $0,1,-1$ and closed under sign changes, a fixed rounding rule, and a certificate-pass budget $C\in {\mathbb{N}}_0$ (default $1$; $C=0$ disables certification). Real and imaginary parts are stored separately. For a real rational result within the format's finite range, $\mathrm{fl}$ returns its specified nearest representable value, with a fixed tie rule, or its specified directed rounding; a result outside that range returns range-failure. Exact zeros are copied as zero. Real addition, subtraction, multiplication, and division mean their exact rational result followed by this rounding; check each divisor before dividing. A real square root of a nonnegative represented number is rounded by comparing that number with squares of the finitely many nonnegative format values: bracket its root, then use the selected directed endpoint or compare with the squared midpoint for nearest rounding. A negative argument or a root beyond the finite range returns range-failure. Thus no exact irrational arithmetic oracle or unbounded precision search is part of this mode. A nonfinite input or any failed required scalar operation immediately returns range-failure with any completed partial data.

Use the following fixed scalar conventions throughout the rounded core. For $z=a+ib$ and $w=c+id$, complex addition/subtraction round their two real components; multiplication returns $\mathrm{fl}(\mathrm{fl}(ac)-\mathrm{fl}(bd))+i\mathrm{fl}(\mathrm{fl}(ad)+\mathrm{fl}(bc))$. For division compute $d_0=\mathrm{fl}(c^2)+\mathrm{fl}(d^2)$ with the addition rounded, require $d_0>0$, compute the rounded product $z\bar{w}$, and divide its two real components by $d_0$, rounding each. Conjugation and sign changes are exact copies. For a vector, compute its squared norm by rounding each real-component square and summing them with rounding in increasing coordinate order, real part before imaginary part; apply the rounded nonnegative square root. Dot products use the above rounded complex products and a left-to-right rounded sum. These rules specify every norm and score calculation below; a nonzero vector whose required computed norm is zero returns range-failure.

The initial permutation is $\mathrm{\Pi }=I_n$, with original labels $1,\dots ,n$. A rounded full QR sweep of a fixed current column order starts from $B=A\mathrm{\Pi }$, $\widehat{Q}=I_m$ and an empty transformation record. During the initial pivoted sweep, at column $j=1,\dots ,n$ compute the rounded norms of $B(j:m,\mathrm{\ell })$ for $\mathrm{\ell }=j,\dots ,n$, choose the largest, breaking ties by smallest original column label, and swap the full columns of $B$ and their labels. A fixed-order sweep skips these pivot searches. In either sweep, if $B(j+1:m,j)=0$ entrywise, record no transformation and continue. Otherwise copy $z=B(j:m,j)$, compute its rounded norm $r>0$, and set $\varphi =1$ if $z_1=0$; if $z_1\mathrm{\ne }0$, compute its rounded modulus $a>0$ and divide its two components by $a$ to form $\varphi$. Set $\alpha =-\mathrm{fl}(r\varphi )$, $v_1=\mathrm{fl}(z_1-\alpha )$ and $v_i=z_i$ for $i>1$. Compute the rounded squared norm $\widehat{d}$ of $v$, require $\widehat{d}>0$, and put $\widehat{\tau }=\mathrm{fl}(2\mathrm{/}\widehat{d})$. Record the active index $j$ and the exact dyadic components of this nonzero stored $v$ before updating anything.

For each column $\mathrm{\ell }=j,\dots ,n$, compute the rounded dot product $a_{\mathrm{\ell }}=v^{\ast }B(j:m,\mathrm{\ell })$ from the old column, set $h_{\mathrm{\ell }}=\mathrm{fl}(\widehat{\tau }{a}_{\mathrm{\ell }})$, then replace each active entry by $\mathrm{fl}(B_{i\mathrm{\ell }}-\mathrm{fl}(v_{i-j+1}{h}_{\mathrm{\ell }}))$. For each row $a=1,\dots ,m$ of $\widehat{Q}$, compute $b_a=\widehat{Q}(a,j:m)v$ from its old row, set $g_a=\mathrm{fl}(\widehat{\tau }{b}_a)$, and replace each active entry by $\mathrm{fl}({\widehat{Q}}_{ai}-\mathrm{fl}(g_a{\bar{v}}_{i-j+1}))$. Finally set $B(j+1:m,j)$ to exact zeros by assignment. At the end return $\widehat{R}=B$, the full $m\times m$ candidate $\widehat{Q}$, the final permutation and the ordered record. The imposed zeros are part of this rounded procedure, not a claim of exact annihilation; the certificate measures their effect as well as all rounding. There are at most $n\le mn$ reflector updates per sweep. Empty active tails are skipped, and no norm or quotient is required for them.

In finite precision step 1 uses this pivoted sweep, and step 2 uses a fresh fixed-order sweep on the represented input $A$ in its current order. Each fresh sweep resets $\widehat{Q}$ and its transformation record. Before the triangular solves, decide singularity of its represented ${\widehat{R}}_{11}$ by the same exact rational elimination as the input rank test; if singular, return singular-selected-block. Otherwise solve ${\widehat{R}}_{11}X={\widehat{R}}_{12}$ and ${\widehat{R}}_{11}Y=I_k$ by rounded back-substitution: for each right-hand-side column, process rows $i=k,\dots ,1$, sum ${\widehat{R}}_{ih}{x}_h$ in increasing $h=i+1,\dots ,k$, subtract from the supplied right-hand side, check the diagonal divisor and use the complex division above. Compute ${\omega }_i$ as the rounded norm of row $i$ of $Y$, ${\gamma }_j$ as the rounded norm of column $j$ of ${\widehat{R}}_{22}$, and $c_{ij}$ as the rounded modulus of $X_{ij}$. Compute $t_{ij}=\mathrm{fl}({\omega }_i{\gamma }_j)$ and $s_{ij}=\mathrm{fl}(\sqrt{\mathrm{fl}(\mathrm{fl}(c_{ij}^2)+\mathrm{fl}(t_{ij}^2))})$. Scan these finite candidate lists for the step-4 maxima and its stated tie rules. The diagnostics $(X,Y,({\omega }_i),({\gamma }_j),(c_{ij}),(s_{ij}))$ are rounded candidates, not exact inverse or score identities. For $k=0$ there are no such solves or score lists.

The numbered instructions below specify exact arithmetic; the preceding rounded instructions replace the corresponding factorization, solve and score operations in finite precision.

1. If $A$ is not full column rank, return rank-deficient-input. Otherwise initialize a column permutation $\mathrm{\Pi }$ and the full column-pivoted QR factorization $A\mathrm{\Pi }=QR$, where $Q\in {\mathbb{F}}^{m\times m}$ is unitary and $R\in {\mathbb{F}}^{m\times n}$ is upper trapezoidal. At each pivot choose a remaining column of maximal residual $2$-norm, resolving ties by the smallest original column index. Partition $R=\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right]$ after the first $k$ rows and columns. Do not return the reduced QR form. The selected first $k$ columns then have a nonsingular leading block. For $k=0$, return this full factorization with success.
2. For $k>0$, refactor the current fixed column order in the full QR form, without further pivoting whenever needed, using unitary Householder or Givens transformations. If the selected leading block is singular, return singular-selected-block.
3. Form $X=R_{11}^{-1}{R}_{12}$. For $1\le i\le k$ let ${\omega }_i=\mathrm{\parallel }{e}_i^{\ast }{R}_{11}^{-1}{\mathrm{\parallel }}_2$; for $1\le j\le n-k$ let ${\gamma }_j=\mathrm{\parallel }{R}_{22}(:,j){\mathrm{\parallel }}_2$. Define the swap score $s_{ij}=(\mathrm{\mid }{X}_{ij}{\mathrm{\mid }}^2+({\omega }_i{\gamma }_j{)}^2{)}^{1\mathrm{/}2}$.
4. If ${\max }_{i,j}\mathrm{\mid }{X}_{ij}\mathrm{\mid }>f$, choose the lexicographically first pair attaining the largest $\mathrm{\mid }{X}_{ij}\mathrm{\mid }$. Otherwise, if ${\max }_{i,j}{s}_{ij}>f$, choose the lexicographically first pair attaining the largest score. If neither inequality holds, return success and the factors, permutation, $X$, and maximum score. The first branch retains the usual coefficient-triggered swap; the second controls the trailing residual needed for the strong spectral comparison.
5. When a violating pair was chosen, if the swap count is already $S$, return swap-budget-exhausted with the current partial factors. Otherwise interchange the chosen selected and trailing columns in the current permutation. In exact arithmetic restore triangular form with Givens rotations; in finite precision discard the old factors and record, since step 2 recomputes them from $A$ in the new fixed order. Increment the swap count and return to step 2.

The statuses are success, rank-deficient-input, singular-selected-block, swap-budget-exhausted, and range-failure. In particular, ordinary column-pivoted QR is only the initialization: the procedure does not assert that its output alone satisfies the strong bounds. Both modes return the full factor shape and available diagnostics. A finite-precision success is a candidate until the following separate certificate is verified. The rounded dense $\widehat{Q}$ is never asserted to be unitary.

For each recorded active dyadic vector $v$, define the embedded exact rational matrix $H=I-\frac{2v{v}^{\ast }}{v^{\ast }v}$ on its active coordinates, with identity on the preceding coordinates. Its denominator is a strictly positive exact rational. It is exactly unitary: writing $P=v{v}^{\ast }\mathrm{/}(v^{\ast }v)$ gives $P^{\ast }=P$, $P^2=P$ and $(I-2P{)}^{\ast }(I-2P)=I$. Thus the ordered sweep record represents the exact unitary product $W=H_1\cdots H_t$; an empty record represents $I_m$. This product corresponds to the final fresh sweep only, and is distinct from the rounded $\widehat{Q}$ and the nonunitary rounded update coefficients $\widehat{\tau }$. Its rational representation requires no normalization by an irrational vector norm.

Certification is finite optional postprocessing, independent of the core status. For a nonsuccess status assign not-applicable. On success with $C=0$, assign uncertified. With $C\ge 1$ perform the following one exact rational pass, using unbounded-size integer numerators and denominators; no retry until success is allowed.

• Check factor dimensions, that $\mathrm{\Pi }$ is a permutation of $1,\dots ,n$, that $\widehat{R}$ is exactly upper trapezoidal, and that each recorded vector has the declared active dimension and is nonzero. Reconstruct each $H$ by the displayed rational formula, multiply them in the recorded order, and verify entrywise $W^{\ast }W=I_m$. Any failed check sets uncertified.
• For $k>0$, apply exact rational elimination to ${\widehat{R}}_{11}$. If singular, the certificate is uncertified. Otherwise compute its exact inverse $Y_0$ by elimination on $[{\widehat{R}}_{11}{I}_k]$, normalizing each nonzero pivot and eliminating above and below it. Form $Z=Y_0{\widehat{R}}_{12}$ and the exact nonnegative rational score squares $q_{ij}=\mathrm{\mid }{Z}_{ij}{\mathrm{\mid }}^2+({\sum }_{a=1}^k\mathrm{\mid }(Y_0{)}_{ia}{\mathrm{\mid }}^2)({\sum }_{b=1}^{m-k}\mathrm{\mid }({\widehat{R}}_{22}{)}_{bj}{\mathrm{\mid }}^2)\mathrm{.}$ Set $b_{ij}=\sqrt{q_{ij}}$, stored as a symbolic nonnegative square root, and check the rational inequalities $q_{ij}\le f^2$ for every $i,j$. These are equivalent to all $b_{ij}\le f$. They certify the exact scores of the represented blocks, independently of rounded $X,Y,\omega ,\gamma ,s$. For $k=0$ this list and every inverse condition are absent.
• Compute exactly $E=A\mathrm{\Pi }-W\widehat{R}$ and $e^2={\sum }_{a,b}\mathrm{\mid }{E}_{ab}{\mathrm{\mid }}^2$. Set $\epsilon =\sqrt{e^2}$, again a symbolic nonnegative square root; then $\mathrm{\parallel }A\mathrm{\Pi }-W\widehat{R}{\mathrm{\parallel }}_2\le \mathrm{\parallel }E{\mathrm{\parallel }}_F=\epsilon$. If all structural and score checks passed, return certified with the rational record, reconstructed $W$, $(q_{ij})$, $(b_{ij})$, $E,e^2,\epsilon$ and their check outcomes; otherwise return uncertified with any completed data. There is no extra residual tolerance implicit in the certificate: its actual residual bound $\epsilon$ is reported, however large, and any application requiring a small bound must check that separately.

This single finite pass uses at most the prescribed positive pass budget; a supplied certificate is accepted only by recomputing these exact checks. Certification does not change the core status or rounded candidates. It gives precisely a nonsingular ${\widehat{R}}_{11}$ when $k>0$, rigorous combined-score bounds $b_{ij}\le f$, and an exact unitary representation $W$ with the reported reconstruction bound. The bounds concern the represented input $A$; a claim about different intended data also needs its separately supplied input-error bound. Ordinary rounded score or residual evaluations do not establish this certificate.