数理解析研究所講究録
第 1606 巻 2008 年 87-104

87

Recent Developments ofthe mdLVs Algorithm for Singular Values
and the I-SVD Algorithm
for Singular Value Decomposition
特異値計算の mdLVs アルゴリズムと特異値分解の I-SVD アルゴリズムにおける最近の進展
京都大学情報学研究科 中村佳正 1 (Yoshimasa Nakamura)
Graduate School of Informatics, Kyoto University
京都府立大学人間環境学部 岩崎雅史 (Masashi Iwasaki)
Faculty of Human Environment, Kyoto Prefectural University
京都大学情報学研究科 木村欣司 (Kinji Kimura)
Graduate School of Informatics, Kyoto University
奈良女子大学人間文化研究科 高田雅美 (Masami Takata)
Graduate School of Humanitics & Sciences, Nara Women’s University

Abstract Recently a new algorithm for computing singular values named the mdLVs (modified
discrete Lotka-Volterra with shift) is designed. The first part of this report is a brief survey ofthe recent developments on the positivity and shifts of the mdLVs algorithm. The second partis an exposition of the I-SVD (Integrable Singular Value Decomposition) algorithm which is acombination of the mdLVs algorithm and the dLV-type transformation for computing singular
vectors of bidiagonal matrices. Because of the separation of computation of singular values
from that of singular vectors the I-SVD algorithm mns in $O(m^{2})$ fiops and is rather faster than
DBDSQR code of LAPACK.

1. Introduction
An algorithm named the $dLV$ (discrete Lotka-Volterra) algorithm for computing bidiagonal
matrix singular values has been discussed in the series ofpapers [11, 12, 13, 14, 15], where suchbidiagonal matrices can be derived from arbitrary nonsingular matrices through the Householder
preconditioning process [91. A general background is fumished in [21, 22]. See also a recent
review paper [3]. The recurrence relation of $dLV$ itself is a discrete-time integrable dynamicalsystem. Convergence of the $dLV$ algorithm to singular values is shown in [11]. See \S 4 ofthis report. The basic fact is that $dLV$ is a deformation equation of orthogonal polynomials(OPs). The parameter of $dLV$ should be positive. Therefore the recurrence relation of $dLV$ islynaka@i.kyoto-u.ac.jp

88
subtraction free. A positivity of Hankel determinants is ensured whose elements are moments
associated with OPs and then all the variables of $dLV$ are also positive [14]. This fact will be
reviewed in \S 3. In a recent work [151 an exponential stability of $dLV$ in a local sense is provedby using the existence of a center manifold around the fixed points. This implies that $dLV$ isrobust and highly credible algorithm. The convergence rate of the $dLV$ algorithm is linear [11].
Therefore some shifted versions of $dLV$ have been formulated in [12, 13].
The mdLVs (modified $dLV$ with shift) algorithm presented in [131 is shown to satisfy thesame positivity of variables and has a higher relative accuracy. Speed of the mdLVs algorithm
depends on the choice of shift of origin. The cost for computing shifts wastes more than 30% of
total execution time of the mdLVs code [291, where the shift is determined as a lower bound of
the ninimal singular value of given upper bidiagonal matrices. Therefore a lighter shift basedon more accurate bound must be important to accelerate the mdLVs algorithm. The Johnson
bound [17], a Gersgorin-type lower bound of symmetric tridiagonal, has been adopted in the
mdLVs code [29]. Recently a new lower bound is found which is called the p-th generalized
Newton bound. In the first part of this report (\S 2 \S 5) we discuss the recent developments on
the positivity and shifts of the mdLVs algorithm.
The second part (\S 6 \S 7) is an exposition of the I-SVD (Integrable Singular Value Decomposition) algorithm[16] which is a combination of the mdLVs algorithm for singular values andthe dLV-type transformation for singular vectors of $m\cross m$ bidiagonal matrices. Namely, computations of singular values and vectors are completely separated in I-SVD. Here the dLV-typetransformation performs an accurate double Cholesky decomposition of a shifted symmetrictridiagonal matrix which gives rise to a twisted factorization of the same matrix. Each singular vector is computed from the twisted matrix within $0(m)$ flops. Then the I-SVD algorithm
solves the bidiagonal SVD problem within $0(m^{2})$ flops. It is shown in [301 for some class oftest matrices that the I-SVD code is rather faster than the standard SVD code of LAPACK.
$\sim$$\sim$

2. Orthogonal Polynomials: Preliminary$(k=1,2, \ldots)$ be a sequences of real
Let us begin with Favard’s theorem [21. Let
numbers. When the bilinear form
is positive for any $m$ , theniscalled positive. It is known thatis positive if and only if the Hankel determinants$\{s_{k}\},$$\sum_{k=0}^{m}\sum_{\ell=0}^{m}s_{k+\ell}\alpha_{k}\alpha_{l}$$\{s_{k}\}$$\{s_{k}\}$$s_{0}$$s_{1}$$s_{1}$$s_{2}$:.

$D_{n+1}:=$:$s_{n}$

are positive for any $n=0,1,$
Theorem (Favard) Let$\{a_{k}\},$$\ldots$$s_{n+1}$

.. .

...
...

$s_{n}$$S_{n+1}$:.

,$(n=0,1,2, \ldots)$$s_{2n}$.

$\{b_{k}\},$$(k=1,1, \ldots)$ be sequences of real numbers. Let $\{p_{k}(\lambda)\}$89
be polynomials of defined by the three terms recuirence relation$\lambda$$Pk+1(\lambda)=(\lambda-b_{k}+i)_{Pk}(\lambda)-a_{k^{2}Pk-1}(\lambda)$$p_{0}(\lambda)=1$,$p_{1}(\lambda)=\lambda-b_{1}$Then there exists a unique linear functional$k,$ $l=0,1,$
for any positive constant$:=J[\lambda^{k}],(k=0,1, \ldots)$ are positive.
$\ldots)$$\ell,$,.

such that $s_{0}=J[1],J\lceil flk(\lambda)p\ell(\lambda)]=0,(k\neq$. Moreover, $a_{k^{2}}>0$ if and only if the moments$J$$s_{0}$$s_{k}$The polynomial$pk(\lambda),$$(k=1,2, \ldots)$ takes the determinant form [27]$Pk( \lambda)=\frac{1}{D_{k}}$$s_{0}$$s_{1}$$s_{1}$$s_{2}$::$s_{k-1}s_{k}$1Then the coefficients$a_{k^{2}}$$\lambda$...

$s_{k}$

...
...

...

$s_{k+1}$:.
$s_{2k-1}$$\lambda^{k}$of the recurrence relation are$a_{k^{2}}= \frac{D_{k-1}D_{k+1}}{D_{k^{2}}}$.

It is shown from $a_{1^{2}}\cdots a_{k^{2}}=D_{k+1}/D_{k}$ that $D_{k}>0$ for any$k$and the corresponding moments

are positive.
Favard’s theorem says that the polynomials
defined by the three terms recurrence
relation with positive coefficientsare orthogonal with respect to the linear functional .
Namely,. In this case the corresponding set of moments
is positive and vice verse.
OPs have some special features. One of them is the position of zeros. It is known that zerosof OPs are mutually distinct real numbers and has an interlacing property[l]. For example,
let $\lambda_{j}^{(n-1)},(j=1,2, \ldots, n-1)$ and $\lambda_{j}^{(n)},(j=1,2, \ldots, n)$ be zeros of,and$\{pk(\lambda)\}$$a_{k^{2}}$$J$$J[pk(\lambda)p\ell(\lambda)]=s_{0}a_{1^{2}}\cdots a_{k^{2}}\delta_{k,l}$$\{s_{k}\}$$p_{n-1}(\lambda)$$p_{n}(\lambda)$

respectively. Then$\lambda_{1}^{(n)}<\lambda_{1}^{(n-1)}<\lambda_{2}^{(n)}<\lambda_{2}^{(n-1)}<\cdots<\lambda_{n-1}^{(n-1)}<\lambda_{n}^{(n)}$This leads to the following statement. The rational function
a partial fraction expansion,$\frac{p_{n-1}(\lambda)}{p_{n}(\lambda)}=\sum_{j=1}^{n}\frac{\rho_{j}^{(n)}}{\lambda-\lambda_{j}^{(n)}}$.

$p_{n-1}(\lambda)/p_{n}(\lambda)$$\rho_{j}^{(n)}=\frac{p_{n-1}(\lambda_{j}^{(n)})}{p_{n}’(\lambda_{j}^{(n)})}$of degree$n$admits.

, the Christoffel coefficients, satisfy
From the interlacing property it follows that the residues.
the positivity condition
For the Hermite, Legendre and Chebyshev polynomials every moments with odd order are
zero, namely, $s_{2k-1}=0$ . In the linear functionals of those cases the measure
and the$\rho_{j}^{(n)}$$\rho_{j}^{(n)}>0$$d\mu(\lambda)$

90
are invariant under the exchange. The linear functional
satisfying$s_{2k-1}=J[\lambda^{2k-1}]=0$ ,$(k=1,2, \ldots)$ is called symmetric and the corresponding OP iscalled symmetric. When $d\mu(\lambda)=w(\lambda)d\lambda$ , the weight function
is an even function overthe interval $(-\xi, \xi)$ . The coefficients of the recurrence relation are zero for symmetric OPs.
Namely, $b_{k}=0,$ $(k=1,2, \ldots)$ . In the following discussions we restrict ourselves to symmetriccontour$C_{\lambda}$$\lambdaarrow-\lambda$$J$$w(\lambda)$$b_{k}$OPs.
Let us consider the three terms recurrence relation of symmetric OPs$pk+1(\lambda)=\lambda_{Pk}(\lambda)-a_{k^{2}}p_{k-1}(\lambda)$$p_{0}(\lambda)=1$,$p_{1}(\lambda)=\lambda$,.

In [14] we obtain the Christoffel-Darboux formula for symmetric OPs as follows.

$\{\begin{array}{ll}a_{1^{2}}\cdots a_{2m-1^{2}}(\sum_{j=1}^{m}\frac{p_{2j-1}(.\lambda)p_{2j-1}(\kappa)}{a_{1^{2}}\cdot\cdot a_{2j-1^{2}}}) =\frac{p_{2m-1}(\lambda)p_{2m+1}(\kappa)-p_{2m+1}(\lambda)p_{2m-1}(\kappa)}{\kappa^{2}-\lambda^{2}} for k=2m-1a_{1^{2}}\cdots a_{2m}^{2}(\sum_{j=1}^{m}\frac{p_{2j}(\lambda)p_{2j}(\kappa)}{a_{1^{2}}\cdots a_{2j^{2}}}+p_{0}(\lambda)p_{0}(\kappa)) =\frac{p_{2m}(\lambda)p_{2m+2}(\kappa)-p_{2m+2}(\lambda)p_{2m}(\kappa)}{\kappa^{2}-\lambda^{2}} for k=2m\end{array}$

In contrast to the case of usual OPs [1, 2, 271, a parity emerges. The Christoffel-Darboux
formula is useful, for example, to discuss the convergence of OP series.

3 Discrete Lotka-Volterra System and Its Positivity
In this section we fist define a kemel polynomial
To this end we assume $pk(\kappa)\neq 0$ .

$p_{k}^{*}(\lambda)$of the original symmetric OP$p_{k}(\lambda)$.

$p_{k}^{*}(\lambda);=\{\begin{array}{ll}-\frac{a_{1^{2}}\cdots a_{2m-1^{2}}}{p_{2m-1}(\kappa)}\sum_{j=1}^{m}\frac{p_{2j-1}(.\lambda).p_{2j-1}(\kappa)}{a_{1^{2}}\cdot a_{2j-1^{2}}} for k=2m-1-\frac{a_{1^{2}}\cdots a_{2m}^{2}}{p_{2m}(\kappa)}(\sum_{j=1}^{m}\frac{p_{2j}(\lambda)p_{2j}(\kappa)}{a_{1^{2}}\cdots a_{2j^{2}}}+p_{0}(\lambda)p_{0}(\kappa)) for k=2m\end{array}$

Then the Christoffel-Darboux formula leads to$p_{k}^{*}( \lambda)=\frac{1}{\kappa^{2}-\lambda^{2}}(Pk+2(\lambda)+A_{kPk}(\lambda))$,$A_{k}$$.=- \frac{p_{k+2}(\kappa)}{pk(\kappa)}$.

is an odd function. When $k=2m,$
When $k=2m-1,$
is even. The poles
is a polynomial of degree . The transformation
are apparent poles. Hence$p_{k}(\lambda)$$p_{k}(\lambda)$$k$$p_{k}^{*}(\lambda)$$\{Pk(\lambda)\}arrow\{p_{k}^{*}(\lambda)\}$$\lambda=\pm\kappa$

91
is just the Christoffel transformation for the symmetric OP
linear functionalby$\{p_{k}(\lambda)\}$. Let us introduce a new$J^{*}$$J^{*}[A(\lambda)]:=J[(\kappa^{2}-\lambda^{2})A(\lambda)]$for any polynomialisfunctional.
$w^{*}(\lambda)$$A(\lambda)$,

and a suitable constant $\kappa<0$ . The corresponding weight function. We can generalize a theorem in $[2|$ on the positivity of linear$:=(\kappa^{2}-\lambda^{2})w(\lambda)$

Theorem [14] Let the linear functional$\xi>0$ . The
is positive, i.e.
$J^{*}$$J$$\{s_{k}^{*}:=J^{*}[\lambda^{k}]$be positive definite over the interval $[-\xi,$is positive over$[-\xi,\xi, ]$$\xi,$$|$

with, if and only if $\kappa\leq-\xi$ .

We now consider a successive use of the Christoffel transformations$p_{k}^{(n+1)}= \frac{1}{(\kappa^{(n)})^{2}-\lambda^{2}}(p_{k+2}^{(n)}+A_{k}^{(n)}p_{k}^{(n)})$,$A_{k}^{(n)}:=- \frac{p_{k+2}^{(n)}(\kappa^{(n)})}{p_{k}^{(n)}(\kappa^{(n)})}$,$(n=0,1, \ldots)$

to generate a sequence of kemel polynomials$\{p_{k}^{(0)}.=p_{k}(\lambda)\}arrow\{p_{k}^{(1)}.=p_{k}^{*}(\lambda)\}arrow\{p_{k}^{(2)}\}arrow\cdots$

,

where
follows from.
As the compatibility conditions of the Christoffel transformation and the recurrence relation
we obtain$p_{k}^{(n)}(\kappa^{(n)})\neq 0$$\kappa^{(n)}>-\lambda_{1}^{(n)}$$p_{k+1}^{(n+1)}=\lambda p_{k}^{(n+1)}-(a_{k}^{(n+1)})^{2}p_{k-1}^{(n+1)}$$(a_{k}^{(n+1)})^{2}$$=$$(a_{k}^{(n)})^{2} \frac{A_{k}^{(n)}}{A_{k-1}^{(n)}}$$=$$(a_{k}^{(n)})^{2} \frac{p_{k+2}^{(n)}(\kappa^{(n)})p_{k-1}^{(n)}(\kappa^{(n)})}{p_{k}^{(n)}(\kappa^{(n)})p_{k+1}^{(n)}(\kappa^{(n)})}$.

Let us set$\hat{u}_{k}^{(n)}:=(a_{k}^{(n)})^{2}\frac{p_{k-1}^{(n)}(\kappa^{(n)})}{p_{k}^{(n)}(\kappa^{(n)})}$

It follows from
that
. Let
that in the partial fraction expansion$p_{-1}^{(n)}=0$$\hat{u}_{0}^{(n)}=0$$\lambda_{j}^{(n)},$.

$(j=1, \ldots, k)$ be

zeros of the OP$p_{k}^{(n)}(\lambda)$. Note$\frac{p_{k-1}^{(n)}(\kappa^{(n)})}{p_{k}^{(n)}(\kappa^{(n)})}=\sum_{j=1}^{k}\frac{\rho_{j}^{(n)}}{\kappa^{(n)}-\lambda_{j}^{(n)}}$the residues
that
Inserting$\rho_{j}^{(n)}$of are positive. While it follows from the positivity of the linear functional. Thusand then.
into the three terms recurrence relation we derive$\kappa^{(n)}-\lambda_{j}^{(n)}<0$$\hat{u}_{k}^{(n)}$$p_{k-1}^{(n)}(\kappa^{(n)})/p_{k}^{(n)}(\kappa^{(n)})<0$$(a_{k+1}^{(n)})^{2}=\hat{u}_{k+1}^{(n)}(\kappa^{(n)}+\hat{u}_{k}^{(n)})$$\hat{u}_{k}^{(n)}<0$.

$J^{*}$

92
Similarly we have.

$(a_{k}^{(n+1)})^{2}=\hat{u}_{k}^{(n)}(\kappa^{(n)}+\hat{u}_{k+1}^{(n)})$We elimuinate$(a_{k}^{(n+1)})^{2}$to have$\hat{u}_{k}^{(n+1)}(\kappa^{(n+1)}+\hat{u}_{k-1}^{(n+1)})=\hat{u}_{k}^{(n)}(\kappa^{(n)}+\hat{u}_{k+1}^{(n)})$.

This equation was first derived by Hirota (1997, [101) and $Sp\ddot{m}donov$ -Zhedanov(1997, [26]),
are arbitrary constant. In our caseindependently, whereshould be less than or equal
functional,toto guarantee the positivity of the linear
sayand then Hankel determi. Sinceis expressed as a ratio of the Hankel determinants, this property is very
nants
important to design stable numerical algorithm. Define$\kappa^{(n)}$$\kappa^{(n)}$$J^{(n)}$$-\xi$$D_{k}^{(n)}$$\hat{u}_{k}^{(n)}$$\delta^{(n)}:=\frac{1}{(\kappa^{(n)})^{2}}>0$,$u_{k}^{(n)}=\kappa^{(n)}\hat{u}_{k}^{(n)}>0$.

We can introduce a scale change
to relax the condition $0<\delta^{(n)}\leq 1/\xi^{2}$
to $0<\delta^{(n)}\leq M$ for some positive constant $M$ . Then we obtain$u_{k}^{(n)}arrow 1/(\xi^{2}M)u_{k}^{(n)}$$u_{k}^{(n+1)}(1+\delta^{(n+1)}u_{k-1}^{(n+1)})=u_{k}^{(n)}(1+\delta^{(n)}u_{k+1}^{(n)})$$u_{k}^{(n)}>0$,$0<\delta^{(n)}\leq M$,,$(n=0,1, \ldots, k=1,2, \ldots)$ .

as the value of $u_{k}=u_{k}(t)$ at the time
Let us regard. Keeping to a$\delta^{(n)}arrow+0$such that
constant we take a limit. We then derive the system ofdifferential equations$u_{k}^{(n)}$$t= \sum_{j=0}^{n-1}\delta^{(j)}$$t$$\delta^{(n+1)}/\delta^{(n)}arrow 1$$\frac{du_{k}}{dt}=u_{k}(u_{k+1}-u_{k-1})$,$u_{0}(t)=0$,$(k=1,2, \ldots)$for the variable $u_{k}=u_{k}(t)$ from the recurrence relation. This process corresponds to the limit
and does not violate the positivity of linear functionals. This system is sometimes
called the Lotka-Volterra (LV) system in mathematical biology. In this section it is shownthat the successive Christoffel transformations of symmetric OPs induce a deformation of theof the three terms recurrence relation. The resulting deformation equation iscoefficients$dLV$
system having a positive explicit solution.
the$\kappa^{(n)}arrow-\infty$$\{a_{k}^{(n)}\}$

4

Convergence of $dLV$ Algorithm

In the series of papers [11, 12, 151 it is shown a solution of $dLV$ converges to the same limit. Namely,as of the LV for any choice of positive$\delta^{(n)}$$\lim_{narrow\infty}u_{2k-1}^{(n)}=\sigma_{k^{2}}$where$\sigma_{k}$are singular vales of$B$,$\lim_{narrow\infty}u_{2k}^{(n)}=0$

such that$\sigma_{1}>\sigma_{2}>\cdots>\sigma_{m}>0$

.

,

93
It is to be remarked that the initial value setting is different from that in the LV case. We shouldchoose$u_{2k-1}^{(0)}= \frac{b_{2k-1^{2}}}{1+\delta(0)_{u_{2k-2}^{(0)}}}$$u_{2k}^{(0)}= \frac{b_{2k^{2}}}{1+\delta(0)u_{2k-1}^{(0)}}$,,$(k=1,2, \ldots, m)$ ,$(k=1,2, \ldots, m-1)$ ,and $u_{2m}^{(0)}=0$ . We named this procedure the $dLV$ algorithm for computing
singular values of bidiagonal matrices.
More important notion is the numerical stability. It is known [25], for example, the qd
algorithm is numerically instable because of division by a small amount. One the other hand
the Demmel-Kahan QR has been standard algorithm for a long time as a stable algorithm inspite of slow speed. The numerical stability of $dLV$ is proved in [12]. The starting point is amatrix representation of $dLV$.

as well as$u_{0}^{(0)}=0$$L^{(n+1)}R^{(n+1)}=R^{(n)}L^{(n)}-( \frac{1}{\delta(n)}-\frac{1}{\delta(n+1)})I$$L^{(n)}:=(\begin{array}{llll}J_{1}^{(n)} O1 J_{2}^{(n)} \ddots .

1 J_{m}^{(n)}\end{array})$,,$R^{(n)}:=(\begin{array}{llll}1 V_{1}^{(n)} 1 \ddots \ddots V_{m-1}^{(n)}O 1\end{array})$$J_{k}^{(n)}:= \frac{1}{\delta(n)}(1+\delta^{(n)}u_{2k-2}^{(n)})(1+\delta^{(n)}u_{2k-1}^{(n)})$,$V_{k}:=\delta^{(n)}u_{2k-1}^{(n)}u_{2k}^{(n)}$is the $m\cross m$ unit matrix. Note thatmatrix. Let us introduce new nonnegative variables

where$1/\delta^{(n)}-1/\delta^{(n+1)}$$I$$R^{(n)}L^{(n)}$$w_{k}^{(n)}$,,gives a shift of origin for thedefined as$w_{k}^{(n)}:=u_{k}^{(n)}(1+\delta^{(n)}u_{k-1}^{(n)})$and a tridiagonal matrix$Y^{(n)}$such that$Y^{(n)}:=L^{(n)}R^{(n)}- \frac{1}{\delta(n)}I$.

We derive from the matrix form of $dLV$$Y^{(n+1)}=R^{(n)}Y^{(n)}(R^{(n)})^{-1}$.

It is not hard to see $w_{k}^{(n)}>0$ providing $u_{k}^{(0)}>0$ andfor $k=1,2,$ $\cdots,2m-1$ . Thus
is nonsingular for any . This similarity transfomation implies that the eigenvalues ofare invariant under the evolution $n\Rightarrow n+1$ of the $dLV$ system. A symmetrization ofis introduced in [11] by using a diagonal matrix
as follows:$\delta^{(n)}>0$$R^{(n)}$$n$$Y^{(n)}$$Y^{(n)}$$G^{(n)}$$A^{(n)}$$:=$$(G^{(n)})^{-1}Y^{(n)}G^{(n)}$,$G^{(n)}$

$:=$

diag

$( \prod_{j=1}^{m-1}\sqrt{w_{2j-1}^{(n)}w_{2j}^{(n)}},$$\ldots,$$\sqrt{w_{2m-3}^{(n)}w_{2m-2}^{(n)}},1)$.

94Note thatis nonsingular for any andtakes the form of similarity transformation$G^{(n)}$$n$$|A^{(n)}|= \prod_{j=1}^{m}w_{2j-1}^{(n)}$$A^{(n+1)}=P^{(n)}A^{(n)}(P^{(n)})^{-1}$,. We see that the $dLV$ system$P^{(n)}:=(G^{(n+1)})^{-1}R^{(n)}G^{(n)}$, which implies that the eigenvalues ofof the positive definite matrix, for all , are$n\Rightarrown+1$
invariant under the evolutionof $dLV$.
Since the eigenvalues ofare identically equal to those of, these eigenvalues are independent from the choice of the variable step size
. Note also thatcan be decomposed
into$A^{(n)}$$A^{(n)}$$A^{(n)}$$n$$Y^{(0)}$$A^{(n)}$$\delta^{(n)}$$\sqrt{w_{3}^{(n)}}\sqrt{w_{2}^{(n)}}..$$A^{(n)}$$=$$(B^{(n)})^{T}B^{(n)}$,.

$B^{(n)}:=(\sqrt{w_{1}^{(n)}}O$.
$\sqrt{w_{2m-2}^{(n)}}\sqrt{w_{2m-1}^{(n)}}$Therefore the singular values ofare equal to the positive square roots of the eigenvalues of. Then the singular values of the upper bidiagonal matrix
are invariant under the time
evolution $n\Rightarrow n+1$ of the $dLV$ system.
Numerical stability of the $dLV$ algorithm is proved as follows. The positivity of the parameterand the variable
play a key role. In Ref. [121 the condition $0<\delta^{(n)}\leq M$ areassumed. As is shown in \S 3 of this paper the condition naturally follows from the positivity ofthe sequence of $kemel0Ps$ . By taking trace of the similarity transformation we see$B^{(n)}$$A^{(n)}$$B^{(n)}$$\delta^{(n)}$$u_{k}^{(n)}$$\sum_{k=1}^{2m-1}w_{k}^{(n)}=\sum_{k=1}^{m}\sigma_{k^{2}}$.

are bounded as well as positive. Consequently
Namelyare also positive and bounded.
$dLV$
Let $k=1$ in thesystem, then we have
. Thisimplies. Since, is monotonically increasing and$n=0,1,$bounded,converges to some positive constant as $narrow\infty$ . Simultaneously,
converges to some positive constant .
Let us assume thatconverges to some positive constant. Let$v_{k}^{(0)}>0$$0<\delta^{(n)}<M$
and. Then, by using, we see that $(v_{2k-1}^{(0)}/pk-1) \prod_{n=1}^{N}(1+$$Narrow\infty$
converges toas. Hence it follows that $0< \prod_{n=0}^{\infty}(1+\delta^{(n)}u_{2k}^{(n)})<M_{3}$for some constant. It is also obvious that, is monoton$N=1,2,$
ically increasing. Therefore it follows that. Simultaneously, we see$\lim_{narrow\infty}u_{2k-1}^{(n)}=v_{2k-1}^{(0)}pk/pk-1>0$, namely,that$w_{k}^{(n)}$$u_{k}^{(n)}$$\lim_{narrow\infty}u_{1}^{(n+1)}=u_{1}^{(0)}\prod_{n=0}^{\infty}(1+\delta^{(n)}u_{2}^{(n)})$$u_{1}^{(0)}\leq u_{1}^{(1)}\leq\cdots u_{1}^{(n)}\leq\cdots$$u_{1}^{(n)},$$\cdots$$u_{1}^{(n)}$$\prod_{n=0}^{\infty}(1+$$c_{1}$$\delta^{(n)}u_{2}^{(n)})$$p_{1}$$\prod_{n=1}^{\infty}(1+\delta^{(n)}u_{2k-2}^{(n)})$$p_{k-1}$$v_{k}^{(0)}=$$u_{k}^{(0)}(1+\delta^{(0)}u_{k+1}^{(0)})$$u_{2k-1}^{(N+1)}$$\delta^{(n)}u_{2k}^{(n)})$$M_{3}$$\prod_{n=1}^{N}(1+\delta^{(n)}u_{2k}^{(n)}),$$\cdots$$\prod_{n=1}^{\infty}(1+\delta^{(n)}u_{2k}^{(n)})=p_{k}$$\lim_{narrow\infty}u_{2k-1}^{(n)}=c_{k}$where is some positive constant.
Note here thatconverges to some constant $s_{k}>0$ if and only if$c_{k}$$\sum_{n=0}^{\infty}\delta^{(n)}u_{2k}^{(n)}$$\prod_{n=1}^{\infty}(1+$95$\delta^{(n)}u_{2k}^{(n)})=p_{k}$for$\delta^{(n)}u_{2k}^{(n)}>0,$

bounded sequence$\delta^{(n)}$, if$n=0,1,$$\cdots$. Moreover$\sum_{n=0}^{\infty}\delta^{(n)}u_{2k}^{(n)}=s_{k}$. Therefore it follows that$\lim_{narrow\infty}u_{2k}^{(n)}=0$

Note here that lim$\lim_{narrow\infty}\delta^{(n)}u_{2k}^{(n)}=0$

for any positive.

diag. This implies that is one of the eigenval, namely, the square of a singular value of
ues of. Singular values ofare equal tothose of $B=B^{(0)}$ . It is concluded that the $dLV$ algorithm converges to singular values of thebidiagonal matrix with nonzero diagonals and sub-diagonals in numerically stable way. The
Christoffel transformation of symmetric OPs gives rise to the positivity and boundedness of theand the variable
parameterof the $dLV$ algorithm. No subtraction appears in $dLV$. Thena higher accuracy follows.
$narrow\infty^{A^{(n)}}=$$(c_{1}, c_{2}, \cdots, c_{m})$$c_{k}$$A^{(n)}$$B^{(n)}$$B^{(n)}$$B$$\delta^{(n)}$5$u_{k}^{(n)}$Recent Developments on Shifts of mdLVs Algorithm

On speed the $dLV$ algorithm has a first order convergence providing that
. It is slow.
$dLV$
The shiftforbrings an ”internal shift” of thealgorithm. To
accelerate the convergence an ”extemal” shift of origin is introduced in [13] through a mapping$:=(B^{(n)})^{T}B^{(n)}-(\theta^{(n)})^{2}I$ . If asummation of shifts
is less
than the squareof the minimal singular value of , then$\delta^{(n)}>0$$1/\delta^{(n)}-1/\delta^{(n+1)}$$R^{(n)}L^{(n)}$$(B^{(n)})^{T}B^{(n)}arrow(\overline{B}^{(n)})^{T}\overline{B}^{(n)}$$(\theta^{(n)})^{2}$$B$$\sigma_{m}^{2}$$\lim_{narrow\infty}u_{2k-1}^{(n)}=\sigma_{k^{2}}-\sum_{n}(\theta^{(n)})^{2}$,$\lim_{narrow\infty}u_{2k}^{(n)}=0$

is shown. The positivity and boundedness of the vaniable
are not violated by the shift.
astableis[17],
shift given by using the Johnson boundSuchfor example. Then a new stablealgorithm with shift for bidiagonal singular value problem results which is named the mdLVs
algorithm. The mdLVs has two types of parameters. One is the intemal parameter. Theother is the extemal shift parameter. The mdLVs algorithm is more accurate than theDemmel-Kahan QR algorithm, the Divide&Conquer algorithm and the dqds algorithm which
are practically used bidiagonal singular value computing algorithms (cf. [4]) through the present
LAPACK codes [20]. The mdLVs algorithm is faster than the Demmel-Kahan QR, Divide&
Conquer the as well as the bisection algorithm. On these established algorithms see the book[4] and references therein.
The Johnson bound has been adopted in the original implementation [28, 29] of the mdLVs
algorithm [13]. The Johnson shift is$u_{k}^{(n)}$$\delta^{(n)}$$\theta^{(n)}$$\Theta^{(m)}$$;=$$\min_{k=1,\ldots m},(\sqrt{w_{2k-1}}-\frac{1}{2}(\sqrt{w_{2k-2}}+\sqrt{w_{2k}})\}$$<$$\sigma_{m}$,where $w_{0}=0,$ $w_{2k-1}+w_{2k-2}$ are the $(k, k)$ -elements and$(k+1, k)$ -elements of the tridiagonal matrix $A=B^{T}B$ .

$\sqrt{w_{2k}w_{2k-1}}$

are the $(k, k+1)$ and96
Recently the following lower bound is found by K. Kimura$\Theta_{p}^{(m)}$

$:=$

$($

trace

$(B^{T}B)^{-p})^{-\frac{1}{2p}}$

1$( \frac{1}{\sigma_{1^{2p}}}+\cdots+\frac{1}{\sigma_{m^{2p}}})^{\iota/2p}$$<$$\sigma_{m}$,$(p=1,2, \ldots)$which is called the p-th generalized Newton bound. The cost for the generalized Newton shift
is shown in [18] to be only $O(m)$ . Y. Yamamoto [191 proves that the generalized Newton shift
performs a weakly $(p+1)$ -th order convergence. The mdLVs code with the generalized Newton
shift where $p=2,3,4$ is faster and more accurate than the mdLVs code with the Johnson shift.
Here$I_{p}$$:=$

trace $(B^{T}B)^{-p}$

are conserved quantities of the $dLV$ system. This strongly suggests thatare also expressed
thisconjecture is proved affirmatively for by T.
only by using positive variables. Recently$O(m)$
formula for computing
is presented in [32]. The
Yamashita and a subtraction freegeneralized Newton shift will be useful for the dqds algorithm [23, 241 for singular values.
$\Theta_{p}^{(m)}$$\Theta_{p}^{(m)}$6 Double Cholesky Decomposition and dLV-type IPansrormation
Let us assume that all of the singular values of an $m\cross m$ upper bidiagonal matnix are
be the computed singular value. In thispositive, simple and are already computed. Let
section we introduce the dLV-type transformation for computing the right singular vector
corresponding to each
. The right singular vector ofand the left singular vector
is a solution vector $v_{j}=$ $(v_{j}(1), v_{j}(2), \cdots , v_{j}(m))^{T}$ of the system of linear equation$B$$\hat{\sigma}_{j}$$v_{j}$$\hat{\sigma}_{j}[22|$$u_{j}$$(B^{T}B-\hat{\sigma}_{i^{2}}I)v_{j}=0$$v_{j}$$B$.

usually contains some errors though the mdLVs algorithm has
Computed singular value
a higher relative accuracy. Conversely, if is a correct singular vector corresponding to thefor an approximant
correct singular value , thenof . Therefore let
us find more accurate singular vector by solving the linear equation$\hat{\sigma}_{j}$$v_{j}$$\sigma_{j}$$(B^{T}B-\hat{\sigma}_{j}^{2}I)v_{j}\neq 0$$\hat{\sigma}_{j}$$\sigma_{j}$$(B^{T}B-\hat{\sigma}_{j^{2}}I)v_{j}=c_{j}$for a suitable constant vector $c_{j}\neq 0$ . A derivation of the residual vector
will be describedbecomes singular. Thus we use ais close to , the coefficient matrix
later. As$c_{j}$$\hat{\sigma}_{j}$$\sigma_{j}$$B^{T}B-\hat{\sigma}_{j}^{2}I$direct method for solving the ill-conditioned linear equation $(B^{T}B-\hat{\sigma}_{j^{2}}I)v_{j}=c_{j}$ .
Any positive definite real symmetric matrix can be decomposed into the product of a lower(or upper) triangular matrix and its transposed. This is called the Cholesky decomposition. Oncethe coefficient matrix of linear equation is decomposed, the linear equation is solved through an

97
inversion of the triangular factor at a lower computational cost than the Gaussian elimination. In
this case we first consider the so-called “double Cholesky decomposition” of a positive definite
matrixas$B^{T}B-\hat{\sigma}_{j^{2}}I$$B^{T}B-\hat{\sigma}_{4^{2}}I=\{\begin{array}{l}(B^{+})^{T}B^{+},(B^{-})^{T}B^{-},\end{array}$.

$b_{3}^{+}b_{2}^{+}$$B^{+}:=(b_{1}^{+}0$,$B^{-}:=(b_{2}^{-}b_{1}^{-}$$b_{3}^{-}b_{2m-2}^{-}b_{2m-1}^{-}0$,$b_{Am-1}^{+}b_{2m-2}^{+}r$

and
where
are upper and lower bidiagonal matrices, respectively. If
isnonsingular but indefinite, namely, if is greater than the minimal singular valueof $B$ , thedouble Cholesky decomposition of complex type can be introduced similarly.
A problem arises. Cholesky decomposition of such an ill-conditioned matnx as
is liable to be numerically unstable. It is difficult to compute accurate triangular factors
and. Note that the Cholesky decomposition takes the form of the shift mapping of the mdLVs
algorithm. Whereas an intemal shift of the matrix representation of $dLV$, wherealgorithm is caused by the differenceandare parameters of $dLV$.
Therefore, we divide
into two$B^{+}$$B^{-}$$B^{T}B-\hat{\sigma}_{j^{2}}I$$\hat{\sigma}_{j}$$\sigma_{m}$$B^{T}B-\hat{\sigma}_{i^{2}}I$$B^{+}$$B^{-}$$B^{T}B-\hat{\sigma}_{j^{2}}I=:\overline{B}^{T}\overline{B}$$1/\delta^{(0)}-1/\delta^{(1)}$$\delta^{(0)}$$\delta^{(1)}$$\hat{\sigma}_{j^{2}}$$\hat{\sigma}_{j^{2}}=\frac{1}{\delta_{+}^{(0)}}-\frac{1}{\delta_{+}^{(1)}}$.

Consequently the target Cholesky decomposition can be divided into three$B^{T}B- \frac{1}{\delta_{+}^{(0)}}I=(\mathcal{W}^{(0)})^{T}\mathcal{W}^{(0)}$$(\mathcal{W}^{(0)})^{T}\mathcal{W}^{(0)}=(\mathcal{W}^{(1)})^{T}\mathcal{W}^{(1)}$,,$( \mathcal{W}^{(1)})^{T}\mathcal{W}^{(1)}+\frac{1}{\delta_{+}^{(1)}}I=(B^{+})^{T}B^{+}$

where$\mathcal{W}_{3}^{(l)}\mathcal{W}_{2}^{(l)}$$\mathcal{W}^{(p)}:=(\mathcal{W}_{1}^{(l)}0$

.

,

,

$\mathcal{W}_{k}^{(\ell)}:=\sqrt{u_{k}^{(p)}(1+\delta_{+}^{(\ell)}u_{k-1}^{(l)})}$$\mathcal{W}_{2m-1}^{(\ell)}\mathcal{W}_{2m-2}^{(l)}$The first equation reads$b_{2k-1^{2}}= \frac{1}{\delta_{+}^{(0)}}(1+\delta_{+}^{(0)}u_{2k-2}^{(0)})(1+\delta_{+}^{(0)}u_{2k-1}^{(0)})$$b_{2k}^{2}=\delta_{+}^{(0)}u_{2k-1}^{(0)}u_{2k}^{(0)}$,$u_{0}^{(0)}\equiv 0$

.

,

,

$(l=0,1)$ .

98From the assumption on the given bidiagonal matrix we see $b_{2k-1}\neq 0$ and $b_{2k}\neq 0$ . Thusand
. The parametershould not be equal to the inverse of any eigenvalueof $B^{T}B$ . Then, since eigenvalues of $B^{T}B-1/\delta_{+}^{(0)}I$ are simple and nonzero,
and 1
. Let
be a lower bound of the minimal singular value of $B$ . We
choose a positive number
so that. Then are computed sequentially. Since$B^{T}B-1/\delta_{+}^{(0)}I$
is not ill-conditioned, in general, the Cholesky factor
can be computed$B^{T}B-1/\delta_{+}^{(0)}I$
stably. It follows fromthatis positive definite. When $B$ has atiny singular values, we can choose a positive number
so that, whereare
singular values of $B$ computed by mdLVs. In this case the elementsof the Cholesky factor
becomes pure imaginary. But there is no additional difficulty.
The second equation is written as$B$$u_{2k-1}^{(0)}\neq 0$$u_{2k}^{(0)}\neq 0$$\delta_{+}^{(0)}$$1+\delta_{+}^{(0)}u_{2k-2}^{(0)}\neq 0$$+\delta_{+}^{(0)}u_{2k-1}^{(0)}\neq 0$$\overline{\sigma}_{m}$$\delta_{+}^{(0)}>1/\sigma_{m}^{2}$$\delta_{+}^{(0)}$$u_{k}^{(0)}$$\mathcal{W}^{(0)}$$1/\delta_{+}^{(0)}<\overline{\sigma}_{m}^{2}$$\delta_{+}^{(0)}\neq 1/\hat{\sigma}_{k^{2}}$$\delta_{+}^{(0)}$$\hat{\sigma}_{k}$$\mathcal{W}_{k}^{(\ell)}$$\mathcal{W}^{(0)}$$u_{k}^{(0)}(1+\delta_{+}^{(0)}u_{k-1}^{(0)})=u_{k}^{(1)}(1+\delta_{+}^{(1)}u_{k-1}^{(1)})$,$\delta_{+}^{(1)}:=\frac{\delta_{+}^{(0)}}{1-\delta_{+}^{(0)}\hat{\sigma}_{j^{2}}}$,$u_{0}^{(1)}\equiv 0$.

toThis is a transformation fromwith the parameter. Sinceand $(1+$in the matrix
are computed sequentially from given
and. If we,settemporarily, thenan identity transformation. Since the transformation$dLV$
is very similar to thesystem, the transformation is named the stationary $dLV$ (stdLV) in$u_{k}^{(1)}$$u_{k}^{(0)}$$\delta_{+}^{(0)}u_{k-1}^{(0)})\neq 0,$$u_{k}^{(0)}\neq 0$$\delta_{+}^{(0)}$$\mathcal{W}^{(1)}$$u_{k}^{(1)}$$u_{k}^{(0)}$$\delta_{+}^{(0)}$$u_{k}^{(1)}=u_{k}^{(0)}$$\delta_{+}^{(1)}=\delta_{+}^{(0)}$[16].

The third equation describes a process to generate bidiagonal matrix
throughwith parameter$u_{k}^{(1)}$$B^{+}$

from the variables$\delta_{+}^{(1)}$$\frac{1}{\delta(1)}(1+\delta_{+}^{(1)}u_{2k-2}^{(1)})(1+\delta_{+}^{(1)}u_{2k-1}^{(1)})=b_{2k-1}^{+2}$$\delta_{+}^{(1)}u_{2k-1}^{(1)}u_{2k}^{(1)}=b_{2k}^{+2}$

,

.

(1)

The left hand side of the flrst equation can be regarded as a shift of origin of the ill-posed. This is because $B^{T}B-1/\delta_{+}^{(0)}I=B^{T}B-\hat{\sigma}_{j^{2}}I-1/\delta_{+}^{(1)}I$ . By a suitable
matrixa possible numerical instability in the Cholesky decomposition
choice of$(B^{+})^{T}B^{+}$ can be avoided.
On the other hand the third equation recovers the factoroffromthe transformthe Cholesky decompositionofby the stdLV transformation. The
is
in
stdLV
following
relationship expressed the
diagram.
$B^{T}B-\hat{\sigma}_{j^{2}}I$$\delta_{+}^{(0)}$$B^{T}B-\hat{\sigma}_{j^{2}}I=$$B^{+}$$u_{k}^{(1)}$$\{b_{k}\}$$u_{k}^{(0)}$Cholesky decomposition–$\delta_{+}^{(0)}\downarrow$$\{u_{k}^{(0)}\}$$\{b_{k}^{+}\}$$\uparrow\delta_{+}^{(1)}$$arrow$
stdLV transformation

Fig. 1: stdLV diagram$\{u_{k}^{(1)}\}$

99
The other Cholesky decomposition $B^{T}B-\hat{\sigma}_{j^{2}}I=(B^{-})^{T}B^{-}$ is also divided into three$B^{T}B- \frac{1}{\delta_{-}^{(0)}}I=(\mathcal{W}^{(0)})^{T}\mathcal{W}^{(0)}$,$(\mathcal{W}^{(0)})^{T}\mathcal{W}^{(0)}=(\mathcal{V}^{(-1)})^{T}\mathcal{V}^{(-1)}$,$( \mathcal{V}^{(-1)})^{T}\mathcal{V}^{(-1)}+\frac{1}{\delta^{\underline{(}-1)}}I=(\mathcal{B}^{-})^{T}\mathcal{B}^{-}$

,

where
バ
$j^{2}= \frac{1}{\delta^{\underline{(}0)}}-\frac{1}{\delta^{\underline{(}-1)}}$,..

$\mathcal{V}^{(-1)}:=(\mathcal{V}_{1}^{(-1)}0\mathcal{V}_{3}^{(-1)}\mathcal{V}_{2}^{(-1)}$,$\mathcal{V}_{k}^{(-1)}=$sgn$(\mathcal{W}_{k}^{(0)})\sqrt{u_{k}^{(-1)}(1+\delta^{\underline{(}-1)}u_{k+1}^{(-1)})}$$\mathcal{V}_{2m1}^{(-}\mathcal{V}_{2m\frac{1)}{1)-}2}^{(-}$The second equation leads to$u_{k}^{(0)}(1+\delta^{\underline{(0})}u_{k-1}^{(0)})=u_{k}^{(-1)}(1+\delta_{-u_{k+1}^{(-1)}}^{(-1)})$

The mapping fromsuitable choice of,$\delta_{-}^{t-1)}:=\frac{\delta^{\underline{(}0)}}{1-\delta_{-}^{(0)}\hat{\sigma}_{j^{2}}}$,$u_{2m}^{(-1)}\equiv 0$.

tois called the reverse-time $dLV$ (rtdLV) transformation, By aa possible numerical instability in the Cholesky decomposition $B^{T}B-$
can be avoided. The relationship between variables is expressed in the rtdLV$u_{k}^{(0)}$$u_{k}^{(-1)}$$\delta_{-}^{(0)}$$\hat{\sigma}_{j^{2}}I=(B^{-})^{T}B^{-}$

diagram.

$\{b_{k}\}$Cholesky decomposition–$\delta_{-}^{(0)}\downarrow$$\{u_{k}^{(0)}\}$$\{b_{k}^{-}\}$$\uparrow\delta^{\underline{(}-1)}$$arrow$$rtdI_{\lrcorner}V$$\{u_{k}^{(-1)}\}$

rransrormation

Fig. 2: rtdLV diagram

We name the pair of stdLV and rtdLV as the dLV-type transformation. This performs the double Cholesky decomposition of a wide class of positive definite symmetric tridiagonal matrices$B^{T}B$ in anumerical stable way [16] by choosing suitable parameters. On the other hand,the qd-type transformation for computing eigenvectors of symmetric tridiagonals, proposed byParlett and Dhillon [5, 23] has no such parameter. Indeed there is anonsingular test matrix
having a tiny eigenvalue. The qd-type transformation make a serious error though the dLV-type
withgives an accurate Cholesky decomposition. Noted that is an eigenvector of the$\delta_{\pm}^{(0)}$$3\cross 3$$\delta_{\pm}^{(0)}=1$$v_{j}$.

100
symmetric positive definite tridiagonal matrix $B^{T}B$ . When a relative gap of eigenvalues is verysmall, both the qd-type and the dLV-type have challenges in orthogonality of resulting eigenvectors and singular vectors. For example, the qd-type fails to compute accurate eigenvectorsof the glued Wilkinson matrix [7].

7Twisted Factorization and Numerical Examples

In this section we explain a procedure for computing accurate right singular vectorsof$B$of double Cholesky decomposition. This part is essentially same asfrom the factors
the twisted factorization method by [5, 6, 231. The residual vector
in the linear system$(B^{T}B-\hat{\sigma}_{j^{2}}I)v_{j}=c_{j}$ is set as$v_{j}$$B^{\pm}$$c_{j}$$c_{j}=\gamma e$$\gamma j,k^{;=b_{2k-1}^{+2}+b_{2k-1}^{-2}-(b_{2k-2^{2}}+b_{2k-1^{2}}-\hat{\sigma}_{j^{2}})\neq 0}$,$e_{\rho}:=(0, \cdots, 0,1,0, \cdots, 0)^{T}$

,

,

where
are the residual parameters, the -th element of is 1 and is a number such thattakes the minimum for $k=\rho$ , namely, indicates the “most accurate” point. Then theis introduced as follows. Set
so-called “twisted matrix”$\gamma j,k$$\rho$$|\gamma_{j,k}|$$e_{\rho}$$\rho$$\rho$$N_{\rho}$$N(k)=\{\begin{array}{l}\frac{b_{2k}^{+}}{b_{2k-1}^{+}}, (k=1,2, \ldots, \rho-1),\frac{b_{2k}^{-}}{b_{2k+1}^{-}}, (k=\rho, \rho+1, \ldots, m-1),\end{array}$Ifis real, then so isreal. Define$b_{2k}^{\pm}$$D^{+}(k)=b_{2k-1}^{+2}$,$(k=1,2, \ldots, \rho-1)$ ,$D^{-}(k)=b_{2k-1}^{-2}$,$(k=\rho+1, \rho+2, \ldots, m)$$b_{2k\mp 1}^{\pm}$. If$b_{2k}^{\pm}$is pure imaginary, so is$b_{2k\mp 1}^{\pm}$.

. Therefore, $N(k)$ are always$1$.
$N(\dot{\rho}.-1)$$N_{\rho}.=[N(1)1$

1$N(m-1)1]$ ,$N(\rho)1^{\cdot}.\cdot$$D_{\rho}:=$diag $(D^{+}(1), \cdots, D^{+}(\rho-1), \gamma_{j_{2}\rho}, D^{-}(\rho+1), \cdots, D^{-}(m))$ .

The coefficient matrix of the linear system$(B^{T}B-\hat{\sigma}_{j^{2}}I)v_{j}=\gamma_{j,\rho}e_{\rho}$$B^{T}B-\hat{\sigma}_{i^{2}}I=N_{\rho}D_{\rho}(N_{\rho})^{T}$.

takes the form

101
Table 1: Accuracy and orthogonality of DBDSQR and I-SVD$||\hat{V}-V\Vert_{sum}$$\Vert B-\hat{U}^{\underline{\hat{\nabla}}}\hat{V}^{T}||_{sum}$$(\cross 10^{-9})$$\Vert\hat{V}^{T}\hat{V}-I||_{sum}$

DBDSQR

I-SVD

DBDSQR

I-SVD

DBDSQR

I-SVD

0.0690

398

102

4.14

0.0831

0.324

Table 2: Execution time of DBDSQR and I-SVD (sec.)$m=1000$$m=2000$$m=6000$

DBDSQR

I-SVD

DBDSQR

I-SVD

DBDSQR

I-SVD

4492

1.13

432.12491

4257360

4373

This is sometimes called the twisted factorization. The substitution for determining $N(k)$ and$D^{\pm}(k)$ can be done from the twisted point $k=\rho$ as $k=\rho,$ $\rho\pm 1,$$\rho\pm 2$ ,. . . . Whilevarious errors may be accumulated if we compute matrix factorizations by the usual one-side
substitution. The twisted factorization for each can be computed by $O(m)$ times of divisions.
Since $D_{\rho}e_{\rho}=\gamma e,$satisfies the linear system $(B^{T}B-$$B^{T}B$
and is an eigenvector ofproviding that$\hat{\sigma}_{j}$$N_{\rho}e_{\rho}=e_{\rho},$ $D_{\rho}N_{\rho}e_{\rho}=N_{\rho}D_{\rho}e_{\rho},$$v_{j}$$\hat{\sigma}_{j^{2}}I)v_{j}=\gamma j,\rho e_{\beta}$$N_{\rho}^{T}v_{j}=e_{\rho}$

The solution vector $v_{j}=(v_{j}(k))$ of$N_{\rho}^{T}v_{j}=e_{\rho}$.

is given by$v_{j}(k)=\{\begin{array}{ll}1, (k=\rho),-N(k)v_{j}(k+1), (k=\rho-1, \rho-2, \ldots, 1),-N(k-1)v_{j}(k-1), (k=\rho+1,\rho+2, \ldots, m),\end{array}$which gives rise to the right singular vector . Since each singular vector can be computed within $O(m)$ flops, the computation of k-singular vectors costs $O(km)$ flops for $k=$
$m$ .
$1,2,$
of $B$ are given through $U=BV\Sigma^{-1}$ , with $U=$
The left singular vectors$(u_{1}, \ldots, u_{m}),$ $V=(v_{1}, \ldots, v_{m}),$, or by solving the system $(BB^{T}-$diagdirectly.
Finally we quote some numerical examples from [301 on an implementation of the I-SVD algorithm. In Table 1 the accuracy of singular value decomposition, the accuracy of right singular
vectors and their orthogonality are considered. Here
indicates the sum of absolutevalues of every elements of the matrix $\hat{V}-V$ , whereinare computed$v_{j}$$\ldots,$$u_{j}$$\Sigma=$$(\hat{\sigma}_{1}, \ldots,\hat{\sigma}_{m})$$\hat{\sigma}_{j^{2}}I)u_{j}=0$$\Vert\hat{V}-V\Vert_{sum}$$\hat{v}_{j}$$\hat{V}=(\hat{v}_{1}, \ldots,\hat{v}_{m})$102
right singular vectors. The numerical numbers in Table 1 are averages of 1001000 $\cross 1000$
bidiagonal test matrices whose singular values are randomly distributed. Here DBDSQR is thestandard LAPACK code for bidiagonal SVD where the Demmel-Kahan QR algorithm is implemented. In the I-SVD code the left singular vectors are computed by
and are-onhogonalization of all of the singular vectors is performed by inverse iterations once. It
takes extra $O(m^{2})$ flops. The parameters
are fixed to 1. Table 1 shows that DBDSQR is$1\sim2$
digits on the accuracy of SVD and the orthogonality ofbetter than the I-SVD code bysingular vectors. However on the accuracy of singular vectors I-SVD is better than DBDSQR.
Table 2 is a comparison of execution time between DBDSQR and the I-SVD code. It is
obvious that the I-SVD algorithm needs only $O(m^{2})$ flops and is rather faster than DBDSQR
of $O(m^{3})$ flops. I-SVD also has a better scalability. The I-SVD code with the Householder
preconditioning to bidiagonal manices and an inverse transformation is still sufficiently fasterthan DBDSQR with Householder [30].
$\hat{U}=B\hat{V}\Sigma^{-1}^{\wedge}$$\delta_{\pm}^{(0)}$8. Concluding Remarks
This report surveys recent developments of the mdLVs algorithm for singular values and the
I-SVD algorithm for bidiagonal SVD. The mdLVs is a shifted version of the $dLV$ algorithm.
We first show how the $dLV$ algorithm has a higher accuracy. The Christoffel transformationof symmetric OPs gives rise to the positivity and boundedness of the parameterand theof the $dLV$ algorithm. No subtraction appears in $dLV$. Positivity is also essential
variablein the formulation of the mdLVs algorithm. Namely, if a shift is less than the minimal singularvalue, then the positivity of mdLVs follows.
The Johnson bound [171 has been adopted in the mdLVs code [291. Recently a new lower
bound is found which is called the p-th generalized Newton bound. The generalized Newton
shift costs only $O(m)$ flops where $m$ is the size of given bidiagonal matrix [18]. Y. Yamamoto[19] proves that the generalized Newton shift perfoims a weakly $(p+1)$ -th order convergence.
The mdLVs code with the generalized Newton shift where $p=2,3,4$ is faster and more accurate than the mdLVs code with the Johnson shift. Though lower and upper bounds of matrix
eigenvalues have been studied fully [311, exploring for new bound is still an important problem
in numerical linear algebra.
The I-SVD algorithm is a combination of the mdLVs algorithm and the dLV-type transformation for singular vectors. Because of the separation of computation of singular values fromthat of singular vectors the I-SVD algorithm runs in $O(m^{2})$ flops. Whereas the Demmel-Kahan
QR algorithm requires $O(m^{3})$ flops. Thus the I-SVD code is rather faster than DBDSQR codeof LAPACK. The I-SVD code has a good orthogonality of singular vectors for the case of random matrices. To improve the orthogonality for clustered matrices the I-SVD algorithm should
be investigated more.
$\delta^{(n)}$$u_{k}^{(n)}$

103
参考文献
[1]

N. I. Akheiezer, The Classical Moment Problem and Some Related Questions in Analysis, Edinburgh: Olver&Boyd, 1965.

[2] T. S. Chihara, An Introduction to Orthogonal Polynomials, New York: Gordon&Breach, 1978.
[3][4]$pearMoody$

T. Chu, “Linear algebra algorithms as dynamical systems,‘’

Acta Numerica, 2008, to ap-

J. Demmel, Applied Numerical Linear Algebm, Philadelphia: SIAM, 1997.

[5] I. S. Dhillon and B. N. Parlett, Orthogonal eigenvectors and relative gaps, SIAM J. Matrix Anal.
Appl., 25(2004), 858-899.
[6] I. S. Dhillon and B. N. Parlett, Multiple representations to compute orthogonal eigenvectors of
symmetric tridiagonal matrices, Lin. Alg. Appl., 387(2004), 1-28.
[7]

I. S. Dhillon, B. N. Parlett and C. Vomel, Glued matrices and the MRRR algorithm, SIAM J. Sci.
Comput., 27(2005), 496-510.

[8]K. V. Femando, B. N. Parlett, Accurate singular values and differential qd algorithms, Numer.
Math., 67 (1994), 191-229.

[9]G. H. Golub and C. F. Van Loan, Matrix Computations 3rd $edn$ , Baltimore: The Johns Hopkins
Univ. Press, 1996.

[10]

R. Hirota, ”Conserved quantities of random-time Toda equation”, J. Phys. Soc. Japan Vol. 66, pp.
283-284, 1997.

[11]

M. Iwasaki and Y. Nakamura, “On a convergence of solution of the discrete Lotka-Volteira system;’Inverse Problems, vol. 18, pp. 1569-1578, 2002.

[12] M. Iwasaki and Y. Nakamura, “An application of the discrete Lotka-Volterra system with variablestep-size to singular value computation,” Inverse Problems, Vol. 20, pp. 553-563, 2004.

[131 M. Iwasaki and Y. Nakamura, “Accurate computation of singular values in terms of shifted inte-

grable schemes, Japan J. Indust. Appl. Math. Vol. 23, pp. 239-259, 2006.
”[14] M. Iwasaki and Y. Nakamura, ”Positivity and stability of the $dLV$ algorithm for computing matrix

singular values,” in:Proceedings of the 2nd Intemational Conference on Informatics Research forDevelopment of Knowledge Society Infrastructure ICKS 2007, $Eds$ . M. Fukushima, K. Tajima andK. Tanaka, IEEE Computer Society Press, pp. 103-110, 2007.

[15]

M. Iwasaki and Y. Nakamura, ”Center manifold approach to discrete integrable systems related toeigenvalues and singular values,’‘ Hokkaido Math. J., Vol. 36, pp. 759-775, 2007.

[16] M. Iwasaki, S. Sakano and Y. Nakamura, ‘Accurate twisted factorizat\’ion of real symmetric tridiagonal matrices and its application to singular value decomposition,” Trans. Japan SIAM, 15(2005),

461-481. (岩崎雅史. 阪野真也, 中村佳正, 実対称 3 重対角行列の高精度ツィスト分解とその特異値分解への応用,[17]

日本応用数理学会論文誌).

C. R. Johnson, ‘A Gersgorin-type lower bound for the smallest singular value;‘ $Lin$ . Alg. Appl.,
Vol. 112, pp. 1-7, 1989.

[18] K. Kimura, T. Yamashita and Y. Nakamura, “On $O(N)$ formula for the diagonal elements of inverse

powers of symmetric positive definite tridiagonal matrices,’‘ in preparation.

[19] K. Kimura, Y Yamamoto, M. Takata and Y Nakamura, “Generalized Newton shifts: An arbitrarily

high order shufting scheme for singular value computation,” in preparation.

[20]

LAPACK. http $;//www$.netlib. $orgnapack/$

104
new approach to numerical algorithms in terms of integrable systems,”
in:Proceedings of Intemational Conference on Informatics Research for Development of KnowlIbaraki, T Inui and K. Tanaka, IEEE Computeredge Society Infrastructure ICKS 2004, $Eds$.
Society Press, pp. 194-205, 2004.

[21] Y. Nakamura, ‘A$T$[221 Y. Nakamura, Functionality of Jntegrable Systems, Tokyo: Kyoritsu Shuppan, 2006, (中村佳正著「可積分系の機能数理」 共立出版).

[23] B. N. Parlett and I. S. Dhillon, Femando’s solution to Wilkinson’s problem:
double factorizatio , Lin. Alg. Appl., 267(1997), 247-279.

An application of

$n$

[24]

B. N. Parlett and $0$ . A. Marques, An implementation of the dqds algorithm (positive case), Lin.
Alg. Appl., 309(2000), 217-259.

[25] H. Rutishauser. “Ein Quotienten-Differenzen-Algorithmus,”233-251, 1954.

$7_{\wedge}$

angew Math. Phys., vol. 5. pp.

[26] V. $Sp\ddot{m}donov$ and A. Zhedanov, ”Discrete-time Volterra chain a $nd$ classical orthogonal polynomials,” J. Phys. $A$ , vol. 30, pp. 8727-8737, 1997.

[27]

G. Szeg\"o, Orthogonal Polynomials, Providence: Amer. Math. Soc., 1939.

[28] M. Takata, M. Iwasaki, K. Kimura and Y. Nakamura, “Implementation and its evaluation of rou-tine for computing singular values with high accuracy,” IPSJ Trans. Adv. Comput. Systems, 46,
No SIG12(2005), 299-311, (高田雅美, 岩崎雅史, 木村欣司, 中村佳正, 高精度特異値計算ルーチンの開発とその性能評価, 情報処理学会論文誌).

[29] M. Takata, M. Iwasaki, K. Kimura and Y. Nakamura. “An evaluation of singular value computationby the discrete Lotka-Volterra system,” in:Proceedings of The 2005 Intemational Conference on
Parallel and Distributed Processing Techniques and Applications (PDPTA2005), Vol. II, pp. 410416, 2005.

[30] M. Takata. K. Kimura. M. IwasA and Y. Nakamura, “Implementation of library for high speed
singular value decomposition2” IPSJ Trans. Adv. Comput. Systems, 47, No.SlG 7(ACS 14), (2006),
91-104, ( ffl% ,,のためのライブラ )Iま $*$ g{F9$\grave$$\ovalbox{\tt\small REJECT}$$\not\equiv$$*\}\backslash \dagger I\overline{\text{ロ}}\urcorner,$$\epsilon ur\ovalbox{\tt\small REJECT}\Re R,$$\mathfrak{c}P\}\backslash \dagger\dagger\not\leqq iE$$\dagger$$\ovalbox{\tt\small REJECT}$$\mathfrak{X},$[31]$’\ovalbox{\tt\small REJECT}\Re\infty g_{\mp\hat{=}p|wX_{\overline{p}}^{-}T)}^{r^{\backslash }3a}$$\beta\not\in$$\ovalbox{\tt\small REJECT}$.

R. Varga, GerSigorin and His Circles, Heidelberg: Springer-Verlag, 2004.

[32] T. Yamashita, K. Kimura and Y. Nakamura, “On subtraction free formula for the diagonal elements

of inverse powers of positive definite tridiagonal matrices,” in preparation.