Computed singular values and rank revealing factorizations

Let $A\in {\mathbb{F}}^{m\times n}$, where $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$. Put $p=\min (m,n)$ and $M=\max (m,n)$. If $m\ge n$, set $C=A$; if $m<n$, set $C=A^{\ast }$. Thus $C\in {\mathbb{F}}^{M\times p}$.

Save the input as $C_0=C$ and initialize $U=I_M$, $V=I_p$, and empty ordered lists of stored local left and right transformations. For $k=1,\dots ,p$, in this order:

1. Let $x=C_{k:M,k}$. If $x=0$, use the identity left transformation. Otherwise choose $\varphi =x_1\mathrm{/}\mathrm{\mid }{x}_1\mathrm{\mid }$ when $x_1\mathrm{\ne }0$ and $\varphi =1$ when $x_1=0$, set $\alpha =-\varphi \mathrm{\parallel }x{\mathrm{\parallel }}_2$, $v=x-\alpha e_1$, and $h=v\mathrm{/}\mathrm{\parallel }v{\mathrm{\parallel }}_2$. Use the Householder reflector $H=I-2h{h}^{\ast }$ on rows $k:M$. Update $C\leftarrow H^{\ast }C$ on those rows and append this local transformation to the left list. In exact arithmetic update $U\leftarrow U\mathrm{diag}(I_{k-1},H)$. In finite precision update columns $k:M$ of the dense accumulator $U$ by right application of the same stored local reflector data, using the separate-block branch of Scaled real and complex Householder and Givens update arithmetic and making no target-zero write to $U$.
2. If $k<p$, let $y=C_{k,k+1:p}^{\ast }$. If $y=0$, use the identity right transformation. Otherwise apply the same phase and normalization rule to $y$ to form the Householder reflector $G$, then apply $C\leftarrow C\mathrm{diag}(I_k,G)$ and append this local transformation to the right list. In exact arithmetic update $V\leftarrow V\mathrm{diag}(I_k,G)$. In finite precision update columns $k+1:p$ of the dense accumulator $V$ by right application of the same stored local reflector data, using the separate-block branch of Scaled real and complex Householder and Givens update arithmetic and making no target-zero write to $V$.

Set $B$ to the leading $p\times p$ block of the final work matrix $C$, and return $B,U,V$ together with the ordered local transformation lists.

In finite precision, each left/right reflector uses Scaled real and complex Householder and Givens update arithmetic, including its real/complex scaled formation, range checks and exact target-column/row zero writes. This is the arithmetic realization of the exact steps above. A zero active vector takes the identity branch. Any helper range-failure, invalid operation, nonfinite value, or nonzero exact intermediate outside the normal range returns range-failure with the partial factors and stored local transformations; otherwise return success, $B,U,V$, the ordered stored local transformations, the reconstruction residual $\mathrm{\parallel }{C}_0-U\left[\begin{array}{c}B\\ 0\end{array}\right]V^{\ast }{\mathrm{\parallel }}_F$, and the two orthogonality residuals $\mathrm{\parallel }{U}^{\ast }U-I{\mathrm{\parallel }}_F$ and $\mathrm{\parallel }{V}^{\ast }V-I{\mathrm{\parallel }}_F$. The reconstruction uses the saved input $C_0$, not the overwritten work matrix. The exact products of the stored unitary representatives are the theoretical factors; the dense finite-precision accumulators are rounded and their residuals are computed quantities, not certified bounds by themselves. Matrix residuals use the Frobenius norm defined in The Frobenius norm $\parallel A{\parallel }_F=({\sum }_{i,j}\mathrm{\mid }{a}_{ij}{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}$ on real or complex matrices.

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.

Let $X\in {\mathbb{F}}^{m\times n}$, $Y\in {\mathbb{F}}^{n\times q}$, and $Z\in {\mathbb{F}}^{r\times m}$ over $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$. For every $1\le i\le \min (m,n,q)$ and $1\le j\le \min (r,m,n)$, ${\sigma }_i(XY)\le {\sigma }_i(X)\mathrm{\parallel }Y{\mathrm{\parallel }}_2,{\sigma }_j(ZX)\le \mathrm{\parallel }Z{\mathrm{\parallel }}_2{\sigma }_j(X)\mathrm{.}$

Let $A,E\in {\mathbb{F}}^{m\times n}$, where $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$, and put $p=\min (m,n)$. Then for every $1\le i\le p$, $\mathrm{\mid }{\sigma }_i(A+E)-{\sigma }_i(A)\mathrm{\mid }\le \mathrm{\parallel }E{\mathrm{\parallel }}_2\mathrm{.}$

Let $m,n\ge 1$, $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, $A\in {\mathbb{F}}^{m\times n}$, $p=\min (m,n)$, and $\tau \ge 0$. For the spectral norm, the numerical rank defined by Numerical rank relative to a declared norm, scale, and tolerance satisfies ${\mathrm{rank}}_{\tau }(A)=\mathrm{\#}\{i\in \{1,\dots ,p\}:{\sigma }_i(A)>\tau \}\mathrm{.}$

Input the original $A\in {\mathbb{F}}^{m\times n}$, the successful output record of Householder bidiagonalization and its finite precision status for $C=A$ if $m\ge n$ and $C=A^{\ast }$ otherwise, a deflation tolerance $\tau \ge 0$, and a finite nonnegative integer sweep budget $N$. The reduction record supplies the upper bidiagonal work matrix $B_0$, dense accumulators $U_B,V_B$, and ordered stored local transformation data. Initialize $B=B_0$, the sweep accumulators $U_{\mathrm{\Sigma }},V_{\mathrm{\Sigma }}$ to identities of order $p=\min (m,n)$, empty ordered left/right transformation lists, an empty deflation record, a pass counter $t=0$, and a sweep counter $s=0$. Write the diagonal as $d_i$ and superdiagonal as $e_i$, and set $e_0=0$.

At the start of each pass, set $t\leftarrow t+1$, set $e_i$ to zero and append $(t,i,e_i^{old})$ to the deflation record whenever $\mathrm{\mid }{e}_i\mathrm{\mid }\le \tau (\mathrm{\mid }{d}_i\mathrm{\mid }+\mathrm{\mid }{d}_{i+1}\mathrm{\mid })$. The corresponding additive matrix change has entry $-e_i^{old}$ at $(i,i+1)$ and zero elsewhere. Record any later non-target sparsity cleanup the same way with its row, column, and old value. The zero entries split the active blocks. If every active block has size one, stop with converged. Otherwise process the active blocks in increasing order of $\mathrm{\ell }$, skipping singleton blocks. Before each block sweep, if $s=N$, stop with budget-exhausted; otherwise form only the tridiagonal entries of $T=B^{\ast }B$ on that block: $T_{ii}=\mathrm{\mid }{d}_i{\mathrm{\mid }}^2+\mathrm{\mid }{e}_{i-1}{\mathrm{\mid }}^2$ and $T_{i,i+1}=\overline{d_i}{e}_i$. Choose $\mu$ to be the eigenvalue of the trailing Hermitian block $\left[\begin{array}{cc}{T}_{r-1,r-1}& T_{r-1,r}\\ T_{r,r-1}& T_{rr}\end{array}\right]$ closest to $T_{rr}$, resolving a tie by the smaller eigenvalue. If the leading pivot $d_{\mathrm{\ell }}$ is zero or has magnitude at most $\tau \max (\mathrm{\mid }{e}_{\mathrm{\ell }-1}\mathrm{\mid },\mathrm{\mid }{e}_{\mathrm{\ell }}\mathrm{\mid })$, use the zero-shift branch $\mu =0$ and choose the initial right Givens factor $G_{\mathrm{\ell }}$ so that $(d_{\mathrm{\ell }},e_{\mathrm{\ell }})G_{\mathrm{\ell }}=(\rho ,0)$, avoiding division by the small or zero diagonal.

Finite-precision shift and comparison schedule. The preceding exact-arithmetic formulas retain their exact meaning. For finite precision use the scalar order, six-real-operation complex product, scaled real-component norm, and first-failing-primitive range policy of Scaled real and complex Householder and Givens update arithmetic also for the following stages. All assignments below are stored; do not fuse operations. Comparisons, maxima, copies, conjugations, sign changes, and real absolute values are exact operations on stored finite values. Evaluate each complex magnitude by that helper's scaled two-component norm, with zero magnitude returned exactly as zero. At each pass visit deflation indices increasingly; at index $i$ evaluate $a_i=\mathrm{\mid }{e}_i\mathrm{\mid }$, $b_i=\mathrm{\mid }{d}_i\mathrm{\mid }$, $c_i=\mathrm{\mid }{d}_{i+1}\mathrm{\mid }$ in that order, then $h_i=\mathrm{fl}(b_i+c_i)$ and $k_i=\mathrm{fl}(\tau h_i)$, and deflate exactly when $a_i\le k_i$.

For each nonsingleton block, form its Gram diagonal entries in increasing index order, followed by its upper off-diagonal entries in increasing index order. To form $T_{ii}$, square $\mathrm{\Re }{d}_i,\mathrm{\Im }{d}_i,\mathrm{\Re }{e}_{i-1},\mathrm{\Im }{e}_{i-1}$ separately with rounding, add the first two rounded squares, add the last two, and add these two rounded sums, rounding each addition. Form $T_{i,i+1}$ by the helper's complex product for $\overline{d_i}{e}_i$ and copy its conjugate to the lower off-diagonal. No other Gram entries are formed. Write the stored trailing block as $\left[\begin{array}{cc}a& z\\ \bar{z}& c\end{array}\right]$. In this order compute $h=\mathrm{fl}(\mathrm{fl}(a+c)\mathrm{/}2)$, $\delta =\mathrm{fl}(\mathrm{fl}(a-c)\mathrm{/}2)$, and $k$ by the helper's scaled norm on the ordered four-real-component list $(\delta ,0,\mathrm{\Re }z,\mathrm{\Im }z)$, corresponding to the complex vector $(\delta ,z)$, returning exact zero for the zero list. Then compute ${\mu }_{-}=\mathrm{fl}(h-k)$, ${\mu }_{+}=\mathrm{fl}(h+k)$, $q_{-}=\mathrm{\mid }\mathrm{fl}({\mu }_{-}-c)\mathrm{\mid }$, and $q_{+}=\mathrm{\mid }\mathrm{fl}({\mu }_{+}-c)\mathrm{\mid }$. Choose $\mu ={\mu }_{-}$ if $q_{-}\le q_{+}$ and $\mu ={\mu }_{+}$ otherwise. In exact arithmetic these are the two Hermitian eigenvalues and the stated nearest-eigenvalue/tie rule; in finite precision they specify the shift estimates without asserting their accuracy.

Next evaluate $\mathrm{\mid }{d}_{\mathrm{\ell }}\mathrm{\mid }$, $\mathrm{\mid }{e}_{\mathrm{\ell }-1}\mathrm{\mid }$, and $\mathrm{\mid }{e}_{\mathrm{\ell }}\mathrm{\mid }$ in that order and form $w=\mathrm{fl}(\tau \max (\mathrm{\mid }{e}_{\mathrm{\ell }-1}\mathrm{\mid },\mathrm{\mid }{e}_{\mathrm{\ell }}\mathrm{\mid }))$. Use the zero-shift branch exactly when $d_{\mathrm{\ell }}=0$ or $\mathrm{\mid }{d}_{\mathrm{\ell }}\mathrm{\mid }\le w$. In the ordinary branch form the first shift-vector entry as $\mathrm{fl}(T_{\mathrm{\ell }\mathrm{\ell }}-\mu )$ and use the stored lower Gram entry as its second entry. If any of these primitives fails the helper's normal-range/finite-value requirements, return the existing range-failure record immediately, before any sweep transformation for that block. This rule permits range failure in shift formation and makes no unconditional progress claim. For $p=0$ or only singleton active blocks none of the block stages is requested.

In the ordinary branch choose the first right Givens factor $G_{\mathrm{\ell }}$ so that $G_{\mathrm{\ell }}^{\ast }(T_{\mathrm{\ell }\mathrm{\ell }}-\mu ,T_{\mathrm{\ell }+1,\mathrm{\ell }}{)}^T=(\rho ,0{)}^T$, where $T_{\mathrm{\ell }+1,\mathrm{\ell }}=\overline{e_{\mathrm{\ell }}}{d}_{\mathrm{\ell }}$. In either branch, apply $B\leftarrow B{G}_{\mathrm{\ell }}$, update columns $\mathrm{\ell },\mathrm{\ell }+1$ of $V_{\mathrm{\Sigma }}$ by right multiplication by this same stored rotation, and append its exact unitary representative to the right list. For $j=\mathrm{\ell },\dots ,r-1$, choose a left Givens factor $H_j$ so that $H_j^{\ast }(B_{j,j},B_{j+1,j}{)}^T=(\rho ,0{)}^T$; apply it to rows $j,j+1$, update columns $j,j+1$ of $U_{\mathrm{\Sigma }}$ by right multiplication by this same stored rotation, and append its exact unitary representative to the left list. If $j<r-1$, choose a right Givens factor $G_{j+1}$ so that $(B_{j,j+1},B_{j,j+2})G_{j+1}=(\rho ,0)$; apply it to columns $j+1,j+2$, update the corresponding columns of $V_{\mathrm{\Sigma }}$ by right multiplication by this same stored rotation, and append its exact unitary representative to the right list. Continue until the bulge reaches the active block endpoint. A zero annihilation vector uses the identity rotation. In finite precision each Givens formation/update uses Scaled real and complex Householder and Givens update arithmetic, conjugate-transposed as needed for the stated left/right orientation; accumulator multiplications use the same stored local data but omit all target-zero writes. The dense accumulator multiplications follow the operationwise arithmetic model, while the ordered lists retain the exact unitary representatives used in analysis. Target zero writes to $B$ are included in their local error bounds. The initial shift-vector rotation does not zero an unrelated entry of $B$. Each Givens factor is chosen using Real Givens rotations and complex Givens transformations and its annihilation guarantee in Householder reflectors and Givens transformations are unitary and can annihilate prescribed entries.

After each processed block, increment $s$ by one. Once all current active blocks have been processed, start a new pass. If any operation, norm or pivot is invalid or nonfinite, or a nonzero exact intermediate is outside the normal range, return range-failure with the partial $B$, accumulators, transformation lists, deflation record, and available residual reports. On budget-exhausted, return the current $B$, accumulators, transformation lists, deflation record, triplet residuals, and orthogonality residuals. On converged, set $B_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}=B$ and return it with both dense accumulators, both transformation lists, the deflation record, triplet/orthogonality residuals, and its tall embedding ${\bar{B}}_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}=\left[\begin{array}{c}{B}_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}\\ 0_{M-p,p}\end{array}\right]$.

For every exit other than range-failure, let $B_{\mathrm{o}\mathrm{u}\mathrm{t}}$ be the current work matrix, and let $d_i=(B_{\mathrm{o}\mathrm{u}\mathrm{t}}{)}_{ii}$. Obtain each candidate value as ${\widehat{\sigma }}_i=\mathrm{\mid }{d}_i\mathrm{\mid }$ by the scaled two-component modulus evaluation (real absolute value is exact). On converged, the active-block condition makes $B_{\mathrm{o}\mathrm{u}\mathrm{t}}=B_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}$ diagonal; after budget-exhausted, the diagonal values are estimates only and carry no accuracy claim.

Let $W_B=\mathrm{diag}(U_{\mathrm{\Sigma }},I_{M-p})$. In exact arithmetic the tall-orientation accumulated factors are $U_C=U_B{W}_B$ and $V_C=V_B{V}_{\mathrm{\Sigma }}$. In finite precision form each dense product by row-by-column multiplication, with six-real-operation complex products and sequential componentwise rounded additions under the operationwise model; return range-failure if any required intermediate is invalid, nonfinite, or outside the normal range. Their exact theoretical factors are the products represented by the ordered local transformation lists together with the bidiagonal-reduction lists. If $m\ge n$, set $\widehat{U}=U_C$ and $\widehat{V}=V_C$; for each nonzero $d_i$, multiply column $i$ of $\widehat{U}$ by $d_i\mathrm{/}\mathrm{\mid }{d}_i\mathrm{\mid }$ in exact arithmetic, using componentwise division by the successful scaled modulus readout in finite precision. If $m<n$, set $\widehat{U}=V_C$ and $\widehat{V}=U_C$; in this wide orientation multiply column $i$ of $\widehat{V}$ by $d_i\mathrm{/}\mathrm{\mid }{d}_i\mathrm{\mid }$ instead. Sort the candidate values in nonincreasing order, breaking ties by their original diagonal index, and permute the corresponding columns of the returned triplet vectors in the same way; leave $B_{\mathrm{o}\mathrm{u}\mathrm{t}}$ and the transformation lists in their algorithmic order. For every candidate triplet $({\widehat{\sigma }}_i,{\widehat{u}}_i,{\widehat{v}}_i)$, report $\mathrm{\parallel }A{\widehat{v}}_i-{\widehat{\sigma }}_i{\widehat{u}}_i{\mathrm{\parallel }}_2$ and $\mathrm{\parallel }{A}^{\ast }{\widehat{u}}_i-{\widehat{\sigma }}_i{\widehat{v}}_i{\mathrm{\parallel }}_2$. Also report $\mathrm{\parallel }{\widehat{U}}^{\ast }\widehat{U}-I{\mathrm{\parallel }}_F$ and $\mathrm{\parallel }{\widehat{V}}^{\ast }\widehat{V}-I{\mathrm{\parallel }}_F$. The reports expose finite-budget and rounding behavior; they do not certify convergence by themselves. Matrix residuals use the Frobenius norm defined in The Frobenius norm $\parallel A{\parallel }_F=({\sum }_{i,j}\mathrm{\mid }{a}_{ij}{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}$ on real or complex matrices. Any additional sparsity cleanup is recorded separately in the deflation record with its discarded-entry perturbation.

For $A\in {\mathbb{F}}^{m\times n}$ over $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$, the exact-arithmetic algorithm Householder bidiagonalization and its finite precision status returns unitary $U,V$ and an upper bidiagonal $B\in {\mathbb{F}}^{p\times p}$, $p=\min (m,n)$, such that if $m\ge n$ then $U^{\ast }AV=\left[\begin{array}{c}B\\ 0\end{array}\right]$, and if $m<n$ then $U^{\ast }{A}^{\ast }V=\left[\begin{array}{c}B\\ 0\end{array}\right]$.

Let $A\in {\mathbb{F}}^{m\times n}$ have full column rank, $m\ge n$, and let $A\mathrm{\Pi }=QR$ be a full QR factorization, with $Q\in {\mathbb{F}}^{m\times m}$ unitary and $R\in {\mathbb{F}}^{m\times n}$ upper trapezoidal. Partition $R=\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right]$ after the first $k$ rows and columns, where $R_{11}\in {\mathbb{F}}^{k\times k}$ is nonsingular and $1\le k<n$. Write $X=R_{11}^{-1}{R}_{12}$, ${\omega }_i=\mathrm{\parallel }{e}_i^{\ast }{R}_{11}^{-1}{\mathrm{\parallel }}_2$, and ${\gamma }_j=\mathrm{\parallel }{R}_{22}(:,j){\mathrm{\parallel }}_2$. Use local one-based QR positions and define the original column labels $\pi (s)$ by $(A\mathrm{\Pi })(:,s)=A(:,\pi (s))$. Let $J=\{\pi (1),\dots ,\pi (k)\}$ be the selected set of $k$ columns and define its volume by $\mathrm{vol}(A(:,J))=\det (A(:,J{)}^{\ast }A(:,J){)}^{1\mathrm{/}2}$. For $1\le i\le k$ and $1\le j\le n-k$, replacing selected column $\pi (i)$ by trailing column $\pi (k+j)$ changes this volume by the factor $\frac{\mathrm{vol}(A(:,(J\setminus \{\pi (i)\})\cup \{\pi (k+j)\}))}{\mathrm{vol}(A(:,J))}=(\mathrm{\mid }{X}_{ij}{\mathrm{\mid }}^2+({\omega }_i{\gamma }_j{)}^2{)}^{1\mathrm{/}2}\mathrm{.}$

Let $n\ge 1$, and let $B,\tilde{B}\in {\mathbb{F}}^{n\times n}$ be upper bidiagonal, $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$, with $0\le ϵ<1$. Suppose each corresponding nonzero diagonal or superdiagonal entry is multiplied by a scalar whose modulus lies in $[1-ϵ,1+ϵ]$, and every zero entry remains zero. Put $q=(1+ϵ)\mathrm{/}(1-ϵ)$. Then $B$ and $\tilde{B}$ have the same number of zero singular values, and for each positive singular value in descending order, $q^{-2n}{\sigma }_i(B)\le {\sigma }_i(\tilde{B})\le q^{2n}{\sigma }_i(B)\mathrm{.}$

This is a relative perturbation theorem. A computed bidiagonal SVD inherits it only when a separate backward-error analysis verifies this componentwise relative assumption; an absolute normwise perturbation estimate by itself is not that verification.