Arithmetical enhancements of the Kogbetliantz method for the SVD of order two

The singular value decomposition (SVD) of a square matrix of order two is a widely used numerical tool. In LAPACK [1] alone, its xLASV2 routine for the SVD of real upper triangular $2\times 2$ matrices is a building block for the QZ algorithm [2] for the generalized eigenvalue problem $Ax=\lambda Bx$, with $A$ and $B$ real and square, and for the SVD of real bidiagonal matrices by the implicit QR algorithm [3]. Also, the oldest method for the SVD of square matrices that is still in use was developed by Kogbetliantz [4], based on the SVD of order two, and as such is the primary motivation for this research.

This work explores how to compute the SVD of a general matrix of order two indirectly, by a careful scaling, a highly relatively accurate triangularization if the matrix indeed contains no zeros, and an alternative triangular SVD method, since the straightforward formulas for general matrices are challenging to be evaluated stably.

Let $G$ be a square real matrix of order $n$. The SVD of $G$ is a decomposition $G=U\Sigma V^{T}$, where $U$ and $V$ are orthogonal¹¹1If $G$ is complex, $U$ and $V$ are unitary and $G=U\Sigma V^{\ast}$, but this case is only briefly dealt with here. matrices of order $n$ of the left and the right singular vectors of $G$, respectively, and $\Sigma=\mathop{\mathrm{diag}}(\sigma_{1},\ldots,\sigma_{n})$ is a diagonal matrix of its singular values, such that $\sigma_{i}\geq\sigma_{j}\geq 0$ for all $i$ and $j$ where $1\leq i<j\leq n$.

In the step $k$ of the Kogbetliantz SVD method, a pivot submatrix of order two (or several of them, not sharing indices with each other, if the method is parallel) is found according to the chosen pivot strategy in the iteration matrix $G_{k}$, its SVD is computed, and $U_{k}$, $V_{k}$, and $G_{k}$ are updated by the transformation matrices $\mathsf{U}_{k}$ and/or $\mathsf{V}_{k}$, leaving zeros in the off-diagonal positions $(j_{k},i_{k})$ and $(i_{k},j_{k})$ of $G_{k+1}$, as in

The left and the right singular vectors of the pivot matrix are embedded into identities to get $\mathsf{U}_{k}$ and $\mathsf{V}_{k}$, respectively, with the index mapping from matrices of order two to $\mathsf{U}_{k}$ and $\mathsf{V}_{k}$ being $(1,1)\mapsto(i_{k},i_{k})$, $(2,1)\mapsto(j_{k},i_{k})$, $(1,2)\mapsto(i_{k},j_{k})$, $(2,2)\mapsto(j_{k},j_{k})$, where $1\leq i_{k}<j_{k}\leq n$ are the pivot indices. The process is repeated until convergence, i.e., until for some $k=K$ the off-diagonal norm of $G_{K}$ falls below a certain threshold.

If $G$ has $m$ rows and $n$ columns, $m>n$, it should be preprocessed [5] to a square matrix $G_{0}$ by a factorization of the URV [6] type (e.g., the QR factorization with column pivoting [7]). Then, $U_{0}^{T}GV_{0}=G_{0}$, where $G_{0}$ is triangular of order $n$ and $U_{0}$ is orthogonal of order $m$. If $m<n$, then the SVD of $G^{T}$ can be computed instead.

In all iterations it would be beneficial to have the pivot matrix $\widehat{G}_{k}$ triangular, since its SVD can be computed with high relative accuracy under mild assumptions [8]. This is however not possible with time consuming but simple, quadratically convergent [9, Remark 6] pivot strategy that chooses the pivot with the largest off-diagonal norm $\sqrt{|g_{j_{k}i_{k}}^{(k)}|^{2}+|g_{i_{k}j_{k}}^{(k)}|^{2}}$, but is possible, if $G_{0}$ is triangular, with certain sequential cyclic (i.e., periodic) strategies [5, 10] like the row-cyclic and column-cyclic, and even with some parallel ones, after further preprocessing $G_{0}$ into a suitable “butterfly” form [11].

Although the row-cyclic and column-cyclic strategies ensure global [10] and asymptotically quadratic [9, 12] convergence of the method, as well as its high relative accuracy [13], the method’s sequential variants remain slow on modern hardware, while preprocessing $G$ to $G_{0}$ (in the butterfly form or not) can only be partially parallelized.

This work is a part of a broader effort [14, 15] to investigate if a fast and accurate (in practice if not in theory) variant of the Kogbetliantz method could be developed, that would be entirely parallel and would function on general square matrices without expensive preprocessing, with full pivots $\widehat{G}_{k}^{[\ell]}$, $1\leq\ell\leq\mathsf{n}\leq\lfloor n/2\rfloor$, that are independently diagonalized, and with $\mathsf{n}$ ensuing concurrent updates of $U_{k}$, $\mathsf{n}$ of $V_{k}$, and $\mathsf{n}$ from each of the sides of $G_{k}$ in a parallel step. This way the employed parallel pivot strategy does not have to be cyclic. A promising candidate is the dynamic ordering [16, 17, 15].

The proposed Kogbetliantz SVD of order two supports a wider exponent range of the elements of a triangular input matrix for which both singular values are computed with high relative accuracy than xLAEV2, although the latter is slightly more accurate when comparison is possible. Matrices of the singular vectors obtained by the proposed method are noticeably more numerically orthogonal. With respect to [14, 15] and a general matrix $G$ of order two, the following enhancements have been implemented:

1. 1.

   The structure of $G$ is exploited to the utmost extent, so the triangularization and a stronger scaling is employed only when $G$ has no zeros, thus preserving accuracy.

2. 2.

   The triangularization of $G$ by a special URV factorization is tweaked so that high relative accuracy of each computed element is provable when no underflow occurs.

3. 3.

   The SVD procedure for triangular matrices utilizes the latest advances in computing the correctly rounded functions, so the pertinent formulas from [14, 15] are altered.

4. 4.

   The left singular vectors are computed by a heuristic when the triangularization is involved, by composing the two plane rotations—one from the URV factorization, and the other from the triangular SVD—into one without the matrix multiplication.

High relative accuracy of the singular values of $G$ is observed, but not proved, when the range of the elements of $G$ is narrower than about half of the entire normal range.

This paper is organized as follows. In Section 2 the floating-point environment and the required operations are described, and some auxiliary results regarding them are proved. Section 3 presents the proposed SVD method. In Section 4 the numerical testing results are shown. Section 5 concludes the paper with the comments on future work.

Let $x$ be a real, infinite, or undefined (Not-a-Number) value: $x\in\mathbb{R}\cup\{-\infty,+\infty,\mathrm{NaN}\}$. Its floating-point representation is denoted by $\mathop{\mathrm{fl}}(x)$ and is obtained by rounding $x$ to a value of the chosen standard floating-point datatype using the rounding mode in effect, that is here assumed to be to nearest (ties to even). If the result is normal, $\mathop{\mathrm{fl}}(x)=x(1+\epsilon)$, where $|\epsilon|\leq\varepsilon=2^{-p}$ and $p$ is the number of bits in the significand of a floating-point value. In the LAPACK’s terms, $\varepsilon_{\text{{x}}}=\text{{xLAMCH('E')}}$, where $\text{{x}}=\text{{S}}$ or D. Thus, $p_{\text{{x}}}=24$ or $53$, and $\varepsilon_{\text{{x}}}=2^{-24}$ or $2^{-53}$ for single (S) or double (D) precision, respectively. The gradual underflow mode, allowing subnormal inputs and outputs, has to be enabled (e.g., on Intel-compatible architectures the Denormals-Are-Zero and Flush-To-Zero modes have to be turned off). Trapping on floating-point exceptions has to be disabled (what is the default non-stop handling from [18, Sect. 7.1]).

Possible discretization errors in input data are ignored. Input matrices in floating-point are thus considered exact, and are, for simplicity, required to have finite elements.

The Fused Multiply-and-Add (FMA) function, $\mathop{\mathrm{fma}}(a,b,c)=\mathop{\mathrm{fl}}(a\cdot b+c)$, is required. Conceptually, the exact value of $a\cdot b+c$ is correctly rounded. Also, the hypotenuse function, $\mathop{\mathrm{hypot}}(a,b)=\mathop{\mathrm{fl}}(\sqrt{a^{2}+b^{2}})$, is assumed to be correctly rounded (unless stated otherwise), as recommended by the standard [18, Sect. 9.2], but unlike²²2https://members.loria.fr/PZimmermann/papers/accuracy.pdf many current implementations of the routines hypotf and hypot. Such a function (see also [19]) never overflows when the rounded result should be finite, it is zero if and only if $|a|=|b|=0$, and is symmetric and monotone. The CORE-MATH library [20] provides an open-source implementation³³3https://core-math.gitlabpages.inria.fr of some of the optional correctly rounded single and double precision mathematical C functions (e.g., cr_hypotf and cr_hypot).

Radix two is assumed for floating-point values. Scaling of a value $x$ by $2^{s}$ where $s\in\mathbb{Z}$ is exact if the result is normal. Only non-normal results can lose precision. Let, for $x=\pm 0$, $e_{x}=0$ and $f_{x}=0$, and for a finite non-zero $x$ let the exponent be $e_{x}=\lfloor\lg|x|\rfloor$ and the “mantissa” $1\leq|f_{x}|<2$, such that $x=\mathop{\mathrm{fl}}(2^{e_{x}}f_{x})$. Also, let $f_{x}=x$ for $x=\pm\infty$, while $e_{x}=0$. Note that $f_{x}$ is normal even for subnormal $x$. Keep in mind that the frexp routine represents a finite non-zero $x$ with $e_{x}^{\prime}=e_{x}+1$ and $f_{x}^{\prime}=f_{x}/2$.

Let $\mu$ be the smallest and $\nu$ the largest positive finite normal value. Then, in the notation just introduced, $e_{\mu}=\lg\mu=-126$ or $-1022$, and $e_{\nu}=\lfloor\lg\nu\rfloor=127$ or $1023$, for single or double precision. Lemma 2.1 can now be stated using this notation.

The claims of Lemma 2.1 and its following corollaries were used and their proofs partially sketched in [21], e.g. They are expanded and clarified here for completeness.

An underlined expression denotes a computed floating-point approximation of the exact value of that expression. Given $\tan(2\phi)$, for $0\leq|\phi|\leq\pi/4$, $\tan\phi$ and $\underline{\tan\phi}$ are

if $\tan(2\phi)$ and $\mathop{\mathrm{fl}}(\tan(2\phi))$ are finite, or $\mathop{\mathrm{sign}}(\tan(2\phi))$ and $\mathop{\mathrm{sign}}(\mathop{\mathrm{fl}}(\tan(2\phi)))$ otherwise.

For any $w\in\mathbb{R}$, let $\mathbf{w}=(e_{w},f_{w})=2^{e_{w}}f_{w}$ and $\mathop{\mathrm{fl}}(\mathbf{w})=\mathop{\mathrm{fl}}(2^{e_{w}}f_{w})\approx w$. Even $w$ such that $|w|>\nu$ or $0<|w|<\check{\mu}$, where $\check{\mu}$ is the smallest positive non-zero floating-point value, can be represented with a finite $e_{w}$ and a normalized $f_{w}$, though $\mathop{\mathrm{fl}}(\mathbf{w})$ is not finite or non-zero, respectively. The closest double precision approximation of $\mathbf{w}$ is $\underline{w}=\text{{scalbn($f_{w},e_{w}$)}}$, with a possible underflow or overflow, and similarly for single precision (using scalbnf). A similar definition could be made with $e_{w}^{\prime}$ and $f_{w}^{\prime}$ instead.

An overflow-avoiding addition and an underflow-avoiding subtraction of positive finite values $x$ and $y$, resulting in such exponent-“mantissa” pairs, can be defined as

and, assuming $x\neq y$ (for $x=y$ let $x\ominus y=(0,0)$ directly),

where $c=e_{\mu}+p-1-\min\{e_{x},e_{y}\}$. In (4), $z\leq\nu$ in both cases, since $\mathop{\mathrm{fl}}(\nu/2+\nu/2)=\nu$.

Several arithmetic operations on $(e,f)$-pairs can be defined (see also [22]), such as

which are unary operations. The binary multiplication and division are defined as

and the relation $\prec$, that compares the represented values in the $<$ sense, is given as

Let, for any $G$ of order $n$, where INT_MIN is the smallest representable integer,

A prescaling of $G$ as $\underline{G^{\prime}}=2^{s}G$, that avoids overflows, and underflows if possible, in the course of computing the SVD of $\underline{G^{\prime}}$ (and thus of $\underline{G}\approx 2^{-s}\underline{G^{\prime}}$), is defined by $s$ such that

where $\mathfrak{s}=0$ for $n=1$. For $n=2$, $\mathfrak{s}$ is chosen such that certain intermediate results while computing the SVD of $\underline{G^{\prime}}$ cannot overflow, as explained in [14, 15] and Section 3, but the final singular values are represented in the $(e,f)$ form, and are immune from overflow and underflow as long as they are not converted to simple floating-point values. If $s\geq 0$, the result of such a prescaling is exact. Otherwise, some elements of $\underline{G^{\prime}}$ might be computed inexactly due to underflow. If for the elements of $G$ holds

then $s\geq 0$, and $g_{ij}^{\prime}=0$ or $\mu\leq|g_{ij}^{\prime}|\leq\nu$, i.e., the elements of $G^{\prime}$ are zero or normal.

This section presents a Kogbetliantz-like procedure for computing the singular values of $G$ when $n=2$, and the matrices of the left ($U$) and the right ($V$) singular vectors.

In general, $U$ is a product of permutations (denoted by $P$ with subscripts and including $I$), sign matrices (denoted by $S$ with subscripts) with each diagonal element being either $1$ or $-1$ while the rest are zeros, and plane rotations by the angles $\vartheta$ and $\varphi$. If $U_{\vartheta}$ is not generated, the notation changes from $U_{\varphi}$ to $U_{\phi}$. Likewise, $V$ is a product of permutations, a sign matrix, and a plane rotation by the angle $\psi$, where

Depending on its pattern of zeros, a matrix of order two falls into one of the 16 types $\mathfrak{t}$ shown in (13), where $\circ=0$ and $\bullet\neq 0$. Some types are permutationally equivalent to others, what is denoted by $\mathfrak{t}_{1}\cong\mathfrak{t}_{2}$, and means that a $\mathfrak{t}_{1}$-matrix can be pre-multiplied and/or post-multiplied by permutations to be transformed into a $\mathfrak{t}_{2}$-matrix, and vice versa, keeping the number of zeros intact. Each $\mathfrak{t}\neq 0$ has its associated scale type $\mathfrak{s}$.

For $\mathfrak{s}=0$, there is one equivalence class of matrix types, represented by $\mathfrak{t}=9$. For $\mathfrak{s}=1$, there are three classes, represented by $\mathfrak{t}=3$, $\mathfrak{t}=5$, and $\mathfrak{t}=13$, while for $\mathfrak{s}=2$ there is one class, $\mathfrak{t}=15$. The SVD computation for the first three classes is straightforward, while for the fourth and the fifth class is more involved. A matrix of any type, except $\mathfrak{t}=15$, can be permuted into an upper triangular one. If a matrix so obtained is well scaled, its SVD can alternatively be computed by xLASV2. However, xLASV2 does not accept general matrices (i.e., $\mathfrak{t}=15$), unlike the proposed method, which is a modification of [15] when $J=I$, and consists of the following three phases:

1. 1.

   For $G$ determine $\mathfrak{t}$, $\mathfrak{s}$, and $s$ to obtain $\underline{G^{\prime}}$. Handle the simple cases of $\mathfrak{t}$ separately.

2. 2.

   If $\mathfrak{t}\cong 13$ or $\mathfrak{t}=15$, factorize $\underline{G^{\prime}}$ as $U_{+}RV_{+}$, such that $U_{+}$ and $V_{+}$ are orthogonal, and $R$ is upper triangular, with $\min\{r_{11},r_{12},r_{22}\}>0$ and all $r_{ij}$ finite, $1\leq i\leq j\leq 2$.

3. 3.

   From the SVD of $\underline{R}$ assemble the SVD of $\underline{G^{\prime}}$. Optionally backscale $\underline{\Sigma^{\prime}}$ by $2^{-s}$.

The phases 1, 2, and 3 are described in Sections 3.1, 3.2, and 3.3, respectively.