Task I.B Basic software tools for structured matrix decompositions and perturbations

A toolbox for structured matrix decompositions has been developed and included in the SLICOT package. The development of this toolbox was the subject of Task I.B of the European Community BRITE-EURAM III Thematic Networks Programme NICONET (contract number BRRT-CT97-5040). The toolbox includes routines for factoring a symmetric positive definite (block) Toeplitz matrix, or its inverse, and for solving associated systems of linear equations. In addition, routines for QR factorization of general (block) Toeplitz matrices, and for solving associated linear systems in a least squares sense are available. Only the first block row, and/or the first block column of the Toeplitz matrix must be given. In the development of the toolbox, computational reliability  and numerical efficiency  were our main concerns. The special structure is exploited, by determining and using the generators of the (block) Toeplitz matrix. The coding of all subroutines uses extensively the linear algebra package LAPACK and is done according to the SLICOT implementation and documentation standards.  All developed subroutines  have been extensively  tested  on various test examples and are fully documented. Their performance has also been evaluated, and compared with standard or structure-exploiting algorithms implemented in MATLAB.  For MATLAB and Scilab users' convenience, mex- and m-functions have been developed so that major functionality of SLICOT routines is available within these user-friendly environments.

Description of the Toolbox for Structured Matrix Decompositions

The toolbox subroutines

The Fortran routines for TASK I.B  are aimed to perform efficiently and reliably factorizations of structured matrices, specifically, block Toeplitz matrices , and solve the associated linear matrix equations, using structure-exploiting algorithms. The Toeplitz matrix could be defined by its first block row and/or its first block column. Previous factorizations can be efficiently updated for additional blocks. A list of the implemented Fortran routines grouped according to their main functionality and with links provided to the associated .html documentation is given in the following table:

 
MB02CD Computes the Cholesky factor and the generator and/or the Cholesky factor of the inverse of a symmetric positive definite block Toeplitz matrix 
MB02CU *Brings the first blocks of a generator in proper form (enhanced version of MB02CX)
MB02CV *Applies the transformations created by the routine MB02CU on other columns or rows of the generator
MB02CX *Brings the first blocks of a generator in proper form
MB02CY *Applies the transformations created by the routine MB02CX on other columns or rows of the generator
MB02DD Updates the Cholesky factor, the generator, and/or the Cholesky factor of the inverse of a symmetric positive definite block Toeplitz matrix, given the information from a previous factorization and additional blocks of its first block row, or its first block column
MB02ED Solves a system of linear equations, T X = B or X T = B, with a symmetric positive definite block Toeplitz matrix T
MB02FD Computes the incomplete Cholesky factor and the generator of a symmetric positive definite block Toeplitz matrix 
MB02GD Computes the Cholesky factor of a banded symmetric positive definite block Toeplitz matrix 
MB02HD Computes the Cholesky factor of the matrix T' T, with T a banded block Toeplitz matrix of full rank 
MB02ID Solves an over- or underdetermined linear system with a full rank block Toeplitz matrix 
MB02JD Computes a full QR factorization of a block Toeplitz matrix of full rank 
MB02JX Computes a low rank QR factorization with column pivoting of a block Toeplitz matrix 
MB02KD Computes the product C = alpha*op( T )*B + beta*C, with T a block Toeplitz matrix 

The documentation of routines is also accessible from the SLICOT Library main index (for the user-callable routines), or from the SLICOT Supporting Routines index (for the auxiliary routines, marked with * in the list). The SLICOT Supporting Routines index is also accessible from the main Library index.

The SLICOT library is available via the SLICOT homepage. All routines in the above list are included in the file slicot.tar.gz (or the PC version, slicotPC.zip), downloadable from the SLICOT homepage.
 

Mex- and m-functions

The mex- and m-functions implemented within TASK I.B  are aimed to provide convenient, efficient and reliable solutions to problems with (block) Toeplitz structure, in the user-friendly environment MATLAB. They are built on the SLICOT Library routines implemented within TASK I.B. Executable SLICOT mexfiles are provided for MATLAB running under WINDOWS 95/98/NT/ME/2000. Note that the Fortran source code of the mexfiles have been written using the Fortran 90 dynamic allocation feature. The following mex-functions have been implemented:
 
fstoep Factors a symmetric positive definite (block) Toeplitz matrix and/or solves associated linear systems (interface to MB02CD, MB02DD, MB02ED)
fstoeq Computes the QR factorization of a (block) Toeplitz matrix and/or solves associated linear least-squares systems (interface to MB02HD, MB02ID, MB02JD, MB02JX, MB02KD)
For user's convenience, several easy-to-use interface m-functions have been additionally implemented in the toolbox, explicitly addressing some of supported features. The following table contains the list of implemented m-functions:
 
fstchol Factors a symmetric positive definite (block) Toeplitz matrix T, and solves associated linear systems, given the first (block) row / column of T
fstgen Factors a symmetric positive definite (block) Toeplitz matrix T, computes the generator of its inverse, inv(T), and/or solves associated linear systems using the Cholesky factor of inv(T), given the first (block) row / column of T
fstsol Solves linear systems X T = B / T X = B, where T is a symmetric positive definite (block) Toeplitz matrix, given the first (block) row / column of T
fstupd Factors and/or updates a factorization of a symmetric positive definite (block) Toeplitz matrix T, and solves associated linear systems, given the first (block) row / column of T
fstqr Computes the orthogonal-triangular decomposition of a (block) Toeplitz matrix T and solves associated linear least-squares problems, given the first (block) row and (block) column of T
fstlsq Solves the linear least-squares problems min(B - T X) or finds the minimum norm solution of T' Y = C, where T is a (block) Toeplitz matrix with full column rank, given the first (block) column and the first (block) row of T
fstmul Computes the matrix-vector products x = T b for a (block) Toeplitz matrix T, given the first (block) column and the first (block) row of T

The executable mexfiles for MATLAB running under WINDOWS 95/98/NT/ME/2000,  together with the source codes, the above m-files and some MATLAB test scripts are available in the compressed archive  fstoep_mex.zip,  from the href="http://www.slicot.de" target="blank"">the SLICOT homepage. Similar interfaces are available in Scilab. Scilab functions are mostly compatible with the MATLAB interface and, in particular, the names of the corresponding mexfiles and m-files are the same.

A demonstration package fstoep_demo.zip for a PC machine is available from the SLICOT homepage. After decompressing that file in the current directory, the demo can be directly tried by typing fstdemo in a MATLAB 6.1 session (from that directory).

Performance Results

Extensive testing has been performed to assess the accuracy and efficiency of the toolbox components. The results have been compared with the available, standard techniques incorporated in MATLAB (namely, the m-function chol, and the linear system solving operators / and \), which do not exploit the problem structure, as well as with the MATLAB implementation of the same fast algorithms for block Toeplitz matrices (the m-function toepinv). A sample of the results can be seen below. Random data sets are used, but positivity is forced. The coefficient matrices are nxk-by-nxk symmetric block Toeplitz matrices BT, with k taking values in [1 2 20 30], and n obtained by dividing 300 to k, hence nxk = 300 ; k is the blocksize, and n is the number of the blocks. Also, the linear systems Y BT = C and BT X = B (C = B') are solved, the matrix B having nrhs columns, with nrhs also taking values in [1 2 20 30].  Click here to see the relative errors of the factorizations for the four problems with k = 1, 2, 20, and 30. The error values e_R, e_L, e_I, and e_C have been computed with the following formulas

e_R = norm(R'*R - BT,1)/norm(BT,1);
e_L = norm(L*BT*L' - BT,1);
e_I = norm(Ti*BT - eye(n*k));
e_C = norm(U'*U - BT,1)/norm(BT,1);

where R and L are the computed upper and lower Cholesky factors of BT and inv(BT), respectively, Ti is an approximate inverse of BT, obtained using its computed generator, G, and U = chol(BT).

To see the relative errors and residuals of the solutions of BT X = B, computed by fstoep and MATLAB \ operator, click here. For each value of k, all values of nrhs have been tried, hence 16 different systems have been solved. The errors and residuals for fstoep are slightly larger than for standard MATLAB functions, but essentially the same.

Click here to see the execution times for computing various factorization results using fstoep and MATLAB toepinv. Similarly, click here to see the execution times for solving the 16 systems mentioned above, BT X = B, or their transposed counterparts, Y BT = C. The speed-up factors for factorizations can be seen by clicking here, and the speed-up factors for solutions can be seen by clicking here.

Larger speed-up factors can be obtained for larger matrices, as can be seen here, for factorizations, and here, for solutions, when n*k = 1000 (k and nrhs take values in [1 2 20 50]). All the timing results have been obtained on a PC machine at 500 Mhz, without optimized BLAS.

Available Reports

1. P. Van Dooren "Selection of Basic Software Tools for Structured Matrix Decompositions and Perturbations", SLICOT Working Note 1999-9, June 1999.
2. D. Kressner and P. Van Dooren "Factorizations and linear system solvers for matrices with Toeplitz structure", SLICOT Working Note 2000-2, June 2000.

Vasile Sima,  February 16, 2001
Ad Van den Boom,  February 16, 2001
Updated: October 10, 2002;  January 17, 2004