Consider a linear time-invariant multi-input, multi-output linear time independent (LTI) dynamical system described by the state-space equations
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
The goal of our model reduction approach consists in defining two orthogonal matrices and (with ) to produce a much smaller order system with the state-space form
and its transfer function is defined by
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 , optimal Hankel norm  and Krylov subspace methods, [8, 9, 14, 21] have been used for large multi-input multi-output (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 large-scale 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 Lanczos-type 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 Lanczos-type algorithm. Section 5 we treat the model reduction of second-order 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
where is called the
-th moments atassociated to the system and defined as follows
In the case where the moments are called Markov parameters and are given by
Problem: Given a full-order model (1) and assume that the following parameters are given:
Left interpolation points & block tangent directions .
Right interpolation points & block tangent directions .
The interpolation points and tangent directions are selected to realize the model reduction goals described later.
3 The block tangential Lanczos-type method
Let the original transfer function be expressed as where and are such that,
Given a system of matrices and where , the approximate solution and of and are computed such that
which gives the following approximate transfer function
where , and . The matrices and are bi-orthonormal, where the . Notice that the residuals can be expressed as
3.1 Tangential block Lanczos-type algorithm
This algorithm consists in constructing two bi-orthonormal bases, spanned by the columns of and , of the following block tangential Krylov subspaces
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 bi-orthogonality conditions for :
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 bi-orthonormal matrices of constructed by Algorithm 1. Then we have the following relations
Let and be the matrices,
then we have
where and . The matrices and are block upper triangular matrices of sizes and are assumed to be non-singular.
From Algorithm 1, we have
multiplying (20) on the left by and re-arranging terms, we get
that written as
where 0 is the zero matrix of size . Then for , we have
now, since , we can deduce from (22), the following expression
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 full-rank, 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  for the non-block case. Define
It is easy to verify that and are projectors. Moreover, for all in a neighborhood of we have
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
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 . For this approach, we extend our subspaces and by adding new blocks and defined as follows
where the new interpolation point , and the new tangent direction , are computed as follows
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,
Then the tangent direction is computed by evaluating (25) at ,
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
The residual can be written as
using Equation (17), we get
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
This means that, solving the problem (25) requires only the computation of matrices of size