Householder bidiagonalization and its finite precision status

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.