Block Krylov subspace methods for large-scale matrix computations in control

Peer review under responsibility of Taibah University.

Abstract

In this paper we show how to improve the approximate solution of the large Sylvester equation obtained by an arbitrary method. Such problems appear in many areas of control theory such as the computation of Hankel singular values, model reduction algorithms and others. Moreover, we propose a new method based on refinement process and weighted block Arnoldi algorithm for solving large Sylvester matrix equation. The numerical tests report the effectiveness of these methods.

MSC:

• 65F30
• 65F50

Keywords:

• Matrix equations
• Sylvester
• Arnoldi
• Refinement

1 Introduction

An important class of linear matrix equations is the Sylvester equation

(1.1) $\mathrm{XA}+\mathrm{BX}=C.$(1.1) Since the fundamental work of Sylvester on the stability of the motion, these matrix equations have been widely used in stability theory of differential equations and play an important role in control and communications theory [1–5]. For small problems, direct methods for solving the Sylvester equation are attractive. The standard methods are based on the Hessenberg–Schur decomposition to transform the original equation into a form that can be easily solved. Iterative methods for solving large Sylvester matrix equations have been developed during the past years [6–8]. A class of classical methods known as the Krylov subspace methods, that include the block Arnoldi and weighted block Arnoldi, etc. have been found to be suitable for sparse matrix computations. In this paper, we extend the idea to propose a new projection method for solving (1.1)(1.1) $\mathrm{XA}+\mathrm{BX}=C.$(1.1) based on weighted block Krylov subspace method. The paper is organized as follows. In Section 2, we introduce the main minimal residual algorithm for iteratively solving the Sylvester matrix equation (1.1)(1.1) $\mathrm{XA}+\mathrm{BX}=C.$(1.1) . Several numerical examples presented in Section 3. Finally, the conclusion is given in the last section.

2 Block refinement Arnoldi method

In this section we propose to show that the obtained approximate solution of the Sylvester equation by any method can be improved, in other words the accuracy can be increased. If this idea is applicable then we have found an iterative method for solving of the Sylvester equation. Therefore let the basis $V_m=[v_1,\dots ,v_m]$ and $W_m=[w_1,\dots ,w_m]$ constructed by the block Arnoldi process, thus we have

$V_m^T{V}_m=I_m{W}_m^T{W}_m=I_m$The square block Hessenberg matrices H_m and ${\hat{H}}_m$ (m = r * l where r and l are the dimensions of blocks) whose nonzero entries are the scalars h_ij and ${\hat{h}}_{\mathrm{ij}}$, constructed by the block Arnoldi process can be expressed as

$H_m=V_m^T{A}^T{V}_m{\hat{H}}_m=W_m^T{\mathrm{BW}}_m$Let X⁽⁰⁾ be an initial approximate solution of the Sylvester equation XA + BX = C where A, B, C, X⁽⁰⁾, ∈ R^n×n.

Also introduce the residual matrix

$R^{(0)}=C-(X^{(0)}A+{\mathrm{BX}}^{(0)})$And let Y_m ∈ R^m×m be the solution of the Sylvester equation:

(2.1) $Y_m{H}_m^T+{\hat{H}}_m{Y}_m=W_m^T{R}^{(0)}{V}_m$(2.1) If set

(2.2) $X^{(1)}=X^{(0)}+W_m{Y}_m{V}_m^T$(2.2) then the corresponding residual R⁽¹⁾ = C − (X⁽¹⁾A + BX⁽¹⁾) satisfies:

$\begin{array}{ccc}{R}^{(1)}\hfill & =\hfill & C-((X^{(0)}+W_m{Y}_m{V}_m^T)A+B(X^{(0)}+W_m{Y}_m{V}_m^T))\hfill \\ \hfill & =\hfill & R^{(0)}-W_m{Y}_m{V}_m^{T}A-{\mathrm{BW}}_m{Y}_m{V}_m^T\hfill \\ \hfill & =\hfill & R^{(0)}-W_m{Y}_m{H}_m^T{V}_m^T-W_m{\hat{H}}_m{Y}_m{V}_m^T\hfill \\ \hfill & =\hfill & R^{(0)}-W_m(Y_m{H}_m^T+{\hat{H}}_m{Y}_m)V_m^T\hfill \end{array}$Since Y_m is the solution of (2.1)(2.1) $Y_m{H}_m^T+{\hat{H}}_m{Y}_m=W_m^T{R}^{(0)}{V}_m$(2.1) we have:

$R^{(1)}=R^{(0)}-W_m{W}_m^T{R}^{(0)}{V}_m{V}_m^T=0$According to the above results we can develop an iterative method for the solving of the Sylvester equation when the matrices A, B and C are large and sparse. For doing this idea if we choose m < n, then instead of solving XA + BX = C we can solve (2.1)(2.1) $Y_m{H}_m^T+{\hat{H}}_m{Y}_m=W_m^T{R}^{(0)}{V}_m$(2.1) . In other words in this method, first we transform the initial Sylvester equation to another Sylvester equation with less dimensions, then in each iteration step solve this matrix equation and extend the obtained solution to the solution of initial equation by (2.2)(2.2) $X^{(1)}=X^{(0)}+W_m{Y}_m{V}_m^T$(2.2) . The algorithm is as follows:

Algorithm 1

Block refinement Arnoldi method (BRA)

Choose an initial solution X⁽⁰⁾, and a tolerance ϵ and select two numbers r and l for dimensions of block and set m = r * l(m < n).

Table

For k = 0… until Convergence.

R^(k) = C − (X^(k)A + BX^(k)).

Construct the orthonormal basis Vm and Wm ∈ R^n×m by the block Arnoldi process

$H_m=V_m^T{A}^T{V}_m{\hat{H}}_m=W_m^T{\mathrm{BW}}_m.$

Solve the reduced Sylvester equation

$Y_m{H}_m^T+{\hat{H}}_m{Y}_m=W_m^T{R}^{(k)}{V}_m.$

$X^{(k+1)}=X^{(k)}+W_m{Y}_m{V}_m^T$.

if $\frac{\parallel X^{(k+1)}-X^{(k)}\parallel }{\parallel X^{(k)}\parallel }\le ϵ$ Stop

else End.

3 Numerical examples

In this section, we present some numerical examples to illustrate the effectiveness of algorithms described in this paper for large and sparse Sylvester equations. All numerical tests are performed in MATLAB software on a PC with 2.20 GHz with main memory 2 GB.

Example 3.1

Consider the Sylvester equation XA + BX = C with n = 100. We apply the Iterative block refinement Arnoldi method for solving this matrix equation and take ϵ = 10⁻⁶. In Table 1, we report the results for different values of m.

$A={\left(\begin{array}{cccccccccc}\hfill 10\hfill & \hfill 1.2\hfill & \hfill .42\hfill & \hfill .8\hfill & \hfill 2.3\hfill & \hfill .8\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 1.8\hfill & \hfill 10\hfill & \hfill 1.2\hfill & \hfill .42\hfill & \hfill .8\hfill & \hfill 2.3\hfill & \hfill .8\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill 1.6\hfill & \hfill 1.8\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 0\hfill & \hfill ⋮\hfill \\ \hfill 1.64\hfill & \hfill 1.6\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .8\hfill & \hfill 0\hfill \\ \hfill 1.3\hfill & \hfill 1.64\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 2.3\hfill & \hfill .8\hfill \\ \hfill 1.61\hfill & \hfill 1.3\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .8\hfill & \hfill 2.3\hfill \\ \hfill 0\hfill & \hfill 1.61\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .42\hfill & \hfill .8\hfill \\ \hfill ⋮\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 1.2\hfill & \hfill .42\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill 0\hfill & \hfill 1.61\hfill & \hfill 1.3\hfill & \hfill 1.64\hfill & \hfill 1.6\hfill & \hfill 1.8\hfill & \hfill 10\hfill & \hfill 1.2\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 1.61\hfill & \hfill 1.3\hfill & \hfill 1.64\hfill & \hfill 1.6\hfill & \hfill 1.8\hfill & \hfill 10\hfill \end{array}\right)}_{n\times n}$$B={\left(\begin{array}{cccccccccc}\hfill 10\hfill & \hfill 2.1\hfill & \hfill .38\hfill & \hfill .7\hfill & \hfill 1.5\hfill & \hfill .4\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 1.21\hfill & \hfill 10\hfill & \hfill 2.1\hfill & \hfill .38\hfill & \hfill .7\hfill & \hfill 1.5\hfill & \hfill .4\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill 1.9\hfill & \hfill 1.21\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 0\hfill & \hfill ⋮\hfill \\ \hfill .64\hfill & \hfill 1.9\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .4\hfill & \hfill 0\hfill \\ \hfill 1.9\hfill & \hfill .64\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 1.5\hfill & \hfill .4\hfill \\ \hfill .87\hfill & \hfill 1.9\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .7\hfill & \hfill 1.5\hfill \\ \hfill 0\hfill & \hfill .87\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .38\hfill & \hfill .7\hfill \\ \hfill ⋮\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 2.1\hfill & \hfill .38\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill 0\hfill & \hfill .87\hfill & \hfill 1.9\hfill & \hfill .64\hfill & \hfill 1.9\hfill & \hfill 1.21\hfill & \hfill 10\hfill & \hfill 2.1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill .87\hfill & \hfill 1.9\hfill & \hfill .64\hfill & \hfill 1.9\hfill & \hfill 1.21\hfill & \hfill 10\hfill \end{array}\right)}_{n\times n}$$C={\left(\begin{array}{cccccccccc}\hfill .1\hfill & \hfill 2.21\hfill & \hfill 1.4\hfill & \hfill 1.5\hfill & \hfill .13\hfill & \hfill 2.62\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 1.3\hfill & \hfill .1\hfill & \hfill 2.21\hfill & \hfill 1.4\hfill & \hfill 1.5\hfill & \hfill .13\hfill & \hfill 2.62\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill 2.6\hfill & \hfill 1.3\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 0\hfill & \hfill ⋮\hfill \\ \hfill 1.7\hfill & \hfill 2.6\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 2.62\hfill & \hfill 0\hfill \\ \hfill 2.3\hfill & \hfill 1.7\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill .13\hfill & \hfill 2.62\hfill \\ \hfill 2.6\hfill & \hfill 2.3\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 1.5\hfill & \hfill .13\hfill \\ \hfill 0\hfill & \hfill 2.6\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 1.4\hfill & \hfill 1.5\hfill \\ \hfill ⋮\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 2.21\hfill & \hfill 1.4\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill 0\hfill & \hfill 2.6\hfill & \hfill 2.3\hfill & \hfill 1.7\hfill & \hfill 2.6\hfill & \hfill 1.3\hfill & \hfill .1\hfill & \hfill 2.21\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 2.6\hfill & \hfill 2.3\hfill & \hfill 1.7\hfill & \hfill 2.6\hfill & \hfill 1.3\hfill & \hfill .1\hfill \end{array}\right)}_{n\times n}$

Table 1 Implementation of iterative BRA method to solve the Sylvester equation with different values of m.

m	r	l	$\parallel V_m^T{A}^T{V}_m-H_m\parallel$	$\parallel W_m^T{\mathrm{BW}}_m-{\hat{H}}_m\parallel$	Iteration	Time
4	2	2	8.63E−014	4.05E−014	278	5.83
8	2	4	1.57E−013	5.49E−014	156	4.68
10	2	5	2.81E−013	1.98E−014	95	3.56
20	2	10	2.84E−013	1.77E−013	58	2.68
30	2	15	3.05E−013	4.46E−014	36	1.74
40	2	20	9.77E−013	8.36E−014	21	0.7811
50	2	25	3.17E−012	4.68E−013	2	0.0918

Example 3.2

Now consider A and B and C are the same matrices that used in Example 3.1. We apply two Iterative methods and Hessenberg–Schur method to solve the Sylvester equation when the dimension of the matrices are large. Results are shown in Table 2.

Table 2 Implementation of new Iterative methods and Hessenberg–Schur method for solving the Sylvester equation.

n	Hessenberg–Schur method		Weighted Krylov method		BRA method		cond(B)
	Error	Time	Error	Time	Error	Time
200	1.16E−010	0.6881	2.50E−012	0.4753	2.38E−014	0.2612	8.53E+003
400	3.89E−007	6.312	4.22E−008	4.642	6.15E−014	3.134	3.17E+004
600	0.0011	28.89	0.0033	65.76	6.95E−014	21.39	7.63E+005
800	8.931	85.74	13.01	174.32	8.37E−014	58.26	2.13E+007
1000	27.35	201.14	48.19	322.11	1.84E−013	121.53	1.50E+008

4 Comments and conclusion

In this paper, we introduced a new method for computing low-rank approximate solutions to large Sylvester matrix equations. Moreover Refinement process presented in Section 2 has the capability of improving the results obtained by an arbitrary method. For example in this paper we apply the refinement process with Hessenberg–Schur method. The experiments illustrate the advantages of the method for large sparse matrices.

Notes

Peer review under responsibility of Taibah University.