Stability analysis of improved Two-level orthogonal Arnoldi procedure

TL;DR

This paper introduces the I-TOAR method, an improved version of the Two-level orthogonal Arnoldi procedure, and provides a stability analysis addressing an open problem related to the numerical stability of second-order Krylov subspace computations.

Contribution

The paper proposes the I-TOAR method and solves the open problem of its stability analysis with respect to quadratic problem coefficients.

Findings

I-TOAR enhances numerical stability over previous methods.

The stability analysis confirms robustness of I-TOAR.

Results demonstrate improved accuracy in second-order Krylov subspace computations.

Abstract

The SOAR method for computing an orthonormal basis of a second-order Krylov subspace can be numerically unstable (see Lu et al. (2016)). In the Two-level orthogonal Arnoldi(TOAR) procedure, an alternative to SOAR, the problem of instability had circumvented. A stability analysis of the second-order Krylov subspace's orthonormal basis in TOAR with respect to the coefficient matrices of a quadratic problem remain open; see Lu et al. (2016). This paper proposes the Improved-TOAR method(I-TOAR) and solves the said open problem for I-TOAR.

Equations

$$r_j=A{r}_{j-1}+B{r}_{j-2}\text{mbox}forj\ge 1$$

$$G_k(A,B;r_{-1},r_0)\equiv \text{mbox}span\{r_{-1},r_0,r_1,\cdots ,r_{k-1}\}$$

$$L=\left[\begin{array}{cc}A& B\\ I& 0\end{array}\right]\text{mbox}and{v}_0=\left(\begin{array}{c}{r}_0\\ r_{-1}\end{array}\right),$$

$$A{Q}_k+B{P}_k=Q_k{T}_k+t_{k+1,k}{q}_{k+1}{e}_k^{\ast },$$

$$Q_k=P_k{T}_k+t_{k+1,k}{p}_{k+1}{e}_k^{\ast },$$

$$\left[\begin{array}{cc}A& B\\ I& 0\end{array}\right]\left[\begin{array}{c}{Q}_k\\ P_k\end{array}\right]=\left[\begin{array}{c}{Q}_k\\ P_k\end{array}\right]H_k+t_{k+1,k}\left(\begin{array}{c}{q}_{k+1}\\ p_{k+1}\end{array}\right)e_k^{\ast }.$$

$$[r_{-1}{r}_0]=Q_{1}X,$$

$$U_{1,1}:=X(:,2)/\Vert v_0\Vert \text{mbox}and{U}_{1,2}:=X(:,1)/\Vert v_0\Vert .$$

$$L{V}_k=V_{k+1}\underline{H_{k+1}},$$

$$\text{mbox}span\{Q_{k+1}\}=\text{mbox}span\{Q_k,r\},$$

$$Q_{k+1}=Q_k,U_{k+1,1}=[U_{k,1}{x}_k],and{U}_{k+1,2}=[U_{k,2}{y}_k];$$

$$Q_{k+1}=[Q_k{q}_{k+1}],U_{k+1,1}=\left[\begin{array}{cc}{U}_{k,1}& x_k\\ 0& {\beta }_k\end{array}\right],and{U}_{k+1,2}=\left[\begin{array}{cc}{U}_{k,2}& y_k\\ 0& 0\end{array}\right].$$

$$\left[\begin{array}{c}{Q}_k{U}_{k,1}\\ Q_k{U}_{k,2}\end{array}\right]=\left[\begin{array}{cc}{Q}_k& \\ & Q_k\end{array}\right]\left[\begin{array}{c}{U}_{k,1}\\ U_{k,2}\end{array}\right].$$

$$\left[\begin{array}{cc}A& B\\ I& 0\end{array}\right]\left[\begin{array}{c}{Q}_k{U}_{k,1}\\ Q_k{U}_{k,2}\end{array}\right]=\left[\begin{array}{c}{Q}_{k+1}{U}_{k+1,1}\\ Q_{k+1}{U}_{k+1,2}\end{array}\right]\underline{H_{k+1}}.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }(A{Q}_k{U}_{k,1}+B{Q}_k{U}_{k,2})+U_{k,2}^{\ast }(Q_k^{\ast }{Q}_k)U_{k,1}=H_k,$$

$$Q_k{U}_{k,1}=Q_{k+1}{U}_{k+1,2}\underline{H_{k+1}}.$$

$$U_{k,1}=[I0]U_{k+1,2}\underline{H_{k+1}}.$$

$$U_{k,1}=U_{k,2}{H}_k+h_{k+1,k}{y}_k{e}_k^{\ast }$$

$$U_{k,2}^{\ast }{U}_{k,1}=U_{k,2}^{\ast }{U}_{k,2}{H}_k+h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }(A{Q}_k{U}_{k,1}+B{Q}_k{U}_{k,2})+U_{k,2}^{\ast }{U}_{k,2}{H}_k+h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }=H_k.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }(A{Q}_k{U}_{k,1}+B{Q}_k{U}_{k,2})-U_{k,1}^{\ast }{U}_{k,1}{H}_k=-h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }{Q}_{k+1}{U}_{k+1,1}\underline{H_{k+1}}-U_{k,1}^{\ast }{U}_{k,1}{H}_k=-h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }.$$

$$U_{k,1}^{\ast }[I0]U_{k+1,1}\underline{H_{k+1}}-U_{k,1}^{\ast }{U}_{k,1}{H}_k=-h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }.$$

$$U_{k,1}^{\ast }[U_{k,1}{x}_k]\underline{H_{k+1}}-U_{k,1}^{\ast }{U}_{k,1}{H}_k=-h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }.$$

$$U_{k,1}^{\ast }{U}_{k,1}{H}_k+h_{k+1,k}{U}_{k,1}^{\ast }{x}_k{e}_k^{\ast }-U_{k,1}^{\ast }{U}_{k,1}{H}_k=-h_{k+1,k}{U}_{k,2}^{\ast }{y}_k{e}_k^{\ast }.$$

$$Q_k^{\ast }(A{Q}_k{U}_{k,1}+B{Q}_k{U}_{k,2})=U_{k,1}{H}_k+h_{k+1,k}{x}_k{e}_k^{\ast }.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }(A{Q}_k{U}_{k,1}+B{Q}_k{U}_{k,2})=U_{k,1}^{\ast }{U}_{k,1}{H}_k.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }A{Q}_k{U}_{k,1}+{\underline{H_{k+1}}}^{\ast }{U}_{k+1,2}^{\ast }{Q}_{k+1}^{\ast }B{Q}_k{U}_{k,2}=U_{k,1}^{\ast }{U}_{k,1}{H}_k.$$

$${\underline{H_{k+1}}}^{\ast }{U}_{k+1,2}^{\ast }{Q}_{k+1}^{\ast }B{Q}_k{U}_{k,2}=H_k^{\ast }{U}_{k,2}^{\ast }{Q}_k^{\ast }B{Q}_k{U}_{k,2}.$$

$$U_{k,1}^{\ast }{Q}_k^{\ast }A{Q}_k{U}_{k,1}=U_{k,1}^{\ast }{U}_{k,1}{H}_k-H_k^{\ast }{U}_{k,2}^{\ast }{Q}_k^{\ast }B{Q}_k{U}_{k,2}.$$

Full text

1 Introduction

Large-scale quadratic problems are ubiquitous in Linear stability analysis, Model order reduction, Dissipative acoustics, and Constraint least squares problems; See Singh et al. (2016); Bermudez et al. (2000); Gander et al. (1989); Sima et al. (2004), and Huitfeldt & Ruhe (1990). We suggest the readers refer Tisseur & Meerbergen (2001) for other applications. Two approaches are well-known to solve quadratic eigenvalue problems and quadratic system of equations. One approach finds an appropriate linearization that results in linear eigenvalue problems or a linear system of equations. Another approach projects larger sparse quadratic problem onto a lower dimensional subspace, and subsequently produce a small, dense QEP or a system of equations. The first approach has a drawback that it increases the condition number due to linear problems of double the size; see, e.g., Hwang et al. (2003). The popular methods such as Residual inverse iteration, Second-Order Arnoldi (SOAR), and Two-level orthogonal Arnoldi (TOAR) methods follow the second approach; See Bai & Su (2005b); Neumaier (1985); Meerbergen (2001), and Lu et al. (2016).

For the given quadratic problem, The SOAR method constructs an orthonormal basis of a second-order Krylov subspace using a recurrence relation an analog to that in the Arnoldi method. It also generates a non-orthonormal basis of a Krylov subspace associated with the corresponding linear problem. To do this, SOAR requires a solution of a triangular linear system, ill-conditioned, in general. Though this causes numerical instability, the SOAR method found applications in Quadratic eigenvalue problems, Structural acoustics analysis, and Model order reduction of second-order dynamical systems, etc. (see Yang (2005), Puri & Morrey (2013), and Bai & Su (2005a) for further information).

Zhu (2005) and Su et al. (2008) proposed the Two-Level orthogonal Arnoldi(TOAR) method to overcome the instability problem in SOAR. As the name suggests, the Gram-Schmidt orthonormalization procedure involved in the two levels of TOAR; in the first and the second levels to construct orthonormal bases for the second-order Krylov subspace and the associated linear Krylov subspace, respectively. Lu et al. (2016) proved under some mild assumptions that TOAR with partial reorthogonalization is backward stable to compute an orthonormal basis of associated linear Krylov subspace. However, similar stability analysis for the second-order Krylov subspace with respect to the coefficient matrices of a quadratic problem left open; see the Concluding Remarks in Lu et al. (2016).

In this paper, we are proposing the Improved TOAR(I-TOAR) method. The I-TOAR method improves TOAR in constructing an orthonormal basis of a second-order Krylov subspace. The proposed improvement is necessary in the TOAR method to solve the said open problem. Using I-TOAR, this paper does the stability analysis for the second-order Krylov subspace with respect to the coefficient matrices of a quadratic problem.

This paper is organized as follows. In Section-2 we briefly discuss the SOAR and TOAR methods and establish relations between the matrix $Q_{k}$ and submatrices of $U_{k}$ those generated by the TOAR method. Then, Section-3 presents theoretical results those motivated to improve TOAR. Section-4 proposes the I-TOAR(Improved TOAR) algorithm and discusses its implementation details. Section-5 does rigorous backward error analysis of I-TOAR in terms of coefficient matrices in quadratic problems for computing an orthonormal basis of a second-order Krylov subspace. Section-6 compares the results of numerical experiments with an application of the TOAR and I-TOAR methods in the Model Order Reduction of second-order dynamical systems. Section-7 concludes the paper.

2 SOAR and TOAR methods

Let $A$ and $B$ be matrices of order $n,$ and $r_{-1},$ $r_{0}$ be vectors of length $n.$ Then, a sequence of vectors $r_{1},r_{2},r_{3},\cdots,$ satisfying the following recurrence relation,

[TABLE]

is called a second-order Krylov sequence. The subspace

[TABLE]

is called a kth second-order Krylov subspace. A second-order Krylov subspace can be embedded in the linear Krylov subspace $K_{k}(L,v_{0}),$ for

[TABLE]

where $I$ is an Identity matrix of order $n,$ for details refer to Lu et al. (2016).

Let a vector $q_{1}$ be a linear combination of $r_{-1},$ $r_{0},$ and $\|q_{1}\|=1.$ If a vector $q_{k+1}$ at $k^{th}$ iteration of SOAR is non-zero then it is orthogonal to the set of unit vectors generated in previous iterations. Further, $\|q_{k+1}\|=1.$ The non-zero column vectors of $Q_{k+1}\equiv[q_{1}~{}q_{2}~{}\cdots~{}q_{k}~{}q_{k+1}]$ form an orthonormal basis for the second-order Krylov subspace $G_{k}(A,B;r_{-1},r_{0}).$ In SOAR, $P_{k}:=[p_{1},p_{2},\cdots p_{k}]$ is the matrix satisfying the following relations:

[TABLE]

and

[TABLE]

where $T_{k}$ is an upper Hessenberg matrix of order $k.$ In compact form these relations can be written as follows:

[TABLE]

Observe that the above equation is similar to the Arnoldi decomposition for the matrix $L$ and the initial vector $v_{0}.$ However, column vectors of a matrix $\begin{bmatrix}Q_{k}\\ P_{k}\end{bmatrix}$ are non-orthonormal even though they form a basis for the linear Krylov subspace $K_{k}(L,v_{0}).$ The SOAR method requires the solution of ill-conditioned triangular linear system in order to avoid explicit computation of $P_{k}.$ Consequently, this makes the SOAR procedure numerically instable. To circumvent the instability an alternative method proposed in Lu et al. (2016), the Two-level orthogonalization Arnoldi(TOAR) method.

The TOAR method starts with randomly chosen non-zero initial vector $v_{0}:=\begin{pmatrix}r_{0}\\ r_{-1}\end{pmatrix}.$ Then, it finds a rank revealing QR decomposition of the matrix $[r_{-1}~{}~{}r_{0}]$:

[TABLE]

where $Q_{1}$ and $X$ are matrices of order $n\times\eta_{1}$ and $\eta_{1}\times 2,$ respectively. If the vectors $r_{-1}$ and $r_{0}$ are linearly independent then $\eta_{1}=2,$ otherwise $\eta_{1}=1.$ Following the MATLAB notation define

[TABLE]

In TOAR the vector $\begin{pmatrix}Q_{1}U_{1,1}\\ Q_{1}U_{1,2}\end{pmatrix}$ forms an orthonormal basis for $K_{1}(L,v_{0}).$ Henceforth, TOAR recursively computes the matrices $Q_{k},$ $U_{k,1},$ and $U_{k,2}$ using the relations in the following lemma; the Lemma-3.1 in Lu et al. (2016).

Lemma 1**.**

Let column vectors of $V_{j}\equiv\begin{bmatrix}Q_{j}U_{j,1}\\ Q_{j}U_{j,2}\end{bmatrix}$ form an orthonormal basis for $K_{j}(L,v_{0}),$ for $j=k,k+1.$ Assume that the matrices $V_{k}$ and $V_{k+1}$ are governed by the following Arnoldi decomposition of order $k:$

[TABLE]

where $V_{k+1}$ is a matrix consisting $V_{k}$ in its first $k$ columns, and $\underline{H_{k+1}}$ is an upper Hessenberg matrix of order $(k+1)\times k.$ Then,

[TABLE]

*where $r=AQ_{k}U_{k,1}(:,k)+BQ_{k}U_{k,2}(:,k).$ Furthermore,

(a) if $r\in\mbox{span}\{Q_{k}\},$ then there exist vectors $x_{k}$ and $y_{k}$ such that*

[TABLE]

(b) otherwise, there exist vectors $x_{k},$ $y_{k},$ and a scalar $\beta_{k}\neq 0$ such that

[TABLE]

Observe from the Lemma-1 that the column vectors of $\begin{bmatrix}U_{k,1}\\ U_{k,2}\end{bmatrix}$ are orthonormal as $Q_{k}$ have orthonormal columns, and

[TABLE]

Also observe from the Lemma-1 that $U_{k,1}$ is an upper Hessenberg matrix if the first column of $U_{k,1}$ has non-zero elements in the rows $1$ and $2.$ But, it need not be an unreduced matrix. On the other hand, $U_{k,1}$ is an upper triangular matrix when the order of a matrix $U_{1,1}$ is $1.$ Similarly, it is easy to see that $U_{k,2}$ is an upper triangular matrix as $U_{1,2}$ is either a matrix of order $1$ or a matrix of order $2\times 1$ with a zero element at the bottom. The following lines state this discussion in the form of a lemma.

Lemma 2**.**

If the vectors $r_{-1}$ and $r_{0}$ are linearly independent, then $U_{k,1}$ is an upper Hessenberg matrix. Otherwise, it is an upper triangular matrix. $U_{k,2}$ is always an upper triangular matrix.

Next, by using the equation (4), the following lemma proves that $U_{k,1}^{\ast}Q_{k}^{\ast}(AQ_{K}U_{k,1}+BQ_{k}U_{k,2})+U_{k,2}^{\ast}U_{k,1}$ is an upper Hessenberg matrix.

Lemma 3**.**

Let $Q_{k},$ $U_{k,1},$ and $U_{k,2}$ be the same as in the Lemma-1. Then, the matrix $U_{k,1}^{\ast}Q_{k}^{\ast}(AQ_{K}U_{k,1}+BQ_{k}U_{k,2})+U_{k,2}^{\ast}U_{k,1}$ is an upper Hessenberg matrix.

Proof.

By considering the equation (4), we have

[TABLE]

Now, apply an inner product on both the sides with the matrix $\begin{bmatrix}Q_{k}U_{k,1}\\ Q_{k}U_{k,2}\end{bmatrix}.$ Furthermore, use the fact from the Lemma- 1 that column vectors of $V_{j}\equiv\begin{bmatrix}Q_{j}U_{j,1}\\ Q_{j}U_{j,2}\end{bmatrix}$ are orthonormal for $j=k,k+1$ to get the following:

[TABLE]

where $H_{k}$ is a principal submatrix of order $k$ from the top left corner of $\underline{H_{k+1}},$ and is also an upper Hessenberg matrix. Since coulmn vectors of $Q_{k}$ are orthonormal, the above equation proves the lemma. Hence the proof is over.∎

The following lemma will be helpful later. It derives a relation between the matrices $U_{k,1},$ $U_{k,2},$ and the vectors $x_{k},$ $y_{k}.$

Lemma 4**.**

Let $U_{i,1}$ and $U_{i,2}$ for $i=k,~{}k+1$ be the same as in the Lemma-1. Then, $U_{k,1}^{\ast}x_{k}=-U_{k,2}^{\ast}y_{k}.$

Proof.

Comparing the last $n$ rows of both the sides of the equation (7) gives

[TABLE]

Multiply both the sides of the above equation with $Q_{k}^{\ast}.$ Since the columns of $Q_{k}$ are orthonormal, this gives

[TABLE]

Moreover, by using the structure of $U_{k+1,2}$ from the equation (6), this implies

[TABLE]

Thus,

[TABLE]

Consequently, substituting this relation in the equation (8) gives the following:

[TABLE]

Further, by using $U_{k,2}^{\ast}U_{k,2}=I-U_{k,1}^{\ast}U_{k,1},$ the above equation becomes as follows:

[TABLE]

Now, compare the first $n$ rows of both sides of the equation (4) to observe $(AQ_{k}U_{k,1}+BQ_{k}U_{k,2})=Q_{k+1}U_{k+1,1}\underline{H_{k+1}}.$ By using this, the above equation gives

[TABLE]

As the column vectors of $Q_{k}$ are orthonormal, this implies

[TABLE]

Because of the structure of $U_{k+1,1}$ in Lemma-1, the above equation gives

[TABLE]

Finally, use $\underline{H_{k+1}}=\begin{bmatrix}H_{k}\\ h_{k+1,k}e_{k}^{\ast}\end{bmatrix}$ to get

[TABLE]

Therefore,we have proved the lemma, since $h_{k+1,k}\neq 0.$ ∎

3 Improved TOAR method

In this section, we present a few results which are the ground for proposing the new algorithm, I-TOAR. In I-ToAR, the computation of an orthonormal basis of a second-order Krylov subspace remains same as in TOAR, except I-TOAR imposes an additional condition of orthogonality on the matrix $U_{k,1}.$ In the following Lemmas, we will prove that imposing such a condition on $U_{k,1}$ prompts $U_{k,2}$ to be a diagonal matrix.

Lemma 5**.**

Let the matrix $Q_{k}U_{k,1}$ have orthogonal columns, and the columns of $Q_{k}U_{k,2}$ are $B-$ orthogonal. If $A$ is symmetric, and $U_{k,1}^{\ast}U_{k,1}+U_{k,2}^{\ast}Q_{k}^{\ast}BQ_{k}U_{k,2}$ is an identity matrix, then a matrix $H_{k}$ also symmetric.

Proof.

From the equation (7) we have

[TABLE]

Since the columns of $Q_{k}$ are orthonormal, the orthogonality of columns of $Q_{k}U_{k,1}$ force the columns of $U_{k,1}$ to be orthogonal. Now, multiply both the sides of the previous equation from the right with $U_{k,1}^{\ast}$ to get

[TABLE]

From the equation (9) note that $Q_{k}U_{k,1}=Q_{k+1}U_{k+1,2}\underline{H_{k+1}}.$ On substituting this the above equation gives:

[TABLE]

Since the columns of $Q_{k}U_{k,2}$ are $B-$ orthogonal, It satisfies the relation:

[TABLE]

Thus, the two previous equations together gives the following:

[TABLE]

Since $A$ is symmetric, $U_{k,1}^{\ast}Q_{k}^{\ast}AQ_{k}U_{k,1}$ also symmetric. This implies

[TABLE]

Now, use the fact that $U_{k,1}^{\ast}U_{k,1}+U_{k,2}^{\ast}Q_{k}^{\ast}BQ_{k}U_{k,2}$ is an Identity matrix to conclude $H_{k}$ is symmetric. Hence, the proof over. ∎

The Lemma-5 envisages parallelizing the Symmetric TOAR (STOAR) procedure in Campos & Roman (2016), by parallelly orthogonalizing columns of $U_{k,1}$ and $B-$ orthogonalizing columns of $Q_{k}U_{k,2},$ such that $U_{k,1}^{\ast}U_{k,1}+U_{k,2}^{\ast}Q_{k}^{\ast}BQ_{k}U_{k,2}$ is an Identity matrix. Though this discussion is interesting, we do not elongate it as latter distracts our attention from the paper. The following lemma generalizes the Lemma-5 for non-symmetric QEP.

Lemma 6**.**

Follow the notation of the Lemma-1. If $U_{k+1,1}$ is an orthogonal matrix then $U_{k,2}$ is a diagonal matrix.

Proof.

As $U_{k+1,1}$ is an orthogonal matrix, its principal submatrix $U_{k,1}$ also orthogonal. Further, orthonormality of the column vectors of $\begin{bmatrix}U_{k,1}\\ U_{k,2}\end{bmatrix}$ implies $U_{k,2}^{\ast}U_{k,2}$ is a diagonal matrix. Recall from the equation (10) that $U_{k,2}^{\ast}U_{k,1}=U_{k,2}^{\ast}U_{k,2}H_{k},$ provided $U_{k,2}^{\ast}y_{k}=0.$ Observe that $U_{k,2}^{\ast}y_{k}=0$ follows from the Lemma-4, and the equations (5), (6) on using the hypothesis that column vectors of $U_{k+1,1}$ are orthogonal. Now, $U_{k,2}^{\ast}U_{k,1}$ is an upper Hessenberg matrix as $H_{k}$ and $U_{k,2}^{\ast}U_{k,2}$ are upper Hessenberg and diagonal matrices, respectively. Further, as $U_{k,1}$ is an upper Hessenberg matrix from the Lemma-1, this is possible only when $U_{k,2}$ is a diagonal matrix. Therefore, the proof is complete. ∎

The above lemma is a base to propose I-TOAR method in the next section for constructing an orthonormal basis of a linear Krylov subspace associated with the given quadratic problem.

4 Implementation

This section includes two subsections. The first subsection derives relations between entries of the matrices those involved in the two successive iterations of I-TOAR. The second subsection will discuss the Improved TOAR procedure to compute the compact Arnoldi factorization for the given QEP.

4.1 Matrices in two successive iterations of I-TOAR

Assume that the initial vectors $r_{-1}~{}\mbox{and}~{}r_{0}$ are chosen randomly such that $v_{0}:=\begin{bmatrix}r_{0}\\ r_{-1}\end{bmatrix}\neq 0.$ The I-TOAR method computes a $QR$ decomposition of the $n\times 2$ matrix $[r_{-1}~{}~{}r_{0}]:$

[TABLE]

where $Q_{1}$ is an orthonormal matrix of order $n\times\alpha$ and $X$ is a matrix of size $\alpha\times 2.$ Here, $\alpha=2,$ when the vectors $r_{-1}$ and $r_{0}$ are linearly independent, otherwise $\alpha=1.$ Now, define a matrix $V_{1}$ as follows:

[TABLE]

where $\gamma=\|v_{0}\|_{2},$ $U_{1,1}=X(:,2)/\gamma,$ and $U_{1,2}=X(:,1)/\gamma$. After the initialization, we have $Q_{1},U_{1}$, and an empty $1\times 0$ Hessenberg matrix $H_{1}=[~{}~{}].$

Assume that we have the compact Arnoldi decomposition for $k=j,~{}j\geq 1.$ Similar to the TOAR method, I-TOAR enlarges this decomposition for $k=j+1$ by computing an orthogonal matrix $Q_{j+1}=[Q_{j}~{}~{}q_{j+1}]$ such that

[TABLE]

Thus, a vector $q_{j+1}$ is given by:

[TABLE]

Here, it is assumed that $\beta\neq 0.$ Note that, if $\beta=0,$ then deflation occurs and $r\in span\{Q_{j}\}.$ We state this discussion in the form of the following lemma.

Lemma 7**.**

Let for $i=j,j+1,$ column vectors of a matrix $Q_{i}$ form an orthonormal basis of a second-order Krylov subspace at the $i^{th}$ iteration of I-TOAR. Then, In case of no deflation, $Q_{j+1}=[Q_{j}~{}~{}q_{j+1}],$ where $q_{j+1}$ is given by the equation (12). Otherwise $Q_{j+1}=Q_{j}.$

The following lemma establishes relations between the matrices $U_{i},~{}i=j,j+1$ in the I-TOAR method.

Lemma 8**.**

Let $U_{i,1},~{}U_{i,2},~{}i=j,j+1$ be matrices at the iterations $j,j+1$ of I-TOAR and of the form described in the equation (6). Then, entries of an upper Hessenberg matrix $\underline{H_{j+1}}$ in I-TOAR satisfies the following relations:

[TABLE]

and

[TABLE]

where $u=U_{j,1}(:,j).$

Proof.

As column vectors of the matrix $Q_{k+1}$ are orthonormal, observe the following relation from the equations (7) and (11):

[TABLE]

The I-TOAR method uses this relation to compute the matrices $U_{j+1,1}$ and $U_{j+1,2},$ which are of the following form:

[TABLE]

From the previous two equations, it is clear that matrix $U_{j,1},$ $x_{j},$ and $\beta_{j}$ satisfy the following equation:

[TABLE]

where $h_{j}={H}_{j+1}(1:j,j)$ and $h_{j+1,j}={H}_{j+1}(j+1,j).$ Since the matrix $U_{j+1,1}$ is orthogonal in I-TOAR, the two previous equations together implies $U_{j,1}^{\ast}x_{j}=0,$ and

[TABLE]

Similarly, comparing the last row of the equation (15) gives the following using an orthonormal property of the matrix $Q_{j+1},$ and the structure of $U_{j+1,2}:$

[TABLE]

where $u\equiv U_{j,1}(:,j).$ By using the Lemma-4 and $U_{j,1}^{\ast}x_{j}=0,$ this gives

[TABLE]

Now adding the equations (17) and (19) based on the fact that column vectors of $\begin{bmatrix}U_{j,1}\\ U_{j,2}\end{bmatrix}$ are orthonormal gives the following relation:

[TABLE]

Similarly, recall the following from the equations (16) and (18):

[TABLE]

By using the fact that $(x_{j+1}~{}\beta_{j+1}~{}y_{j+1}~{}0)^{\prime}$ is a column vector of an orthonormal matrix $\begin{bmatrix}U_{j+1,1}\\ U_{j+1,2}\end{bmatrix},$ It gives the following:

[TABLE]

Hence, the equation (14) is proved by observing from the equation (16) that $h_{j+1,j}\beta_{j+1}=\beta.$ Similarly, observe that the equation (13) follows from the equations (17), (19), and (20). Therefore, the proof is complete. ∎

The Lemma-8 gives the relations to transit from $j^{th}$ to $(j+1)^{th}$ iteration in I-TOAR, provided there is no deflation at the $(j+1)^{th}$ iteration. In the following, we derive similar expressions for computing the matrix $\underline{H_{j+1}}$ in I-TOAR, in case of deflation.

Lemma 9**.**

Let $\mbox{for}~{}i=j,j+1,$ $U_{i,1}$ and $U_{i,2}$ be matrices same as in the Lemma-8, but of the form described in the equation (5). Then, entries of an upper Hessenberg matrix $\underline{H_{j+1}}$ in I-TOAR satisfies the equation (13), and also the following one:

[TABLE]

Proof.

Recall from the Lemma-1 that in the case of deflation, $Q_{j+1}=Q_{j},$ and $U_{j+1,1,}$ $U_{j+1,2}$ are of the following form:

[TABLE]

These equations together with the equation (15) gives the following relations:

[TABLE]

and

[TABLE]

These relations are similar to the equations (16) and (18) in the previous lemma. As $\begin{pmatrix}x_{j+1}\\ y_{j+1}\end{pmatrix}$ is a column vector of an orthonormal matrix $\begin{bmatrix}U_{j+1,1}\\ U_{j+1,2}\end{bmatrix},$ we have $\|x_{j+1}\|^{2}+\|y_{j+1}\|^{2}=1.$ Using this, the previous two equations gives: $h_{j+1,j}^{2}=\|s-U_{j,1}h_{j}\|^{2}+\|u-U_{j,2}h_{j}\|^{2}.$ Now, see the following relations to observe this is equal to $\|(I-U_{j,1}U_{j,1}^{\dagger})s\|^{2}+\|(I-U_{j,2}U_{j,2}^{\dagger})u\|^{2}.$

[TABLE]

The above equation follows from the equations (23) and (24) by using the identities $U_{j,1}^{\ast}x_{j}=0=U_{j,2}^{\ast}y_{j}.$ These identities follows from the fact that the matrix $U_{j+1,1}$ is orthogonal, and the Lemma-4. Therefore, we proved (22). Now, it is required to prove the equation (13).

As the matrix $\begin{bmatrix}U_{j,1}\\ U_{j,2}\end{bmatrix}$ has orthonormal columns, $U_{j,1}^{\ast}U_{j,1}+U_{j,2}^{\ast}U_{j,2}$ is an Identity matrix. Using this, the equation (25) gives $U_{j,1}^{\ast}s+U_{j,2}^{\ast}u=h_{j}.$ In turn, this equation together with (25) proves the equation (13). Hence, the proof is over. ∎

The equation (22) in the Lemma-9 shows that $h_{j+1,j}^{2}=0,$ that means, the I-TOAR algorithm break down when the vectors $s=Q_{j}^{\ast}r$ and $u$ are in the range space of matrices $U_{j,1}$ and $U_{j,2},$ respectively. Otherwise, since the column vectors of $Q_{j}$ are orthonormal, the equation (25) shows that the vector $h_{j}$ is the least squares approximation to the vectors $s$ and $u$ from the range space of $Q_{j}U_{j,1}$ and $Q_{j}U_{j,2},$ respectively.

4.2 I-TOAR implementation

This subsection uses the results of the subsection-4.1 to discuss computational details of I-TOAR. Then, it proposes the I-TOAR algorithm. Though I-TOAR follows the TOAR for constructing an orthonormal basis of the second-order Krylov subspace, I-TOAR can build an orthonormal basis for the associated linear Krylov subspace in various ways. For example, the following is the one which requires the least computation compared to all other procedures.

Procedure-1: Use the second equality relation of the equation (13) to compute the vector $h_{j}.$ It requires at most $3(j+1)$ flops. Now, use first relation in the equation (21) to compute the vector $h_{j+1,j}x_{j+1},$ and then comparing only the last element in the second relation of the equation (21) gives the vector $h_{j+1,j}y_{j+1}.$ Note that, $y_{j+1}$ is a column vector of a diagonal matrix $U_{j+1,2}$ and has only one non-zero entry. This approach takes at most $(j+1)^{2}+2(j+1)$ flops to compute $h_{j+1,j}x_{j+1},$ and $h_{j+1,j}y_{j+1}.$ After that, as already $\beta$ is known from the equation (16), it requires another $2(j+1)$ flops to compute $h_{j+1,j}.$ Then, scaling of the vectors $h_{j+1,j}x_{j+1},$ and $h_{j+1,j}y_{j+1}$ require $2(j+1)$ flops. Thus, overall computational cost is at most $(j+1)^{2}+9(j+1)$ flops for this procedure.

However, when diagonal elements of $U_{j,2}$ are smaller in magnitude, the division operation in the equation (19) introduce large floating point arithmetic errors into this procedure.

Procedure-2: Instead of the second relation as in the Procedure-1, use the first relation of the equation (13) for computing $h_{j}.$ Then, follow the Procedure-1 to compute $x_{j+1},y_{j+1},$ and $h_{j+1,j}.$ The computational cost of this procedure is equal to that of the Procedure-1.

Observe that the Procedures- 1 and 2 requires the less computation compared to $6(j+1)^{2}+3(k+1)$ flops required by TOAR to complete the same task(See, Lu et al. (2016)). However, the computed column vectors of $U_{j,1}$ are may not orthogonal, as both the procedures do not orthogonalize these vectors explicitly. The following algorithm uses the Modified Gram-Schmidt(MGS) process to do explicit orthogonalization. In the steps (k)-(n), it orthogonalizes the vector $s=Q_{j}^{\ast}r$ against all the columns of $U_{j,1}$ to compute the vector $x_{j+1}.$ Though MGS increase the computation cost, the overall computation cost in the following algorithm to perform the task that of the Procedure-1 is $5(j+1)^{2}+4(j+1)$ flops. Still, it is less compared to the TOAR method.

Similar to the TOAR algorithm in Lu et al. (2016), the Algorithm-1 also uses the MGS process to maintain orthonormality of the matrix $Q_{j}.$ The computed $Q_{j}$ may not orthonormal up to machine precision. To keep the level of orthonormality of $Q_{j}$ as close to machine precision as possible it is required to reorthogonalize a vector $r$ in the step-(g) against the columns of $Q_{j}$ by inserting the following code segments between the steps (e) and (h) of the Algorithm-1.

for i = 0, . . . , $\eta_{j}$

$\tilde{s}_{i}=q^{\ast}_{i}r~{}~{}~{}~{}$

$r=r-\tilde{s}_{i}q_{i}$

$s_{i}=\tilde{s}_{i}+s_{i}~{}$

end for

$\alpha=\|r\|_{2}~{}.$

To end this section, note that it would be possible to apply a similar reorthogonalization procedure between the steps (k) and (n) of the Algorithm-1, to keep the orthogonality of columns of the matrix $U_{j,1},$ close to machine precision. In the next section, we provide a rigorous backward error analysis of the Algorithm-1.

5 Backward Error Analysis

This section provides backward error analysis of the I-TOAR algorithm, in the presence of finite precision arithmetic. The backward error analysis for the associated linear Krylov subspace in TOAR presented in Lu et al. (2016), is also valid for I-TOAR with insignificant changes. So, this section provides backward error analysis only for the second-order Krylov subspace in I-TOAR, in terms of the matrix pair $(A,B).$ That means, this section proves that the computed basis $\hat{Q}_{j}$ in I-TOAR is an exact basis matrix of ${G}_{k}(A+\triangle A,B+\triangle B;r_{-1},r_{0})$ with small $\|[\triangle A~{}\triangle B]\|_{2}.$ Note that, for TOAR this is an open problem; See Lu et al. (2016).

Let us assume that by taking the floating point errors into account, the Compact Arnoldi decomposition computed by I-TOAR satisfies the following:

[TABLE]

where matrices with $~{}~{}\bm{\hat{}}{}~{}~{}$ on the top are the computed matrices counterpart to the matrices in the exact arithmetic, and $E$ is the error matrix. Now, we introduce two matrices $F_{mv}$ and $F$ that represent floating point error of matrix-vector product and orthogonalization process, respectively.

[TABLE]

[TABLE]

where $\bm{\hat{}}{R}_{k-1}=[\hat{r}_{1},\hat{r}_{2},\cdots,\hat{r}_{k-1}],$ and $\hat{r}_{j}=fl(A\bm{\hat{}}{Q}_{j}\bm{\hat{}}{U}_{j,1}(:,j))+B\bm{\hat{}}{Q}_{j}\bm{\hat{}}{U}_{j,2}(:,j)).$

Note that the matrices $\triangle A_{1}:=-F_{mv}\alpha(\bm{\hat{}}{Q}_{k-1}\bm{\hat{}}{U}_{k-1,1})^{\dagger}$ and $\triangle B_{1}:=-F_{mv}(1-\alpha)(\bm{\hat{}}{Q}_{k-1}\bm{\hat{}}{U}_{k-1,2})^{\dagger}$ satisfy the relation

[TABLE]

where $\alpha$ is some non-zero scalar. Then, using the equations (26) and (27) the matrices $\triangle A_{1}$ and $\triangle B_{1}$ satisfy the following relations, respectively:

[TABLE]

and

[TABLE]

Thus, from these two equations, it is easy to see that $F_{1}+F_{mv}$ is the overall error matrix in the orthogonalization process. Furthermore, on the introduction of two matrices $\triangle A_{2}:=F_{1}\alpha(\bm{\hat{}}{Q}_{k-1}\bm{\hat{}}{U}_{k-1,1})^{\dagger}$ and $\triangle B_{2}:=F_{1}(1-\alpha)(\bm{\hat{}}{Q}_{k-1}\bm{\hat{}}{U}_{k-1,2})^{\dagger},$ the previous equation becomes as follows:

[TABLE]

Now, define $[\triangle A~{}~{}\triangle B]:=[\triangle A_{1}+\triangle A_{2}~{}~{}\triangle B_{1}+\triangle B_{2}].$ Then, the equation (28) shows that the basis matrix $\bm{\hat{}}{Q}_{k}$ computed in the I-TOAR algorithm is an exact basis matrix of a second-order Krylov subspace ${G}_{k}(A+\triangle A,B+\triangle B;r_{-1},r_{0}).$ Next, to prove the I-TOAR algorithm is backward stable, it is required that the relative backward error $\frac{\|[\triangle A~{}~{}\triangle B]\|}{\|[A~{}~{}B]\|}$ is of the order of the machine precision $\varepsilon.$ To verify this, we need to derive upper bounds for $\|F_{mv}\|_{F}$ and $\|F_{1}\|_{F}.$

We adopt the following standard model for rounding errors in the floating point arithmetic. Let $\alpha$ and $\beta$ be any two real scalars. Then,

[TABLE]

where $fl(x)$ denotes the computed quantity, and $\varepsilon$ denotes the machine precision. We will use the following lemma also in the backward error analysis.

Lemma 10**.**

*(a). For $x,y\in\mathcal{R}^{n},~{}fl(x+y)~{}=~{}x+y+f,$ where $\|f\|_{2}\leq(\|x\|_{2}+\|y\|_{2}).$

(b). For $X\in\mathcal{R}^{n\times k}$ and $y\in\mathcal{R}^{k}$, $fl(Xy)~{}=~{}Xy+w,$ where $\|w\|_{2}\leq k\|X\|_{F}\|y\|_{2}\varepsilon+\mathcal{O}(\varepsilon^{2}).$

(c). For $X\in\mathcal{R}^{n\times k},$ $y\in\mathcal{R}^{k},~{}b\in\mathcal{R}^{n},$ and $\beta\in\mathcal{R},$ $\hat{c}\equiv fl((b-Xy)/\beta)$ satisfies*

[TABLE]

For the Proof of the Lemma-10, see the Lemma-4.1 in Lu et al. (2016). This lemma holds true for I-TOAR as well. Next, the following lemma derives an upper bound for $\|F_{mv}\|_{F},$ where $F_{mv}$ is a matrix defined as in the equation (26).

Lemma 11**.**

Let $\bm{\hat{}}{Q}_{k-1}$ and $\bm{\hat{}}{U}_{k-1}$ be orthonormal matrices computed by the I-TOAR procedure. Then,

[TABLE]

Proof.

From the definition of $F_{mv}$ in the equation (26), we have

[TABLE]

where \bm{\hat{}}{r}_{j}\equiv fl\Big{(}[A~{}~{}B]\begin{bmatrix}\bm{\hat{}}{Q}_{j}\bm{\hat{}}{U}_{j,1}(:,j)\\ \hat{Q}_{j}\bm{\hat{}}{U}_{j,2}(:,j)\end{bmatrix}\Big{)}. Here, we used the fact that $\bm{\hat{}}{U}_{k-1,1},$ and $\bm{\hat{}}{U}_{k-1,2}$ are upper Hessenberg and diagonal matrices, respectively and $j\leq k-1$. Now, by the repeated application of the Lemma-10, we have

[TABLE]

[TABLE]

where $w_{j}^{(1)},$ $w_{j}^{(2)},$ and $w_{j}$ are the floating point error vectors satisfying the following relations, respectively.

[TABLE]

[TABLE]

and

[TABLE]

Now, combine all the three previous inequalities, and use the facts $\bm{\hat{}}{U}_{j}=\begin{bmatrix}\bm{\hat{}}{U}_{j,1}\\ \bm{\hat{}}{U}_{j,2}\end{bmatrix},$ and $\|\bm{\hat{}}{U}_{j,i}(:,j)\|_{2}\leq\|\bm{\hat{}}{U}_{j}(:,j)\|_{2},~{}\mbox{for}~{}i=1,2.$ It gives

[TABLE]

Then, using $\|\bm{\hat{}}{Q}_{j}\|_{F}\leq\|\bm{\hat{}}{Q}_{k-1}\|_{F}$ for $j\leq k-1,$ the above inequality gives the following:

[TABLE]

Observe that in matrix terms this equation can be written as follows:

[TABLE]

Further, using the inequalities $\|\bm{\hat{}}{Q}_{j}\|_{F}\leq\|\bm{\hat{}}{Q}_{k-1}\|_{F}\leq\sqrt{k}\|\bm{\hat{}}{Q}_{k-1}\|_{2},$ and $\|\bm{\hat{}}{U}_{k-1}\|_{F}\leq\sqrt{k}\|\bm{\hat{}}{U}_{k-1}\|_{2}$ this gives

[TABLE]

Since the I-TOAR Algorithm uses the MGS process to generate orthonormal matrices $\bm{\hat{}}{Q}_{k-1}$ and $\bm{\hat{}}{U}_{k-1}$, we have $\|\bm{\hat{}}{Q}_{k-1}\|_{2}=\|\bm{\hat{}}{U}_{k-1}\|_{2}=1+\mathcal{O}(\varepsilon).$ Therefore, neglecting the term $2k^{2}\|\bm{\hat{}}{Q}_{k-1}\|_{2}\|\bm{\hat{}}{U}_{k-1}\|_{2}\varepsilon$ in the above equation completes the proof of the lemma. ∎

Next, to derive an upper bound for $\|F_{1}\|_{F},$ the following lemma is required.

Lemma 12**.**

Let $f$ be an error vector resulting from the computation in the steps (k)-(n) of the Algorithm-1. Then

[TABLE]

Proof.

The elements of computed vector $\alpha_{1}$ in the step-(l) of the Algorithm-1 satisfy the relation:

[TABLE]

where $f_{i}$ is the error resulting from the computation of an inner product in the numerator and the division. By using the fact that $\bm{\hat{}}{U}_{j,1}^{\ast}\bm{\hat{}}{U}_{j,1}$ is a diagonal matrix, the above element-wise computation can be written in the vector form as follows:

[TABLE]

where an error vector $f$ satisfies the following relation:

[TABLE]

Therefore, using $\|\bm{\hat{}}{U}_{j,1}\|_{F}=\|\bm{\hat{}}{U}_{j,1}^{\ast}\|_{F},$ this gives the inequality

[TABLE]

Let $\sigma_{max}$ and $\sigma_{min}$ be the largest and the smallest singular values of $\bm{\hat{}}{U}_{j,1}^{\ast},$ respectively. Then, we have $\|\bm{\hat{}}{U}_{j,1}\|_{F}^{2}\leq(j+1)\sigma_{max}^{2}$ and $\|(\bm{\hat{}}{U}_{j,1}^{\ast}\bm{\hat{}}{U}_{j,1})^{-1}\|_{F}\leq\frac{j+1}{\sigma_{min}^{2}}.$ Substitute these inequalities in the previous equation and use $\mathcal{K}_{2}(\bm{\hat{}}{U}_{j,1})^{2}=\frac{\sigma_{max}^{2}}{\sigma_{min}^{2}}$ to complete the proof. ∎

In the next lemma, we use the Lemma-12 to bound the norm of error vector in the second orthogonalization process in the steps (k) - (n) of the Algorithm-1. Moreover, we are assuming that $\mathcal{K}_{2}(\bm{\hat{}}{U}_{j,1})$ is moderately small.

Lemma 13**.**

Let $g_{1,j}$ be a overall floating point error vector resulting from the second level orthogonalization process in the steps (k)-(n) of the Algorithm-1. Then

[TABLE]

Proof.

On applying Lemma-10(c) to the second orthogonalization process in the steps (k) - (n), the computed $(j+1)$th column of $U_{k+1,1}$ satisfies,

[TABLE]

where $g_{1,j}$ is a floating point error vector satisfying the inequality,

[TABLE]

where $f$ is the error vector same as in the previous lemma. In the above equation, we used the fact $h=\gamma_{1}=(\bm{\hat{}}{U}_{j,1}^{\ast}\bm{\hat{}}{U}_{j,1})^{-1}\bm{\hat{}}{U}_{j,1}^{\ast}\bm{\hat{}}{s}.$ Therefore, the proof is complete by using the fact from the Lemma-12 that $\|\bm{\hat{}}{U}_{j,1}\|_{F}\|f\|_{2}\varepsilon$ in the above inequality is of $\mathcal{O}(\varepsilon^{2}).$ ∎

In the following, we use the Lemma-13 to derive an upper bound for $\|F_{1}\|_{F}^{2}.$

Lemma 14**.**

Let $\bm{\hat{}}{Q}_{k},$ $\bm{\hat{}}{U}_{k}$ and $\underline{\hat{H}_{k}}$ be computed in the k-step of the I-TOAR procedure. Then

[TABLE]

where $\varphi=(k+1)(2k+1).$

Proof.

Note that, in the $j$th iteration of I-TOAR, the computed quantity at step-(d) is

[TABLE]

Now, apply the Lemma-10(c) to the orthogonalization and normalization processes in the steps (e) - (h) and (x) of the Algorithm-1. Then, the computed column $\bm{\hat{}}{q}_{j+1}$ of $Q_{k}$ and $\bm{\hat{}}{s},$ $\bm{\hat{}}{\beta}$ computed in the steps (m) and (i) satisfy the following equation

[TABLE]

where $\tilde{f}_{j}$ is the error vector, which satisfies the following:

[TABLE]

Now, from the equation (27), consider the $j$th column of $F_{1},$

[TABLE]

where, for the second equality, we exploited the upper Hessenberg structure of $\bm{\hat{}}{H}_{k}.$ Moreover, from the equation (6), we have

[TABLE]

Further, on left multiplying the equation (30) with $\bm{\hat{}}{Q}_{j},$ we have

[TABLE]

Hence, the substitution of the equations (32) and (35) in the equation (34) will give

[TABLE]

Thus,

[TABLE]

In the second, third and fourth inequalities, we used the equations (33), (30) and (29), respectively, whereas the following were used to obtain the last inequation:

[TABLE]

Now, from the equation (5), $\|F_{1}\|_{F}$ is given by

[TABLE]

[TABLE]

Therefore, the proof will be complete by converting the Frobenius norm to 2-norm. ∎

Remark-1: When the complete reorthogonalization has applied, the above result holds true with little change in the coefficients. In this case, following a similar procedure as in the remark-1 in Lu et al. (2016), it is easy to see that the norm of the error vector $\tilde{f}_{j}$ satisfies the following relation:

[TABLE]

where $c$ is a small constant.

We now present the main theorem for an upper bound of relative backward error in the I-TOAR procedure. Moreover, we are using the following assumption in the proof:

[TABLE]

where $\sigma_{min}(X)$ denotes the smallest singular value of a matrix $X.$

Theorem 1**.**

Let $\bm{\hat{}}{Q}_{k},~{}\bm{\hat{}}{U}_{k,1}$ and $\bm{\hat{}}{U}_{k,2}$ be matrices of full column rank. Let

[TABLE]

where for any matrix $X,$ $\mathcal{K}_{2}(X)$ denotes its 2-norm condition number. If $(k+1)(2k+1)\mathcal{K}^{4}(\zeta_{1})/\zeta_{2}\varepsilon<1,$ then

[TABLE]

where $\varphi_{2}=(k+1)(2k+1)\zeta_{1}/\zeta_{2}.$

Proof.

We have

[TABLE]

Using the definition of $\triangle A_{1}$ and $\triangle B_{1}$, we have the following for $\alpha=1/2.$

[TABLE]

Similarly, using the definition of $\triangle A_{2}$ and $\triangle B_{2},$ we have

[TABLE]

Further, by using the equation (37) we have the following inequality:

[TABLE]

[TABLE]

In addition, with the bound for $\|F_{mv}\|_{F}$ in Lemma-11, this gives the following upper bound for $\|[\triangle A_{1}~{}~{}\triangle B_{1}]\|_{F}.$

[TABLE]

Now, recall from the equation (28) that

[TABLE]

and repeatedly apply the inequality, $\|XY\|_{F}\leq\|X\|_{2}\|Y\|_{F}$ to obtain the following:

[TABLE]

Then, use the following inequalities in the above equation,

[TABLE]

to get the following:

[TABLE]

By using this equation and the result in the Lemma-14, we have

[TABLE]

where $\varphi_{2}=(2k+1)(k+1)\zeta_{1}/\zeta_{2}.$ The two inequalities in the equation (37) were used for the first and second inequalities, respectively . Now, by combining the equations (5) and (5), we have

[TABLE]

As we assumed $(\varphi_{2}/2)\mathcal{K}^{4}\varepsilon<1,$ the theorem is proven by omitting the denominator, since it can be covered by the term $\mathcal{O}(\varepsilon^{2}).$ ∎

Remark-2: Similar to the Theorem-2.5 in Lu et al. (2016), the previous lemmas were assumed that the matrices A and B are known explicitly so that the standard error bound for the matrix-vector multiplication applicable. The stability analysis of the Arnoldi method in (Stewart, 2001, Theorem-2.5) also used the same assumption. Otherwise, an error bound for matrix-vector multiplication depends on the specific formulation of $A$ and $B.$

6 Numerical examples

In this section, we apply the I-TOAR procedure to the application of Model order reduction of second-order dynamical systems. A continuous time invariant dynamical system is of the following form:

[TABLE]

where $M,~{}D,~{}K\in\mathcal{R}^{n\times n}$ are mass, damping, and stiffness matrices, respectively. $F\in\mathcal{R}^{n\times m},$ $C_{p},~{}C_{v}\in\mathcal{R}^{q\times n}$ are constant matrices. In this paper, we are assuming $C_{v}=0$ and $m=q=1.$ Thus, we considering the following single input-single output dynamical system of the form:

[TABLE]

where $x(t)$ is the state vector, $u(t)$ is the input vector and $y(t)$ is the output vector. Here, $f$ is a input distribution array and $c$ is the outer measurement array. Moreover, for the convenience, we are assuming that $x(0)=\dot{x}(0)=0.$

Applying the Laplace transform on both sides of the previous equation will give

[TABLE]

where $X(s),Y(s),$ and $U(s)$ are Laplace transforms of $x(t),y(t),$ and $u(t),$ respectively. Thus, we have

[TABLE]

$h(s):=c(s^{2}M+sD+K)^{-1}f$ is called as the Transfer function. Using the definition of $h(s),$ the previous equation can be written as

[TABLE]

The Model order reduction method produces a lower order dynamical system that closely resembles the characteristics of the original system. Though model order reduction is possible in many ways, in this paper, we are using the Galerkin projection based reduction method. This method defines a projection operator using an orthonormal basis of a subspace generated by the mass, damping and stiffness matrices. Then, it projects the system (6) onto a subspace of smaller dimension. The following is the resulting reduced model of a system in the equation (6):

[TABLE]

The above system will be solved for $X_{k}(s).$ Note that $X_{k}(s)$ has a lesser number of elements compared to $X(s)$ in (6). Now, the approximation to $Y(s)$ is given by

[TABLE]

Here, $h_{k}(s)$ is an approximation to the Transfer function $h(s)$ and is given by

[TABLE]

where $c_{k}=cV,~{}M_{k}=V^{\ast}MV,~{}D_{k}=V^{\ast}DV$ and $K_{k}=V^{\ast}KV.$

The main objective of model order reduction techniques is to compute $h_{k}(s),$ as an accurate approximation of $h(s)$ over a wide range of frequency intervals around a prescribed shift $s_{0}.$ As in Lu et al. (2016), to meet this objective, we rewrite the transfer function $h(s)$ by including the shift $s_{0}$ as follows:

[TABLE]

Using I-TOAR, we compute an orthonormal basis matrix $Q_{k}\in\mathcal{R}^{n\times\eta_{k}}$ of the second order Krylov subspace

[TABLE]

where $\tilde{D}=2s_{0}M+D$ and $\tilde{K}=s_{0}^{2}M+s_{0}D+K.$ Then, we used $Q_{k}$ in place of $V$ in the equation (43) for computing $h_{k}(s),$ an approximation to $h(s).$

Our numerical experiments compare the accuracy of reduced dynamical systems defined using the orthonormal basis matrices in the TOAR and I-TOAR methods. In both the TOAR and I-TOAR methods, we apply the reorthogonalization to ensure that the computed basis is orthonormal up to the machine precision. For I-TOAR, we used the complete reorthogonalization whereas, for the TOAR method we used the same setup as in Lu et al. (2016). All algorithms are implemented using MATLAB and were run on a machine Intel(R) Core(TM)i7-4770 CPU@3.40GHz with 8GB RAM. For the convenience, we use the same examples as in Lu et al. (2016). The author would acknowledge D.LU, an author of Lu et al. (2016) for providing the data of these examples.

Example 1**.**

This example is a finite element model of a shaft on bearing supports with a damper in MSC/NASTRAN. It is a second-order system and of dimension 400. The mass, damping matrices are symmetric, and the stiffness matrix is symmetric positive definite. We use the expansion point $s_{0}=150\times 2\pi$ to approximate the Transfer function $h(s)$ over the frequency interval $[0,3000].$

The left plot of Figure-1 shows the magnitudes of the transfer function $h(s)$ of the full-order system, and the transfer functions $h_{k}(s)$ of the reduced systems generated by the I-TOAR and TOAR procedures for $k=40.$ The relative errors of the transfer functions in the I-TOAR and TOAR procedures for $k=10,20$ and $40$ are shown in the middle and right plots of the Figure-1 respectively. As we can see that, the transfer functions $h_{k}(s)$ in the I-TOAR and TOAR methods produces almost the same accuracy in the frequency interval $[0,2000],$ for $k=10,20,$ and $40.$ In the frequency interval $[2000,3000],$ the transfer function $h_{k}(s)$ by I-TOAR is a more accurate approximation than the one produced by the TOAR method. Like the TOAR method, in the I-TOAR method also, the approximation accuracy of $h_{k}(s)$ is improved, when increasing $k$ from $10$ to $40.$

Example 2**.**

This example is the butterfly gyroscope problem from the Oberwolfach collection. The full dynamical system is of the order $n=17361.$ The mass and stiffness matrices $M$ and $K$ are symmetric. The damping matrix is of the form $D=\alpha M+\beta K.$ This second-order system have 1 input vector and 12 output vectors. Following the experiments in Lu et al. (2016), we considered the first output vector as the output vector ’c.’ The damping parameters $\alpha$ and $\beta$ are chosen same as in Lu et al. (2016); Li et al. (2012), $\alpha=0$ and $\beta=10^{-7}.$ The expansion point $s_{0}$ also same as in Lu et al. (2016), $s_{0}=1.05\times 10^{5}.$

The magnitudes of the transfer functions shown in the left plot of the Figure-2. The relative errors of transfer functions in I-TOAR and TOAR show in the middle and right plots of the Figure-2. From the figure, it is easy to observe the advantage of I-TOAR over TOAR for the frequency range of $10^{5}-10^{6}Hz$ for $k=90$ and $200.$ For $k=20$ and $30$, both I-TOAR and TOAR produced nearly the same accuracy. We have observed the stagnation in both TOAR and I-TOAR from $k=200$ onwards. From the middle and right plots of the Figure-2, it is clear that the transfer function in I-TOAR is more accurate than the transfer function in the TOAR procedure. Moreover, for this example, we found that the quantity $\zeta_{1}$ of Theorem-1 is of order $10^{11}$ and $\zeta_{2}=1.000000000001526e+000.$ Further, found that condition number of the computed matrices $\bm{\hat{}}{Q}_{k}$ and $\bm{\hat{}}{U}_{k}$ are equal to $1+\mathcal{O}(\varepsilon),$ and these quantities satisfy the condition $(\varphi_{2}/2)\mathcal{K}^{4}\varepsilon<1.$

7 Conclusion

In this paper, we have proposed a new TOAR procedure. It imposes an extra condition on the orthogonality of the matrices in the second-level orthogonalization of TOAR. Imposing the new condition gives orthonormal basis of an associated linear Krylov subspace without any extra computation. A rigorous stability analysis has done on the proposed method. The backward analysis is in terms of the matrix $[A~{}~{}B]$ of the quadratic problem. It has shown that in the proposed method the second-order Krylov subspace of $[A~{}~{}B]$ embedded in that of $[A+\triangle A~{}~{}B+\triangle B]$ for sufficiently small $\|[\triangle A~{}~{}\triangle B]\|.$ This problem was left open for TOAR in Lu et al. (2016). Numerical experiments have shown that the basis matrices in I-TOAR are as accurate as ones in TOAR in the application of dimension reduction in second-order dynamical systems. The method proposed in this paper may help us to improve the methods for solving polynomial eigenvalue problems.

Acknowledgements

The author would like to acknowledge D.LU, author of Lu et al. (2016) for providing the data used in the numerical examples. Further, this work was supported by the National Board of Higher Mathematics, India under Grant number 2/40(3)/2016/R&D-II/9602.