The Quadratic Eigenvalue Problem

TL;DR: This work surveys the quadratic eigenvalue problem, treating its many applications, its mathematical properties, and a variety of numerical solution techniques.

Abstract: We survey the quadratic eigenvalue problem, treating its many applications, its mathematical properties, and a variety of numerical solution techniques. Emphasis is given to exploiting both the structure of the matrices in the problem (dense, sparse, real, complex, Hermitian, skew-Hermitian) and the spectral properties of the problem. We classify numerical methods and catalogue available software.

Summary

2. Applications of QEPs.

• A wide variety of applications require the solution of a QEP, most of them arising in the dynamic analysis of structural mechanical, and acoustic systems, in electrical circuit simulation, in fluid mechanics, and, more recently, in modeling microelectronic mechanical systems (MEMS) [28] , [155] .
• QEPs also have interesting applications in linear algebra problems and signal processing.
• The list of applications discussed in this section is by no means exhaustive, and, in fact, the number of such applications is constantly growing as the methodologies for solving QEPs improve.

C

• Figure 2 .2 illustrates the behavior of the solution for different values of the viscous damping factor ζ.
• The properties of systems of second-order differential equations (2.1) for n ≥ 1 have been analyzed in some detail by Lancaster [87] and more recently by Gohberg, Lancaster, and Rodman [62] .
• If an eigenvalue has positive real part, q(t) can grow exponentially toward infinity.
• For the Millennium Bridge it is unwanted.
• By this means, eigenvalues corresponding to unstable modes or yielding large vibrations can be relocated or damped.

2.2. Vibration Analysis of Structural Systems-Modal Superposition

• The futuristic footbridge, spanning the Thames by St Paul's Cathedral, will be fitted with shock absorbers to reduce the alarming 10cm swaying first noticed when a huge crowd surged across at its opening in June.
• The bridge, designed by Lord Foster and closed three days after its opening, will be fitted with "viscous dampers" and "tuned mass dampers"-likened to shock absorbers-by the bridge's engineers, Ove Arup.

2.3. Vibro-Acoustics.

• This application raises the problem of truncation error due to the discretization of operators in an infinite-dimensional space.
• The authors will not treat this aspect in this paper but will restrict ourselves to the study of QEPs of finite dimension.

2.4. Fluid Mechanics.

• When studying the temporal stability, β is a real and fixed wavenumber.
• The imaginary part of ω represents the temporal growth rate.

3.8.1. Spectrum Location.

• Symmetry alone does not guarantee that the eigenvalues are real and that the pencil is diagonalizable by congruence transformations.
• Veselić [147] shows that the overdamping condition (3.17) is equivalent to the definiteness of the symmetric linearization obtained from L1 or L2.
• Hence, for overdamped problems (see section 3.9), the symmetric linearization fully reflects the spectral properties of the QEP (real eigenvalues).

3.8.3. Self-Adjoint Triple, Sign Characteristic, and Factorization.

• All the eigenvalues lie in the left half-plane as shown in Figure 3 .3.
• Figure 3 .4 displays the eigenvalue distribution, with the characteristic gap between the n smallest eigenvalues and the n largest.

4.1. Conditioning.

• Condition numbers for the eigenvalues can be derived without assuming that the eigenvalue is finite.
• Moreover, by working in projective spaces the problem of choosing a normalization for the eigenvectors is avoided.

4.2. Backward Error. A natural definition of the normwise backward error of an approximate eigenpair (

• It is often not appreciated that in finite precision arithmetic these two choices are not equivalent (see the example in section 5.1).
• One approach, adopted by MATLAB 6's polyeig function for solving the polynomial eigenvalue problem and illustrated in Algorithm 5.1 below, is to use whichever part of ξ yields the smallest backward error (4.7).

4.3. Pseudospectra.

• The eigenvalues of the unperturbed QEP, marked by dots, all lie in the left half-plane, so the unperturbed system is stable.
• The plot shows that relative perturbations of order 10 −10 (corresponding to the yellow region of the plot) can move the eigenvalues to the right half-plane, making the system unstable.

5. Numerical Methods for Dense Problems.

• Is the backward error of the GEP solution and, as expected, is of the order of the unit roundoff because the QZ algorithm is a backward stable algorithm for the solution of the GEP.
• For the smallest eigenvalue in modulus, the second choice ξ 2 yields a smaller backward error.
• This example shows that even if the algorithm chooses the part of the vector ξ that yields the smallest backward error, these backward errors can be much larger than the unit roundoff.
• Even though backward stability is not guaranteed, Algorithm 5.1 is the preferred method if M , C, and K have no particular structure and are not too large.

5.2. Symmetric Linearization.

• If the pencil is indefinite, the HR [21] , [23] , LR [124] , and Falk-Langemeyer [45] algorithms can be employed in order to take advantage of the symmetry of the GEP, but all these methods can be numerically unstable and can even break down completely.
• There is a need to develop efficient and reliable numerical methods for the solution of symmetric indefinite GEPs.

5.3. Hamiltonian/Skew-Hamiltonian

• Figure 5 .2 compares the spectrum of the QEP associated with (5.3) when computed by the MATLAB 6 function polyeig and their implementation of Van Loan's square-reduced algorithm [146] .
• Not surprisingly, the eigenvalues computed by polyeig (which uses a companion linearization and the QZ algorithm) do not have a Hamiltonian structure, and some of them have positive real part, suggesting incorrectly that the system described by (5.3) is unstable.
• In contrast, Van Loan's squarereduced algorithm together with the Hamiltonian/skew-Hamiltonian linearization preserves the Hamiltonian structure, yielding pure imaginary computed eigenvalues that confirm the system's stability.

5.4. Sensitivity of the Linearization.

• This leads to the Lanczos method when S is Hermitian and the Arnoldi or non-Hermitian Lanczos method when S is non-Hermitian.
• The authors aim is to explain how existing Krylov subspace methods can be applied to the QEP.

6.2. Projection Methods Applied

• The goal is to add a new vector v k+1 so that the Ritz pair becomes more accurate.
• The subspace is extended by a quadratic Cayley transform applied to the Ritz vector.

7.3. Hamiltonian/Skew-Hamiltonian Linearization-Direct Methods.

• As the authors have stressed throughout this survey, QEPs are of growing importance and many open questions are associated with them.
• The authors hope that their unified overview of applications, theory, and numerical methods for QEPs will stimulate further research in this area.

Figures

Fig. 2.1 Left: Single mass-spring system. Right: A simple electric circuit.

Fig. 2.2 Response of the single degree of freedom mass-spring system of Figure 2.1 for different values of the viscous damping factor ζ = C/2Mω and initial conditions q(0) = 0, q̇(0) = 10.

Fig. 2.4 Spacecraft with deployable antenna and two solar panels attached to a rigid base undergoing angular motion of angular velocity Ω.

Table 5.2 Comparison, for three different linearizations, of the condition number and relative error of a small and a large eigenvalue in absolute value of the QEP associated with the massspring problem.

Fig. 3.2 An n-degrees-of-freedom damped mass-spring system.

Fig. 3.5 Location of the spectrum of the gyroscopic system (3.21) as the gyroscopic parameter Ω increases.

Table 5.1 Nuclear power plant problem. Backward errors of the computed eigenpairs corresponding to the smallest and largest eigenvalues in modulus.

Fig. 2.3 Nuclear power plant simplified into an eight-degrees-of-freedom system, from [79]. The system is composed of four elastically interconnected different types of rigid structures: core, prestressed concrete pressure vessel (PCPV), building, and basement. Each structure has two degrees of freedom that correspond to sway direction u and rocking direction φ.

Fig. 1.1 The Millennium footbridge over the river Thames.

Fig. 4.1 Spectral portrait of the eight-degrees-of-freedom nuclear power plant system illustrated in Figure 2.3. The bar shows the mapping of colors to log10 of the reciprocal of the scaled resolvent norm.

Fig. 3.1 Spectrum of the eight-degrees-of-freedom simplified nuclear power plant system illustrated in Figure 2.3. Left: M , C, and K are real. Right: Hysteretic damping is added so that K is complex.

Fig. 3.5 Location of the spectrum of the gyroscopic system (3.21) as the gyroscopic parameter Ω increases.

Fig. 2.6 Coordinate system for attachment-line boundary layer flow.

Table 1.1 Matrix properties of QEPs considered in this survey with corresponding spectral properties. The first column refers to the section where the problem is treated. Properties can be added: QEPs for which M , C, and K are real symmetric have properties P3 and P4 so that their eigenvalues are real or come in pairs (λ, λ̄) and the sets of left and right eigenvectors coincide. Overdamped systems yield QEPs having the property P6. In P6, γ(M,C,K) = min{(x∗Cx)2 − 4(x∗Mx)(x∗Kx) : ||x||2 = 1}. Gyroscopic systems yield QEPs having at least the property P7.

Fig. 5.1 A schematic view of a traveling band.

Fig. 2.5 Fluid in a cavity with one absorbing wall.

Fig. 3.4 Eigenvalue distribution of the QEP for the overdamped mass-spring system with n = 50. Now all the eigenvalues are real, so they are plotted against the index rather than as points in the complex plane.

Fig. 3.3 Eigenvalues in the complex plane of the QEP for the nonoverdamped mass-spring system with n = 50.

Fig. 5.2 Spectrum of the moving band illustrated in Figure 5.1. polyeig uses a companion linearization and the QZ algorithm. srmh uses a Hamiltonian linearization and the square-reduced algorithm for the Hamiltonian eigenproblem.

Full text

The quadratic eigenvalue problem
Tisseur, Françoise and Meerbergen, Karl
2001
MIMS EPrint: 2006.256
Manchester Institute for Mathematical Sciences
School of Mathematics
The University of Manchester
Reports available from: http://eprints.maths.manchester.ac.uk/
And by contacting: The MIMS Secretary
School of Mathematics
The University of Manchester
Manchester, M13 9PL, UK
ISSN 1749-9097

SIAM REVIEW
c

2001 Society for Industrial and Applied Mathematics
Vol. 43, No. 2, pp. 235–286
The Quadratic Eigenvalue
Problem
∗
Fran¸coise Tisseur
†
Karl Meerbergen
‡
Abstract. We survey the quadratic eigenvalue problem, treating its many applications, its mathe-
matical properties, and a variety of numerical solution techniques. Emphasis is given to
exploiting both the structure of the matrices in the problem (dense, sparse, real, com-
plex, Hermitian, skew-Hermitian) and the spectral properties of the problem. We classify
numerical methods and catalogue available software.
Key words. quadratic eigenvalue problem, eigenvalue, eigenvector, λ-matrix, matrix polynomial,
second-order differential equation, vibration, Millennium footbridge, overdamped sys-
tem, gyroscopic system, linearization, backward error, pseudospectrum, condition num-
ber, Krylov methods, Arnoldi method, Lanczos method, Jacobi–Davidson method
AMS subject classification. 65F30
PII. S0036144500381988
1. Introduction. On its opening day in June 2000, the 320-meter-long Millen-
nium footbridge over the river Thames in London (see Figure 1.1) started to wobble
alarmingly under the weight of thousands of people; two days later the bridge was
closed. To explain the connection between this incident and the quadratic eigenvalue
problem (QEP), the subject of this survey, we need to introduce some ideas from
vibrating systems. A natural frequency of a structure is a frequency at which the
structure prefers to vibrate. When a structure is excited by external forces whose
frequencies are close to the natural frequencies, the vibrations are amplified and
the system becomes unstable. This is the phenomenon of resonance. The exter-
nal forces on the Millennium Bridge were pedestrian-induced movements. On the
opening day, a large crowd of people initially walking randomly started to adjust
their balance to the bridge movement, probably due to high winds on that day. As
they walked, they became more synchronized with each other and the bridge started
to wobble even more. The lateral vibrations experienced on the bridge occurred
because some of its natural modes of vibration are similar in frequency to the side-
ways component of pedestrian footsteps on the bridge.
1
The connection with this
survey is that the natural modes and frequencies of a structure are the solution of
∗
Received by the editors November 3, 2000; accepted for publication (in revised form) December
4, 2000; published electronically May 2, 2001.
http://www.siam.org/journals/sirev/43-2/38198.html
†
Department of Mathematics, University of Manchester, Manchester M13 9PL, England (ftisseur@
ma.man.ac.uk, http://www.ma.man.ac.uk/˜ftisseur/). This author’s work was supported by Engi-
neering and Physical Sciences Research Council grant GR/L76532.
‡
Free Field Technologies, 16, place de l’Universit´e, 1348 Louvain-la-Neuve, Belgium (Karl.
Meerbergen@fft.be).
1
More details on the Millennium Bridge, its construction, and the wobbling problem can be found
at http://www.arup.com/MillenniumBridge.
235

236 FRANC¸ OISE TISSEUR AND KARL MEERBERGEN
Fig. 1.1 The Millennium footbridge over the river Thames.
an eigenvalue problem that is quadratic when damping effects are included in the
model.
The QEP is currently receiving much attention because of its extensive applica-
tions in areas such as the dynamic analysis of mechanical systems in acoustics and
linear stability of flows in fluid mechanics. The QEP is to find scalars λ and nonzero
vectors x, y satisfying
(λ
2
M + λC + K)x =0,y
∗
(λ
2
M + λC + K)=0,
where M, C, and K are n ×n complex matrices and x, y are the right and left eigen-
vectors, respectively, corresponding to the eigenvalue λ. A major algebraic difference
between the QEP and the standard eigenvalue problem (SEP),
Ax = λx,
and the generalized eigenvalue problem (GEP),
Ax = λBx,
is that the QEP has 2n eigenvalues (finite or infinite) with up to 2n right and 2n left
eigenvectors, and if there are more than n eigenvectors they do not, of course, form a
linearly independent set.
QEPs are an important class of nonlinear eigenvalue problems that are less fa-
miliar and less routinely solved than the SEP and the GEP, and they need special
attention. A major complication is that there is no simple canonical form analogous
to the Schur form for the SEP or the generalized Schur form for the GEP.
A large literature exists on QEPs, spanning the range from theory to applications,
but the results are widely scattered across disciplines. Our goal in this work is to
gather together a collection of applications where QEPs arise, to identify and classify
their characteristics, and to summarize current knowledge of both theoretical and
algorithmic aspects. The type of problems we will consider in this survey and their
spectral properties are summarized in Table 1.1.
The structure of the survey is as follows. In section 2 we discuss a number
of applications, covering the motivation, characteristics, theoretical results, and algo-
rithmic issues. Section 3 reviews the spectral theory of QEPs and discusses important
classes of QEPs coming from overdamped systems and gyroscopic systems. Several

THE QUADRATIC EIGENVALUE PROBLEM 237
Tab le 1 . 1 Matrix properties of QEPs considered in this survey with corresponding spectral proper-
ties. The first column refers to the section where the problem is treated. Properties can
be added: QEPs for which M, C, and K are real symmetric have properties P 3 and P 4
so that their eigenvalues are real or come in pairs (λ,
¯
λ) and the sets of left and right
eigenvectors coincide. Overdamped systems yield QEPs having the property P 6.InP 6,
γ(M, C,K) = min{(x
∗
Cx)
2
− 4(x
∗
Mx)(x
∗
Kx):||x||
2
=1}. Gyroscopic systems yield
QEPs having at least the property P 7.
Matrix properties Eigenvalue properties Eigenvector properties
P1
§3.1
M nonsingular 2n finite eigenvalues
P2
§3.1
M singular
Finite and infinite eigen-
values
P3
§3.1
M, C, K real Eigenvalues are real or
come in pairs (λ,
¯
λ)
If x is a right eigenvector of λ
then ¯x is a right eigenvector
of
¯
λ
P4
§3.8
M, C, K Hermitian Eigenvalues are real or
come in pairs (λ,
¯
λ)
If x is a right eigenvector of
λ then x is a left eigenvector
of
¯
λ
P5
§3.8
M Hermitian positive
definite, C, K Hermitian
positive semidefinite
Re(λ) ≤ 0
P6
§3.9
M, C symmetric positive
definite, K symmetric
positive semidefinite,
γ(M, C,K) > 0
λs are real and negative,
gap between n largest and
n smallest eigenvalues
n linearly independent
eigenvectors associated with
the n largest (n smallest)
eigenvalues
P7
§3.10
M, K Hermitian,
M positive definite,
C = −C
∗
Eigenvalues are purely
imaginary or come in
pairs (λ, −
¯
λ)
If x is a right eigenvector of
λ then x is a left eigenvector
of −
¯
λ
P8
§3.10
M, K real symmetric
and positive definite,
C = −C
T
Eigenvalues are purely
imaginary
linearizations are introduced and their properties discussed. In this paper, the term
linearization means that the nonlinear QEP is transformed into a linear eigenvalue
problem with the same eigenvalues. Our treatment of the large existing body of the-
ory is necessarily selective, but we give references where further details can be found.
Section 4 deals with tools that offer an understanding of the sensitivity of the problem
and the behavior of numerical methods, covering condition numbers, backward errors,
and pseudospectra.
A major division in numerical methods for solving the QEP is between those
that treat the problem in its original form and those that linearize it into a GEP of
twice the dimension and then apply GEP techniques. A second division is between
methods for dense, small- to medium-size problems and iterative methods for large-
scale problems. The former methods compute all the eigenvalues and are discussed
in section 5, while the latter compute only a selection of eigenvalues and eigenvectors
and are reviewed in section 6. The proper choice of method for particular problems
is discussed. Finally, section 7 catalogues available software, and section 8 contains
further discussion and related problems.
We shall generally adopt the Householder notation: capital letters A, B, C, ...
denote matrices, lower case roman letters denote column vectors, Greek letters denote
scalars, ¯α denotes the conjugate of the complex number α, A
T
denotes the transpose
of the complex matrix A, A
∗
denotes the conjugate transpose of A, and ·is any
vector norm and the corresponding subordinate matrix norm. The values taken by

238 FRANC¸ OISE TISSEUR AND KARL MEERBERGEN
any integer variable are described using the colon notation: “i =1:n” means the
same as “i =1, 2,...,n.” We write A>0(A ≥ 0) if A is Hermitian positive definite
(positive semidefinite) and A<0(A ≤ 0) if A is Hermitian negative definite (negative
semidefinite). A definite pair (A, B) is defined by the property that A, B ∈ C
n×n
are
Hermitian and
min
z∈C
n
z
2
=1

(z
∗
Az)
2
+(z
∗
Bz)
2
> 0,
which is certainly true if A>0orB>0.
2. Applications of QEPs. A wide variety of applications require the solution
of a QEP, most of them arising in the dynamic analysis of structural mechanical,
and acoustic systems, in electrical circuit simulation, in fluid mechanics, and, more
recently, in modeling microelectronic mechanical systems (MEMS) [28], [155]. QEPs
also have interesting applications in linear algebra problems and signal processing.
The list of applications discussed in this section is by no means exhaustive, and, in
fact, the number of such applications is constantly growing as the methodologies for
solving QEPs improve.
2.1. Second-Order Differential Equations. To start, we consider the solution
of a linear second-order differential equation
M ¨q(t)+C ˙q(t)+Kq(t)=f(t),(2.1)
where M, C, and K are n × n matrices and q(t)isannth-order vector. This is the
underlying equation in many engineering applications. We show that the solution can
be expressed in terms of the eigensolution of the corresponding QEP and explain why
eigenvalues and eigenvectors give useful information.
Two important areas where second-order differential equations arise are the fields
of mechanical and electrical oscillation. The left-hand picture in Figure 2.1 illustrates
a single mass-spring system in which a rigid block of mass M is on rollers and can
move only in simple translation. The (static) resistance to displacement is provided
by a spring of stiffness K, while the (dynamic) energy loss mechanism is represented
by a damper C; f(t) represents an external force. The equation of motion governing
this system is of the form (2.1) with n =1.
As a second example, we consider the flow of electric current in a simple RLC
circuit composed of an inductor with inductance L, a resistor with resistance R, and
a capacitor with capacitance C, as illustrated on the right of Figure 2.1; E(t)is
the input voltage. The Kirchhoff loop rule requires that the sum of the changes in
potential around the circuit must be zero, so
L
di(t)
dt
+ Ri(t)+
q(t)
C
− E(t)=0,(2.2)
where i(t) is the current through the resistor, q(t) is the charge on the capacitor, and
t is the elapsed time. The charge q(t) is related to the current i(t)byi(t)=dq(t)/dt.
Differentiation of (2.2) gives the second-order differential equation
L
d
2
i(t)
dt
2
+ R
di(t)
dt
+
1
C
i(t)=
dE(t)
dt
,(2.3)
which is of the same form as (2.1), again with n =1.

Frequently asked questions

Q1. What can be done to solve the projected problem?

As an orthogonal projection preserves symmetry and skew-symmetry, the numerical methods of sections 5.2 and 5.3 can be used to solve the projected problem. 

Q2. What is the projection method for L()?

The projection method approximates an eigenvector x of L(λ) by a vector x̃ = V ξ ∈ Kk with corresponding approximate eigenvalue λ̃. AsW ∗L(λ̃)x̃ = W ∗L(λ̃)V ξ = Lk(λ̃)ξ = 0, the projection method forces the residual r = L(λ̃)x̃ to be orthogonal to Lk. 

Q3. What is the eigenvalue of the skew-Hamiltonian matrix?

The reduction to a Hamiltonian eigenproblem uses the fact that when the skew-Hamiltonian matrix B is nonsingular, it can be written in factored form asB = B1B2 = [ The author12C 0 M ] [ M 12C 0 The author] with BT2 J = JB1.(5.2)Then H = B−11 AB −1 2 is Hamiltonian. 

Q4. What is the simplest way to construct a Krylov subspace?

It produces a non-Hermitian tridiagonal matrix Tk and a pair of matrices Vk and Wk such that W ∗kVk = The authorand whose columns form bases for the Krylov subspaces Kk(S, v) and Kk(S∗, w), where v and w are starting vectors such that w∗v = 1. 

Q5. What is the simplest way to measure the robustness of a numerical algorithm?

To measure the robustness of the system, one can take as a global measureν2 = 2n∑ k=1 ω2kκ(λk) 2,where the ωk are positive weights. 

Q6. What is the advantage of working with a linearization of Q() and a?

One important advantage of working with a linearization of Q(λ) and a Krylov subspace method is that one can get at the same time the partial Schur decomposition of the single matrix S that is used to define the Krylov subspaces. 

Q7. What is the spectral transformation used to approximate eigenvalues?

The shift-and-invert spectral transformation f(λ) = 1/(λ − σ) and the Cayley spectral transformation f(λ) = (λ− β)/(λ− σ) (for β = σ), used to approximate eigenvalues λ closest to the shift σ, are other possible spectral transformations that are discussed in [7], for example. 

Q8. What is the way to deflate a Ritz eigenpair?

If (λ̃, x̃) is a converged Ritz eigenpair that belongs to the set of desired eigenvalues, one may want to lock it and then continue to compute the remaining eigenvalues without altering (λ̃, x̃). 

Q9. What is the main drawback of the pseudo-Lanczos algorithm?

Parlett and Chen [114] introduced a pseudo-Lanczos algorithm for symmetric pencils that uses an indefinite inner product and respects the symmetry of the problem. 

Q10. What is the distribution of the eigenvalues of G() in the complex?

G(λ)∗ = G(−λ̄),(3.19)the distribution of the eigenvalues of G(λ) in the complex plane is symmetric with respect to the imaginary axis. 

Q11. What is the common method for finding eigenpairs of large QEPs?

Most algorithms for large QEPs proceed by generating a sequence of subspaces {Kk}k≥0 that contain increasingly accurate approximations to the desired eigenvectors. 

Q12. What is the way to solve the eigenvalue problem?

One approach, adopted by MATLAB 6’s polyeig function for solving the polynomial eigenvalue problem and illustrated in Algorithm 5.1 below, is to use whichever part of ξ̃ yields the smallest backward error (4.7). 

Q13. What can be done to solve the symmetric indefinite GEP?

If the pencil is indefinite, the HR [21], [23], LR [124], and Falk–Langemeyer [45] algorithms can be employed in order to take advantage of the symmetry of the GEP, but all these methods can be numerically unstable and can even break down completely. 

Q14. What is the disadvantage of the Arnoldi and Lanczos methods?

(Note that the matrix Hk is the Galerkin projection of S and not of A− λB.)A major disadvantage of the shift-and-invert Arnoldi and Lanczos methods is that a change of shift σ requires building a new Krylov subspace: all information built with the old σ is lost. 

Q15. What is the way to deflate a converged eigenvalue?

If the converged (λ̃, x̃) does not belong to the set of wanted eigenvalues, one may want to remove it from the current subspace Kk. 

Q16. What is the oblique projection of L() onto Kk?

In this case Lk is the orthogonal projection of L(λ) onto Kk. When W = V , Lk is the oblique projection of L(λ) onto Kk along Lk. 

Q17. What is the reorthogonalization of columns of Vk?

This process does not guarantee orthogonality of the columns of Vk in floating point arithmetic, so reorthogonalization is recommended to improve the numerical stability of the method [32], [135].