Finite budget implicit shifted QR for bidiagonal singular values

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.