Table Of Content. .
ARTICLES doi:
A new shift strategy for the implicitly restarted
7
1 generalized second-order Arnoldi method
0
2
n GONG FangHui & SUN YuQuan∗
a
J
1 School of Mathematics and SystemsScience, BeiHang University,Beijing 100191, China;
1
Email: [email protected], [email protected]
] ReceivedJune24,2015; accepted June24,2015
A
N
. Abstract In this paper, a new shift strategy for the implicitly restarted generalized second-order Arnoldi
h
(GSOAR) method is proposed. In implicitly restarted processes, we can get a k-step GSOAR decomposition
t
a from a m-step GSOAR decomposition by performing p = m−k implicit shifted QR iterations. The problem
m of the implicitly restarted GSOAR is the mismatch between the number of shifts and the dimension of the
[ subspace. There are 2p shifts for p QR iterations. We use the shifts to filter out the unwanted information
in the current subspace; when more shifts are used, one obtains a better updated subspace. But, if we use
1 morethan pshifts,thestructure of theGSOAR decomposition willbedestroyed. Wepropose anovel method
v whichcanuseall2pcandidatesandpreservethespecialstructure. Thenewmethodvastlyenhancestheoverall
2
efficiencyofthealgorithm. Numericalexperimentsillustratetheefficiencyofeveryrestartprocess.
4
0
Keywords QEP,GSOARmethod, implicitlyrestart,shifts,Ritzvector
3
0 MSC(2010) 65F15, 65Y20
.
1
0 Citation: GongFH,SunYQ. Anewshiftstrategyfor theimplicitlyrestartedgeneralizedsecond-orderArnoldi
7
method (in Chinese). Sci Sin Math, 2017, 47: 116, doi: 10.1360/012016-22. You can see the Chinese
1 versiononline: http://engine.scichina.com/doi/10.1360/012016-22
:
v
i
X 1 Introduction
r
a
The quadratic eigenvalue problem (QEP)
Q(λ)x=(λ2M +λC+K)x=0 (1.1)
hasvariousapplicationsintheengineeringfields,likedynamicanalysisofacousticsystems,fluidmechan-
ics and microelectronicmechanicalsystems [3,12,15]. There are two majorclasses of numericalmethods
to solve large QEPs. One is to linearize the QEP into an equivalent generalized eigenvalue problem
(GEP) such as
−C −K λx M 0 λx
=λ . (1.2)
" I 0 #" x # " 0 I #" x #
Currently,theGEPtechniqueshavebeenquitemature[1,7]. Butthelinearizationmethodshavesuffered
some disadvantages: losing the spectral property of original QEP and doubling the dimension of the
originalone. AnotherclassofmethodsworkontheQEPdirectly. Numericallystableprojectionmethods
can keep the spectral properties as well as the structure of original QEP.
∗Correspondingauthor
(cid:13)c
2
GONGFangHuietal. SciChinaMathforReview
The Krylov subspace method plays a significant role in numerical techniques for solving large-scale
GEPs. Bai and Su developed the standard Krylov subspace to second-order and proposed the second-
order Arnoldi (SOAR) method for large-scale sparse QEPs [2]. They propose a SOAR procedure that
computes an orthonormal basis of a second-order Krylov subspace. The SOAR method then projects
(1.1) ontothis subspace andcomputes the Ritz pairs to approximatethe desiredeigenpairsof (1.1). The
generalconvergencepropertyofRayleigh-Ritzhasbeendemonstratedin[4]andtherefinedRayleigh-Ritz
orthonormal projection method also has been proposed.
On account of the limitation of high computational cost and undue storage requirement, explicit or
implicit restartingis usuallynecessarywhen we do large-scalematricescalculating. Implicit restartingis
a powerful and extensively used technique in eigenvalue problem calculations [13]. Although the SOAR
methodplaysanimportantroleinsolvingQEPs,wecannotdirectlyapplytheimplicitrestartingscheme
to the SOAR procedure due to the special structure of the initial vector [10,11]. Based on Bai and Su’s
work, many researchers work specifically to implicitly restart the SOAR procedure. Otto proposed a
generalized SOAR (GSOAR) method which substitutes general vector for the original special-structure
starting vector. The generalized SOAR procedure can directly apply implicit restarting technique under
thehypothesisofnodeflation[11]. AnexplicitrestartinggeneralKrylovsubspacemethodwasputforward
after making a similar modification by Zhou [16]. Huang brought forward the implicit restarting semi
orthogonal generalized Arnoldi (SGA) method [5]. At the same time, other researchers are constantly
improving the SOAR method in all aspects [9].
Among these modified restarting SOAR methods, they all encountered the problem that implicit
restarting failed to work when shifts did not match up with the dimensions of subspace. From recent
researches, it can be seen that the quantity and quality of shifts are of great significance [9,11]. Jia
and Sun explored more properties and features of the GSOAR and gave an efficient method to compute
refined Ritz vectors. They proposeda refined GSOAR (RGSOAR) method. Then they advancedcertain
exact shifts and refined shifts for respective use within the implicitly restarted GSOAR and RGSOAR
algorithms. TherefinedshiftsarebasedontherefinedRitzvectors. Theypresentedanefficientalgorithm
to compute the exact and refined shift candidates reliably. Unlike the implicitly restarted algorithms for
the linear eigenvalue problem, both exact and refined shift candidates are more than the shifts allowed.
So the authors showed how to reasonably select the desired shifts among them. In addition, they pro-
posed an effective approach to cure deflation in implicit restarts, so that implicit restarting is useable
unconditionally.
Based on Jia and Sun’s work, we give a new shift strategy for the implicitly restarted GSOAR and
RGSOARmethods whichcan use allshifts. Inthe implicit restarting,we use the unwantedapproximate
eigenvalues as shifts to filter out the corresponding eigenvector information from the current subspace.
A general result is that with more shifts approximating the unwanted eigenvalues, we can generate a
better updated subspace. For the QEPs, one eigenvector can correspond to two different eigenvalues.
If the eigenvalue pairs corresponding to the same eigenvector are present in exact or refined shifts, this
means some essential shifts will be abandoned. This can be avoided by using all of the shift candidates.
Obviously this will destroy the original structure. We use a subtle transformation making the damaged
structure back to the Hessenberg form. Due to all shifts, the unwanted information is excluded more
thoroughly and we can get a better subspace in each restarting. Because the implicit restarting in [9] is
useable unconditionally, our algorithm is useable inherently.
Therestofthispaperisarrangedasfollows. Insection2,weintroducetheSOARmethodandGSOAR
method, implicit restarting, selection of exact and refined shifts. In section 3, we present the new strat-
egy: using all shifts for projection method to enhance the efficiency of the method. In section 4, several
numerical experiments are presented to illustrate the efficiency of the new implicitly restarted GSOAR
method.
3
GONGFangHuietal. Sci ChinaMathforReview
2 The implicitly restarted GSOAR and RGSOAR methods
The SOAR method was proposedto solvethe large-scalesparseQEPs basedon the second-orderKrylov
subspace by Bai and Su [2].
Definition 2.1. Let matrices A, B ∈ Cn×n, non-zero vector u ∈ Cn, and define the sequence
r ,r ,r ,··· ,r based on A, B and u, where
0 1 2 n−1
r = u
0
r = Au
1
r = Ar +Br for j >2.
j j−1 j−2
The space
G (A,B;u)=span{r ,r ,r ,··· ,r } (2.1)
n 0 1 2 n−1
is called a nth second-order Krylov subspace.
If we apply standard Krylov subspace technique to (1.2), then the corresponding Krylov subspace is:
κ (H,ν)=span{ν,Hν,H2ν,··· ,Hn−1ν}, (2.2)
n
u
with ν = . After some simple calculations, we note that there are some associations between the
"0#
standard Krylov vectors Hjν of length 2n in (2.2)
r
Hjν = j ,j >1, (2.3)
"rj−1#
and the vectors of second-order Krylov subspace
r =Ar +Br . (2.4)
j j−1 j−2
The subspace (2.2) based on matrix H and starting vector v can be generalized by vector sequence
r ,r ,r ,··· ,r . In the meantime, (2.3) indicates that the second-order Krylov subspace G (A,B;u)
0 1 2 n−1 j
can be used as the projection subspace of QEPs while we only need to generate an orthonormal basis
q n :
jj=1
G (A,B;u)=span{q ,q ,··· ,q },j >1. (2.5)
j 1 2 j
Bai and Su presented a procedure for generating an orthonormal basis of the second-order Krylov
subspace G (A,B;u).
j
Algorithm 1. SOAR procedure with deflation remedy
1: q =u /ku k
1 1 1 2
2: p =0
1
3: for j =1:m
4: r=Aq +Bp
j j
5: s=q
j
6: for i=1:j
7: t =qTr
ij i
8: r =r−q t
i ij
4
GONGFangHuietal. SciChinaMathforReview
9: s=s−p t
i ij
10: end for
11: t =kr k
j+1,j 2
12: if t =0,
j+1,j
13: if s∈span{p |i:q =0,i6i6j}
i i
14: break
15: else deflation
16: t =1
j+1,j
17: q =0
j+1
18: p =s
j+1
19: end if
20: else
21: p =s/t
j+1 j+1,j
22: q =r/t
j+1 j+1,j
23: end if
24: end for
At stepj, if r , i=0,1,...,j are linearlydependent but [rT,rT ]T, i=0,...,j with r =0 are not,
i i i−1 −1
we call this situation deflation; if both {r } and {[rT,rT ]T} are linearly dependent at step j, we call
i i i−1
this situation breakdown. If t = 0 at a certain step j (j < k), that means deflation occurs. In this
j+1,j
case the second-order Krylov subspace G (A,B;u) does not contain any exact eigenvector of (1.1) [9].
j
Therefore, the measure to remedy deflation is to reset t a non-zero value, here simply set it to be
j+1,j
one, let q = 0, and continue the algorithm. The non-zero vectors of sequence {q } still maintain
j+1 j
orthonormal and span the subspace G (A,B;u).
j
For the quadratic eigenvalue problem (1.1) of n dimension , it corresponds to a second-order Krylov
subspaceG (A,B;u), withm6n,andA=−M−1C, B =−M−1K. LetQ =[q ,q ,··· ,q ]∈Cn×m,
m m 1 2 m
T
P = [p ,p ,··· ,p ] ∈ Cn×m and unreduced upper Hessenberg matrix Tˆ = [t ] = m ∈
m 1 2 m m ij
"tm+1,meTm#
C(m+1)×m. If Algorithm 1 does not stop before m, then we have
span{Q }=G (A,B;u), (2.6)
m m
and define
Q Q
H m = m+1 Tˆ (2.7)
m
"Pm# "Pm+1#
A B
as m-step SOAR decomposition, where Q =[Q ,q ], P =[P ,p ], H = .
m+1 m m+1 m+1 m m+1
"I 0#
Rayleigh-Ritz method can be adapted in (1.1) to seek the approximate eigenpairs (θ,y) satisfying the
Galerkin restriction:
(θ2M +θC+K)y⊥G (A,B;u), (2.8)
m
5
GONGFangHuietal. Sci ChinaMathforReview
where y ∈ G (A,B;u), θ ∈ C. The columns of n × n matrix Q generated by Algorithm 1 are an
m m
orthonormal basis of G (A,B;u). Let θ, g satisfy the small-scale QEP:
m
(θ2M +θC +K )g =0, (2.9)
m m m
whereM =QTMQ ,C =QTCQ ,K =QTKQ ,thenwecanobtaintheapproximateeigenpairs
m m m m m m m m m
θ and y = Q g of (1.1) . The small-scale QEP (2.9) is produced by projecting the original QEP (1.1)
m
onto G (A,B;u) explicitly. This method is named second-order Arnoldi (SOAR) method since it is an
m
Arnoldi-like procedure. From the theory of Arnoldi method, we can find that the SOARmethod has the
advantages of fast convergence rate and simultaneous convergence of a group of eigenvalues. Moreover,
it can still keep the special structure of original QEP after projection.
Adisadvantageofthis method is thatimplicit restartingscheme cannotbe directly adapted. The new
starting vector p will no longer be zero after truncating. In order to satisfy the implicitly restarted
1
conditions, researchers suggested the GSOAR method which replaces p = 0 by p = u / k u k
1 1 2 2 2
and meets the requirement that starting vector p = u / k u k is non-zero, where u is a non-zero
1 2 2 2 2
vector. Combining the deflation remedy SOAR procedure with GSOAR procedure, we get the following
algorithm:
Algorithm 2. GSOAR procedure with deflation remedy
1: q =u /ku k , p =u /ku k .
1 1 1 2 1 2 2 2
2: for j =1,2,...,k do
3: r =Aq +Bp
j j
4: s=q
j
5: for i=1,2,...,j do
6: t =q∗r
ij i
7: r =r−t q
ij i
8: s=s−t p
ij i
9: end for
10: t =krk
j+1j
11: if t =0
j+1j
12: if s∈span{p |i:q =0,16i6j}
i i
13: break
14: else deflation
15: reset t =1
j+1j
16: q =0
j+1
17: p =s
j+1
18: end if
19: else
20: q =r/t
j+1 j+1j
21: p =s/t
j+1 j+1j
6
GONGFangHuietal. SciChinaMathforReview
22: end if
23: end for
IfAlgorithm2 does notstopbefore m, itgivesthe followinggeneralizedsecond-orderKrylovsequence
and subspace; see [9].
Let A and B be n×n matrices and for vectors u ,u ∈Cn, and define
1 2
r = u ,
0 1
r = Ar +Bu ,
1 0 2
r = Ar +Br for j >2.
j j−1 j−2
Then r ,r ,r ,...,r is called a generalized second-order Krylov sequence based on A,B and u ,u ,
0 1 2 m−1 1 2
and G (A,B;u ,u )=span{r ,r ,r ,...,r } the m-th generalized second-order Krylov subspace.
m 1 2 0 1 2 m−1
T
Let Q = [q ,q ,··· ,q ] ∈ Cn×m, P = [p ,p ,··· ,p ] ∈ Cn×m and Tˆ = [t ] = m ∈
m 1 2 m m 1 2 m m ij
"tm+1,meTm#
C(m+1)×m. Then Q is an orthonormal basis of the m-th generalized second-order Krylov subspace
m
G (A,B;u ,u ) and we get the m-step GSOAR decomposition:
m 1 2
Q Q
H m = m+1 Tˆ , (2.10)
m
"Pm# "Pm+1#
A B
where Q =[Q ,q ], P =[P ,p ], H = .
m+1 m m+1 m+1 m m+1
"I 0#
WecanprojecttheoriginalQEP(1.1)ontothegeneralizedsecond-orderKrylovsubspaceG (A,B;u ,u )
m 1 2
explicitly. We still get the small-scale QEP:
(θ2M +θC +K )g =0, (2.11)
m m m
where M = QTMQ , C = QTCQ , K = QTKQ . We can get the Ritz pairs θ and y = Q g
m m m m m m m m m m
of (1.1) by computing the eigenpairs of (2.11). Then we get the GSOAR method. Suppose that we
have computed the Ritz values θ by the GSOAR method. For each θ, we seek a unit length vector
u˜∈G (A,B;u ,u ) satisfying the optimal requirement
m 1 2
k(θ2M +θC+K)u˜k=arg min k(θ2M +θC+K)uk (2.12)
u∈Gm(A,B;u1,u2)
kuk=1
and use it as an approximate eigenvector, called the refined Ritz vector. The pairs (θ,u˜) are also called
therefinedRayleigh–RitzapproximationsandthemethodiscalledrefinedgeneralizedSOAR(RGSOAR)
method. The refinedvector can be computed by seeking a unit length vectorz˜∈Cm suchthat u˜=Q z˜
m
with
z˜=arg min k(θ2M +θC+K)Q zk, (2.13)
m
z∈Cm
kzk=1
the right singular vector of the matrix (θ2MQ +θCQ +KQ ) associated with its smallest singular
m m m
value σ (θ2MQ +θCQ +KQ ). It was shown in [9] that the singular vectors can be computed
min m m m
with moderate cost.
2.1 Implicit restarting
For the SOAR method and GSOAR method, in order to guarantee the convergence of the approximate
eigenpairs, m should be as large as possible theoretically. Nevertheless, the storage and computing
7
GONGFangHuietal. Sci ChinaMathforReview
requirements become massive as the dimension m of the subspace increases. So restarting is necessary
for a practical method. If no deflation occurs, it is direct to adapt the implicit restarting scheme [13] to
the GSOAR procedure. Given p shifts µ ,µ ,...,µ , we can apply p implicit shifted QR algorithm to
1 2 p
matrix T in (2.10)
m
(T −µ I)···(T −µ I)=V R, (2.14)
m 1 m p m
where V is a m×m orthogonal matrix and R is a m×m upper triangular matrix. By relation (2.10),
m
we have
Q Q q
H m V = m V (VTT V )+t m+1 eTV . (2.15)
m m m m m m+1,m m m
"Pm# "Pm# "pm+1#
Q˜ Q
Withsimplecomputation,itisnotedthatVTT V isaHessenbergmatrix. Ifwetake m = m V
m m m "P˜m# "Pm# m
and take VTT V =T˜ , then we have
m m m m
Q˜ Q˜ q
H m = m T˜ +t m+1 bT, (2.16)
"P˜m# "P˜m# m m+1,m"pm+1# m
where bTm =(0 ··· 0 bm−p ···bm−1 bm) which contains at least m−p−1 zeros.
Given a m-step GSOAR decomposition(2.16), it can be truncatedas a k-step GSOAR decomposition
while k 6m−p
Q˜ Q˜ q˜
H k = k T˜ +t k+1 eT. (2.17)
"P˜k# "P˜k# k k+1,k"p˜k+1# k
Meanwhile,wegetak-threstartedsubspaceG (A,B;u ,u )spannedbyanorthonormalvectorssequence
k 1 2
{q˜,q˜,··· ,q˜}.
1 2 k
BecausetheimplicitlyrestartedGSOARmethodhasastrictrequirementforthedimensionofrestarted
subspace,wecantruncate(2.16)intoanewk-threstartedsubspaceifandifonlythesubspacedimension
k 6 m−p. Specifically, for (2.17), when k does not exceed m−p, the last term in the right looks like
eT after truncating and (2.17) will have the same formas (2.10). The k-stepGSOAR decompositioncan
k
be expanded to m-step GSOAR decomposition again via Algorithm 2. Then looking for new Ritz pairs,
computing their deviation and the process repeats.
Equation (2.17) is equivalent to a k-step GSOAR procedure mathematically. If deflations occur at
steps m ,m ,...,m 6 k, then the corresponding j columns q˜ of Q˜ are zeros. Jia and Sun gave a
1 2 j mj k
method to cure this problem, so the implicit restarting is unconditional??.
2.2 Exact shifts and refined shifts
We considerthe computationofshifts for GSOARandRGSOARmethods in this subsection. Two kinds
of shifts were given in [9]: the exact shifts and refined shifts. We first show how to compute the refined
shifts as example.
The RGSOAR method computes the refined Ritz vectors u˜ , which can be much more accurate than
i
the Ritz vectors y .There were certain refined shifts for the refined Arnoldi method and the refined
i
harmonic Arnoldi method for the linear eigenvalue problem [7,8]. In the same spirit, the refined shifts
were proposed for the RGSOAR method.
Let span{uˆ ,...,uˆ } be the orthogonalcomplement of the refined Ritz vectors u˜ , i=1,2,...,k with
1 p i
respect to G (A,B;q ,p ). Project QEP (1.1) onto the span{uˆ ,...,uˆ }. Then we obtain a p= m−k-
k 1 1 1 p
dimensional projected QEP and compute its 2p eigenvalues. As a result, the 2p eigenvalues are the
unwanted approximate eigenvalues. We call them refined shift candidates. We select p = m−k refined
shift candidates as shifts, called the refined shifts, for use within the implicitly restarted RGSOAR
algorithm.
8
GONGFangHuietal. SciChinaMathforReview
For GSOAR method, we can take the 2(m−k) unwanted Ritz values as the shift candidates directly.
But it is just possible that one wanted Ritz value and one unwanted Ritz value correspond to the same
Ritz vector. When we take the unwanted Ritz value as shift, that will filter out the wanted Ritz vector
by restarting. To avoid this situation, we can compute the exact shift candidates in the orthogonal
complement of the Ritz vectors y , i = 1,2,...,k with respect to G (A,B;q ,p ). We select p = m−k
i k 1 1
exact shift candidates as shifts, called the exact shifts, for GSOAR method.
Both methods have 2p shift candidates, but if we use more than p shifts, equation (2.16) can not be
truncated at the first k columns. So we can select the p shifts farest from the target as shifts.
3 Shift strategy for implicitly restarted GSOAR and RGSOAR methods
In both methods, all refined vectors u˜ ,...,u˜ and Ritz vectors y ,...,y are approximating the de-
1 k 1 k
sired eigenvectors x ,...,x . The exact and refined shift candidates are approximations to some of
1 k
the unwantedeigenvalues of QEP (1.1) because the information on x ,...,x has been removedfrom
k+1 m
G (A,B;q ,p ). Afterrestartingin(2.17),RGSOARobtainsspan{Q˜ }=span{u˜ ,...,u˜ }andGSOAR
m 1 1 k 1 k
obtains span{Q˜′ }=span{y ,...,y }. Both methods achieve this goal through removing the unwanted
k 1 k
information by shifts. As mentioned above, if two different shifts corresponding to same unwanted
eigenvalue, at least one unwanted eigenvector can not be removed from the subspace G (A,B;q ,p ).
m 1 1
The result of restarting can not be guaranteed. In order to achieve the restarting goal, all shift candi-
dates must be used. The relationship between the number of shifts and the number of zero elements
of bTm = (0 ··· 0 bm−p ···bm−1 bm) in (2.16) is shown in Figure 1. If the number of shifts p is
more than m−k, the first k columns of (2.16) will not be a GSOAR decomposition. Because the last
term in the right side is not only with eT but also with eT , eT until eT . So the implicit restarting
k k−1 k−2 m−p
technique can not be applied directly. Our main work is to propose an implicitly restarted generalized
second-order Arnoldi strategy which can use all shift candidates.
It was shown in [14], that Krylov decomposition and Arnoldi decomposition are equivalent. If we
can take Figure 1 (d) as a generalized Krylov decomposition, then it can be transformed to a GSOAR
decomposition by orthogonal transformation. The scheme is written as the following result.
(a) GSOARdecomposition (b) p=m−k
(c) p>m−k (d) p>m
Figure1 ImplicitlyrestartedGSOARmethodschematic
Theorem 3.1. For equation (2.16), there exists an orthonormal matrix W, that transforms bT into
m
αeT, where α is a non-zero constant, and WTT˜ W is a Hessenberg matrix.
m m
9
GONGFangHuietal. Sci ChinaMathforReview
Proof. Let W be a Householder matrix, such that
1
W b =αe , (3.1)
1 m m
with α as a non-zero constant. Transpose of above equation,
bTWT =αeT, (3.2)
m 1 m
with W =WT, we have
1 1
bTW =αeT. (3.3)
m 1 m
Then multiply matrix W on the right side of equation (2.16)
1
Q Q q
H m W = m W (WTT W )+t m+1 bTW . (3.4)
1 1 1 m 1 m+1,m m 1
"Pm# "Pm# "pm+1#
e e
e
Q e Q e
Take m W as a new m and denote WTT W by new matrix B . Note that B loses the
1 m 1 m m
"Pm# "Pm#
e
Hessenberg form. Combining the above two equations, we have
e
e
Q Q q
H m = m B +b m+1 eT, (3.5)
m m+1,m m
"Pm# "Pm# "pm+1#
where b =t α. At this moment, if we want to do implicit restarting, we must transform B
m+1,m m+1,m m
into Hessenberg again and maintain the last term eT.
m
Householdertransformationandorthogonaltransformationtechniquesareadoptedonceagaintotrans-
formB intoupperHessenbergmatrix. Topreventthenexttransformationdamagingthepreviousresult,
m
the Householder transformations are performed on B from the bottom row to the top row.
m
When l=1, denote B by matrix block form:
m
b b ··· b
11 12 1m
b21 b22 ··· b2m B11 B12
Bm = ... ... ... ... =" c1 bmm#, (3.6)
bm1 bm2 ··· bmm
where c1 = bm,1 ··· bm,m−1 , B11 = Bm(1 : m − 1,1 : m −1) is a (m −1)×(m − 1) matrix,
B12 =Bm(1(cid:16):m−1,m) is a colum(cid:17)n vector of length m−1. Let R1 be a (m−1)×(m−1) Householder
matrix such that
c R =α eT , (3.7)
1 1 1 m−1
R 0
1
where α is a non-zero constant. Let U = , we perform orthogonaltransformation on B ,
1 1 m
" 0 1#
B(1) =UTB U , (3.8)
m 1 m 1
then B(1) has the following structure
m
B(1) B(1)
11 12
B B
U1T 11 12 U1 = c2 , (3.9)
" c1 bmm# 0 ··· 0 B2(12)
where B(1) is a (m−2)×(m−2) matrix, B(1) is a (m−2)×2 matrix, B(1) is a 2×2 matrix and
11 12 22
c2 = bm−1,1 ··· bm−1,m−2 .
(cid:16) (cid:17)
10
GONGFangHuietal. SciChina Mathfor Review
When l=2, we can build a (m−2)×(m−2) Householder matrix R , such that
2
c R =α eT , (3.10)
2 2 2 m−2
R 0
where α is a non-zero constant and use U = 2 to transform B(1) into
2 2 m
" 0 I2×2#
B(2) B(2)
B˜(1) B˜(1) 11 12
11 12
Bm(2) =U2TBm(1)U2 =0 ··· α2 =0 ·c·3· 0 B(2), (3.11)
0 ··· 0 B(1) 22
22 0 ··· 0
where
∗ ∗,∗
B(2) = α (3.12)
22 2 (1)
B
22
0
is a 3×3 Hessenberg matrix.
Suppose B(k−1) is the following matrix block form:
m
(k−1) (k−1)
B B
11 12
c
k
B(k−1) =0 ··· 0 , (3.13)
m
... ··· ... B(k−1)
22
0 ··· 0
whereB(k−1) is a(m−k)×(m−k)matrix, B(k−1) is a(m−k)×k matrix, B(k−1) isa k×k Hessenberg
11 12 22
matrix and ck = bm−k+1,1 ··· bm−k+1,m−k . We can build a (m−k)×(m−k) Householder matrix
Rk, such that (cid:16) (cid:17)
c R =α eT , (3.14)
k k k m−k
R 0
where α is a non-zero constant. Let U = k , we obtain B(k):
k k m
" 0 Ik×k#
B(k) B(k)
B˜(k−1) B˜(k−1) 11 12
11 12
c
k+1
0 ··· αk
Bm(k) =UkTBm(k−1)Uk =0... ······ 0... B2(k2−1)=00... ········· 00... B2(k2), (3.15)
with B(k) is a (m −k −1)×(m −k −1) matrix, B(k) is a (m −k −1)×(k +1) matrix, B(k) is a
11 12 22
(k+1)×(k+1) Hessenberg matrix and ck = bm−k,1 ··· bm−k,m−k−1 .
Obviously, when l=m−2, B(m−2) is an up(cid:16)per Hessenberg matrix. (cid:17)
m
Let U = U U ···U , obviously U is an orthonormal matrix as U ,U ,··· ,U are orthonormal
1 2 m−2 1 1 m−2
matrices. Multiply U on equation (3.5) from right:
Q Q q
H m U = m U(UTB U)+b m+1 eTU. (3.16)
m m+1,m m
"Pm# "Pm# "pm+1#
