The huge amount of theoretical results available in systems and control theory has led to a variety of methods and techniques used both in industry and academia. Although based on sound theoretical results, these methods often fail when applied to real-life problems, which frequently tend to be ill-posed or highly dimensional. Such failures are often due to the lack of numerical robustness when implemented in finite-precision computations.
SLICOT Library (Subroutine Library in Control Theory) has been designed and developed with the main aim of achieving the robustness and efficiency needed in solving possibly ill-conditioned, ill-scaled, and/or large-scale real-life problems. A special emphasis has been put on structure-preserving algorithms. The main advantage of such algorithms is that the structural properties of a problem are preserved during computations. This allows the computed result to be interpreted as the exact solution of the original problem with perturbed input data, which may otherwise not be the case. This not only increases the reliability of the returned results, but often also improves their accuracy.

The following results illustrate the performance (in terms of efficiency, reliability, and accuracy) which can be obtained using SLICOT components, in comparison with equivalent computations performed by some MATLAB 5.3 functions included in the MATLAB nucleus or in the MATLAB Control Toolbox. Results obtained on newer MATLAB versions are published in related NICONET reports and papers authored by the NICONET team. The results show that some SLICOT computations are several times faster than MATLAB computations, at comparable or better accuracy; moreover, less memory is required by SLICOT routines, because the problem structure is fully exploited whenever possible. The calculations have been performed on a SUN Ultra 2 Creator 2200 workstation with 128 M RAM under SunOS 5.5, by calling from MATLAB the SLICOT gateways produced by the NAGWare Gateway Generator of the Numerical Algorithms Group Ltd. (NAG). SLICOT routines have been compiled with f77 compiler using options -O4 -native -u. SUN performance libraries have not been used for obtaining the results. The data have been randomly generated from a uniform distribution in the interval (0,1). The timing values reported are the execution times in seconds exclusively spent in the corresponding SLICOT or MATLAB code portions.

Fast Fourier transforms

Tables 1 and 2 display comparative results for the SLICOT Fast Fourier transform routines DG01MD and DG01ND (for complex and real sequences, respectively), and the corresponding MATLAB function fft. Random sequences X of length n were generated. Besides better efficiency, the accuracy of SLICOT routines shows an improvement which can be two orders of magnitude. The accuracy was measured by computing the distance between the original sequence and the inverse Fourier transform of the transformed sequence, norm(X - ifft( fft( X ) ))/norm(X), which should theoretically be zero, and similarly for the SLICOT routines. Note that DG01ND is 2-4 times faster than MATLAB.

Matrix exponential

Tables 3 and 4 display comparative results for the SLICOT subroutines MB05MD and MB05OD, for the computation of the matrix exponential, and the equivalent MATLAB functions expm3 and expm1, respectively. The relative errors reported are those corresponding to the differences between the results obtained from SLICOT and MATLAB. SLICOT codes are faster than MATLAB calculations, at a comparable accuracy. An increase of over 30 times in speed has been achieved for n = 256 for upper triangular matrices. This remarkable speeding up was possible because the balancing strategy implemented in MB05OD can reduce the matrix norm for triangular matrices and save matrix multiplications, while expm1 cannot.

A = triu(rand(n,n),n/2)

Controllability staircase form

Tables 5 and 6 demonstrate the behaviour of the SLICOT routine AB01ND for computing the controllability staircase form of a linear system in state-space form given by a matrix pair (A,B), with A n-by-n, B n-by-m, and MATLAB function ctrbf provided by the Control Toolbox. In Table 6, m is taken as 1. It should be mentioned that ctrbf uses singular value decompositions for rank determination while AB01ND employs rank-revealing QR factorizations. The relative residuals have been computed using the following formulas
max(norm(z'*A*z - a)/norm(A), norm(z'*B - b)/norm(B));
max(norm(A - z*a*z')/norm(A), norm(B - z*b)/norm(B));
for SLICOT and MATLAB, respectively, where A and B are the given system matrices, a and b are the computed matrices in controllability staircase form, and z is the computed orthogonal transformation matrix. The SLICOT routine is up to ten times faster than ctrbf in the multi-input case and up to 60 times faster in the single-input case.

Lyapunov and Riccati equations

Tables 7 and 8 give performance results for SLICOT's general Lyapunov solver SB03MD, which solves both continuous- and discrete-time equations, and MATLAB functions lyap and dlyap. The SLICOT routine can solve equations defined by
op(A)' X + X op(A) = sigma C,(1)
op(A)' X op(A) - X = sigma C,(2)
where op(A) = A or A' (the transpose of A), and sigma is a scale factor, set less than or equal to 1 to avoid overflow in computing X. Moreover, an estimate of the condition number is returned. Specifically, Table 7 reports results for continuous-time problems with op(A) = A, generated using the MATLAB statements
A = rand(n,n); C = rand(n,n); C = C + C';
so the solutions were not known. Therefore, only relative residuals are given in Table 7.

Similar results have been obtained for the transposed case, or for discrete-time equations. SLICOT residuals were usually less then MATLAB residuals in our tests. The speed-up has been even larger when compared with the previous MATLAB 4.2 versions of lyap and dlyap.
Table 8 reports results for discrete-time problems with op(A) = A'. The problems were generated using the MATLAB statements
A = rand(n,n); X0 = rand(n,n) + n*eye(n); X0 = X0 + X0';
C = A*X0*A' - X0; C = (C + C')/2;
so the relative errors could be computed.

Table 9 gives performance results for the SLICOT Riccati solver SB02OD, which solves both continuous- and discrete-time equations, and the MATLAB functions care and dare, all based on the generalized Schur vector approaches. The SLICOT solver can also compute the anti-stabilizing solutions. The problems solved were generated using the MATLAB statements
A = rand(n,n); B = rand(n,m); G = B*inv(R)*B'; G = (G + G')/2;
(with m = n), and the weighting matrices R and Q have been computed to have given random (non-negative) eigenvalues, using a gateway for the test routine DLATMS from LAPACK. Table 9 reports the timings and relative residuals for computing stabilizing solutions for continuous-time problems.

Table 10 gives performance results for the SLICOT Riccati solver SB02MD, which solves both continuous- and discrete-time equations, using the Schur vectors approach. For discrete-time problems, SB02MD includes an option to use the inverse of the symplectic 2n-by-2n matrix, which is always faster than the standard approach.

Discrete-time problems can be solved even faster by SB02MD, but the accuracy depends on the conditioning of the matrix A.

Structure-Preserving Algorithms

Hankel singular values

The use of structure-preserving algorithms not only increases the reliability of the returned results, but often also improves their accuracy, as shown below for an important computational problem.
The controllability and observability Gramians Pc, Po, respectively, of a stable state-space realization (A,B,C) of a continuous-time linear time-invariant system, where A is n-by-n, B is n-by-m, and C is p-by-n, are given by the solutions of the stable Lyapunov equations
A Pc + Pc A' = - B B' (3)
A' Po + Po A = - C' C (4)
which are non-negative definite. The nonnegative square roots of the eigenvalues of the product PcPo, called the Hankel singular values, play a fundamental role in finding balanced realizations and in model reduction. To guarantee the symmetry and definiteness of the computed Gramians, the special problem structure should be taken into account. Besides this, exploiting the structure usually results in a reduction of necessary computational operations and memory requirements. The SLICOT Library routine SB03OD solves for X = op(U)' op(U) either continuous-time or discrete-time Lyapunov equations,
op(A)' X + X op(A) = -sigma*sigma op(M)' op(M) (5)
op(A)' X op(A) - X = -sigma*sigma op(M)' op(M) (6)
where A is square and stable (in the continuous- or discrete-time sense, respectively), M is rectangular, U is upper triangular, and sigma is a scale factor, set less than or equal to 1 to avoid overflow in computing X. The routine uses the same real Schur form of A for all different computations, and optionally A can be given in such a form on input.
The data have been generated as
A = rand(n,n) - n*eye(n); B = rand(n,m); C = rand(l,n);
Table 11 gives comparative results using the gateway for the SLICOT routine SB03OD and MATLAB function lyap, for solving the two Lyapunov equations (3) and (4). It is apparent that the speed-up compared with MATLAB is about 5. Besides improved efficiency, the accuracy of the SLICOT routine is sometimes better. More important, the Hankel singular values can be obtained as the singular values of the triangular matrix V U, where U and V are the Cholesky factors computed by SB03OD for the solutions of the equations (5), for op(K) = K' and op(K) = K, respectively. On the contrary, MATLAB calculations need to be based on sqrt(eig(X*Y)), where X = Pc and Y = Po are the solutions of the Lyapunov equations (3) and (4), respectively.

Similar results have been obtained for the relative residuals in Y.
For one sample of the largest example, two computed eigenvalues of the matrix X Y were negative (close to the machine accuracy), and four were complex conjugate. In other words, in extreme cases, MATLAB could give physically meaningless negative or complex Hankel "singular values". This, of course, does not happen for well-conditioned problems.

Eigenvalues of Hamiltonian matrices

Table 12 displays timing results for the SLICOT Hamiltonian eigenvalue solver based on the square-reduced form (MB04ZD and MB03SD) and MATLAB function eig. There is no MATLAB function exploiting the Hamiltonian structure.

For random matrices, MATLAB often gave eigenvalues very close to the imaginary axis, and which do not not appear in pairs lambda, and -lambda, as the structure ensures. An example is given below. Let G and Q be symmetric matrices, and consider the MATLAB statements
A = [ ...
8.156785069508218e-01 9.583801041682488e-01 3.592539466930961e-01
3.786356003299222e-01 3.034683538610545e-01 7.876353002957823e-01
4.773517592106414e-01 8.360131178618799e-01 7.242888156356292e-01 ];

Q = [ ...
8.890570332721912e-02 9.990740343443876e-01 8.732331409873815e-01
9.990740343443876e-01 6.033984620815241e-01 5.861126862288001e-01
8.732331409873815e-01 5.861126862288001e-01 7.904409030391265e-01 ];

G = [ ...
5.951660134153334e-01 9.187384627530358e-01 6.419067121426090e-01
9.187384627530358e-01 2.943561392390958e-01 4.831226444677942e-01
6.419067121426090e-01 4.831226444677942e-01 7.635187432100927e-01 ];

H = [ A G; Q -A'];

eig(H)

ans =
2.734666984279399e+00
-2.734666984279397e+00
7.747532918205485e-01
-7.747532918205493e-01
6.938893903907228e-17 + 2.317956004861914e-01i
6.938893903907228e-17 - 2.317956004861914e-01i

Condition estimates for Lyapunov and Riccati equations

Many linear systems analysis and design problems (for instance, stability analysis, transfer-function matrix factorization, linear-quadratic or robust optimization problems), involve the solution of either continuous-time or discrete-time matrix Lyapunov and Riccati equations. In such applications, it is very useful to have easily computable estimates of the problem conditioning and error bounds for the solutions. Such estimators have been developed and are included in the SLICOT Library.
Table 13 displays estimates for the "separation" of the matrices A' and A, sepd(A',A), which is needed to compute the condition number for discrete-time matrix Lyapunov and Riccati equations. Similarly, the separation of the matrices A' and -A, sep(A',-A), is used for continuous-time equations. For A an n-by-n matrix with real elements, sep(.,.) and sepd(.,.) are defined in terms of the smallest singular value of a large matrix, by

sep(A', -A) = sigma_min (P), P = kron(I_n, A') + kron(A', I_n),
sepd(A', A) = sigma_min (P), P = kron(A', A') - I_n*n,

where kron(X,Y) denotes the Kronecker product of the matrices A and B, kron(X,Y) = (x_{ij} Y). The evaluation of the norms above by computing the minimal singular value sigma_min is not practical for medium-size or large-scale problems, due to the excessive computational effort and memory requirements. SLICOT routines exploit the problem structure.
The results in Table 13 have been computed by two approaches:
(a) solution of original Lyapunov equations;
(b) solution of reduced Lyapunov equations,
as well as by MATLAB function svd. The best results were obtained for approach (b), which can be two-four times faster than approach (a), and possibly thousands times faster than MATLAB svd.

Conclusions

SLICOT computations are several times faster than MATLAB computations, at comparable or better accuracy. Moreover, the problem structure is better exploited by SLICOT routines. Even better efficiency is attainable by using machine-tailored BLAS libraries, which, together with LAPACK subroutines provide the basic algorithmic level of SLICOT, and enable to exploit the potential of modern high-performance computer architectures.

This email address is being protected from spambots. You need JavaScript enabled to view it., April 7, 1999