1 Introduction
Consider a linear timeinvariant multiinput, multioutput linear time independent (LTI) dynamical system described by the statespace equations
(1) 
where
denotes the state vector,
and are the input and the output signal vectors, respectively. The matrix is assumed to be large and sparse and . The transfer function associated to the system in (1) is given as(2) 
The goal of our model reduction approach consists in defining two orthogonal matrices and (with ) to produce a much smaller order system with the statespace form
(3) 
and its transfer function is defined by
(4) 
where , and ,
such that the reduced system will have an output as close as possible to the one of the original system to any
given input , which means that for some chosen norm, is small.
Various model reduction techniques, such as Padé approximation
[17, 27], balanced truncation [28], optimal Hankel norm [16] and Krylov subspace methods, [8, 9, 14, 21] have been used for large multiinput multioutput (MIMO) dynamical systems, see [25, 16, 20]. Balanced Truncation Model Reduction (BTMR) method is a very popular method; [1, 15]; the method preserves the stability and provides a bound for the approximation error. In the case of small to medium systems, (BTMR) can be implemented efficiently. However, for largescale settings, the method is quite expensive to implement, because it requires the computation of two Lyapunov equations, and results in a computational complexity of and a storage requirement of , see [1, 5, 18]. In this paper, we project the system (1) onto the following block tangential Krylov subspaces defined as,in order to obtain a small scale dynamical system. The and are the right and left interpolation points, the and are the right and left blocks tangent directions with with . Later, we will show how to choose these tangent interpolation points and directions.
The paper is organized as follows: In section 2 we give some definition used later and we introduce the tangential interpolation. In Section 3, we present the tangential block Lanczostype method and the corresponding algorithm. Section 4 is devoted to the selection of the interpolation points and the tangential directions that are used in the construction of block tangential Krylov subspaces, and we present briefly the adaptive tangential block Lanczostype algorithm. Section 5 we treat the model reduction of secondorder systems. The last section is devoted to numerical tests and comparisons with some well known model order reduction methods.
2 Moments and interpolation
We first give the following definition. Given the system , its associated transfer function can be decomposed through a Laurent series expansion around a given (shift point), as follows
(5) 
where is called the
th moments at
associated to the system and defined as follows(6) 
In the case where the moments are called Markov parameters and are given by
Problem: Given a fullorder model (1) and assume that the following parameters are given:

Left interpolation points & block tangent directions .

Right interpolation points & block tangent directions .
The main problem is to find a reducedorder model (3) such that the associated transfer function, in (4) is a tangential interpolant to in (2), i.e.
(7) 
The interpolation points and tangent directions are selected to realize the model reduction goals described later.
3 The block tangential Lanczostype method
Let the original transfer function be expressed as where and are such that,
(8) 
Given a system of matrices and where , the approximate solution and of and are computed such that
(9) 
and
(10) 
(11) 
where , , and are the th columns of , , and , respectively. If we set and , then from (9), (10) and (11) , we obtain
which gives the following approximate transfer function
where , and . The matrices and are biorthonormal, where the . Notice that the residuals can be expressed as
(12) 
(13) 
3.1 Tangential block Lanczostype algorithm
This algorithm consists in constructing two biorthonormal bases, spanned by the columns of and , of the following block tangential Krylov subspaces
(14) 
and
(15) 
where and are the right and left interpolation points respectively and , are the right and left tangent directions with . We present next the Block Tangential Lanczos (BTL) algorithm that allows us to construct such bases. It is summarized in the following steps.
Let and . Then, we should have the biorthogonality conditions for :
(16) 
Here we suppose that we already have the set of interpolation points , and the tangential matrix directions and . The upper block upper Hessenberg matrices and are obtained from the BTL algorithm, with
The matrices and constructed in Step 3 of Algorithm 1 are of size and 0
is the zero matrix of size
. We define also the following matrices,where and . With all those notations, we have the following theorem. Let and be the biorthonormal matrices of constructed by Algorithm 1. Then we have the following relations
(17) 
and
(18) 
Let and be the matrices,
then we have
(19) 
where and . The matrices and are block upper triangular matrices of sizes and are assumed to be nonsingular.
From Algorithm 1, we have
(20) 
multiplying (20) on the left by and rearranging terms, we get
which gives
that written as
(21) 
where 0 is the zero matrix of size . Then for , we have
(22) 
now, since , we can deduce from (22), the following expression
which ends the proof of (17). The same proof can be done for the relation (18).
For the proof of (19), we first use (20) to obtain
which gives
it follows that
which ends the proof of the first relation of (19). In the same manner, we can prove the second relation.
Let be such that is invertible for . Let and have fullrank, where the . Let be a chosen tangential matrix directions. Then,

If for , then

If for , then

If both (1) and (2) hold and in addition we have , then,
1) We follow the same techniques as those given in [2] for the nonblock case. Define
and
It is easy to verify that and are projectors. Moreover, for all in a neighborhood of we have
and
Observe that
(23) 
Evaluating the expression (23) at and multiplying by from the right, yields the first assertion, and evaluating the same expression at and multiplying by from the left, yields the second assertion
2) Now if both (1) and (2) hold and , notice that
and
Therefore, evaluating (23) at , multiplying by and , from the left and the right respectively, for , we get
Now notice that since , we have
which proves the third assertion.
In the following theorem, we give the exact expression of the residual norms in a simplified and economical computational form.
4 An adaptive choice of the interpolation points and tangent directions
In the section, we will see how to chose the interpolation points , and tangential directions , , where . In this paper we adopted the adaptive approach, inspired by the work in [10]. For this approach, we extend our subspaces and by adding new blocks and defined as follows
(24) 
where the new interpolation point , and the new tangent direction , are computed as follows
(25) 
(26) 
Here is defined as the convex hull of where
are the eigenvalues of the matrix
. For solving the problem (25), we proceed as follows. First we compute the next interpolation point, by computing the norm of for each in , i.e we solve the following problem,(27) 
Then the tangent direction is computed by evaluating (25) at ,
(28) 
We can easily prove that the tangent matrix direction is given as
where the ’s are the right singular vectors corresponding to the largest singular values of the matrix . This approach of maximizing the residual norm, works efficiently for small to medium matrices, but cannot be used for large scale systems. To overcome this problem, we give the following proposition. Let and be the residuals given in (12) and (23), where and . Then we have the following new expressions
(29) 
and
(30) 
The residual can be written as
using Equation (17), we get
which gives,
The expression of given in (29) allows us to significantly reduce the computational cost while seeking the next pole and direction. In fact, applying the skinny QR decomposition
we get the simplified residual norm
(31) 
This means that, solving the problem (25) requires only the computation of matrices of size
Comments
There are no comments yet.