MINISTERO DELL'UNIVERSITÀ E DELLA RICERCA SCIENTIFICA E TECNOLOGICA
PROGRAMMI DI RICERCA ANNO 2004
COMPITI E SUDDIVISIONE DELLE UNITÀ DI RICERCA

 



 

 

 Purpose 

This project is the continuation of a previous MIUR project 2002-2003, its goal is to enlarge the range of structured tools and applications under investigation by expanding and integrating the results obtained in 2002-2003 (see http://bezout.dm.unipi.it). The aim of the project consists of

A- Analysis of algebraic, analytic and computational properties of classes of structured matrices encountered in theoretical and in applied problems, including linear and nonlinear structures. In particular, displacement structures (Toeplitz, Hankel, Cauchy, Bezout, etc), semi-separable matrices, multi-level (multi-index) matrices, locally structured matrices, parametric systems.

B- Design of fast and numerically reliable (stable and robust) algorithms for solving computational problems involving structured matrices, like total least squares, linear systems, eigenvalue computations. Adaptation of known techniques like multigrid, PCG, GMRES, QR iteration, to specific structures. Use of structured tools for locally structured or even unstructured (minimization) problems.

C- Application of the results to the solution of computational problems from scientific computing, including: image processing, nonlinear inverse problems, queueing models, matrix equations, integral and differential equations, inverse scattering problems, approximating polynomial zeros, computer aided geometric design, human metabolism models, information retrieval, neural networks.

D- Implementation of the algorithms in a prototype code, numerical comparisons and validation.

There are several international research groups which are very active in the field. The Italian group is actively working in the area since many years. There is a strong theoretical and algorithmic basis made up by many tools and powerful techniques developed in the years by many researchers. There are different research lines where the interest is growing up at the international level. Some new lines stemmed from the previous MIUR project.

 

Evaluation

1) The evaluation of the scientific production will be based on the number of contributions as well as on the importance of the journals where the contributions will appear.

2) The evaluation of the international relevance of the results will be based on the number as well as on the importance of the conferences where the results will be presented.

3) Algorithmic improvements. The evaluation of the improvements is given by the reduction of complexity and/or by the reduction of the cpu time required by the algorithms once they are implemented. Numerical stability is also a decisive issue.

 

Research Unities
 

 

    Unity of PISA

    Leader DARIO ANDREA BINI

 

 

    Task


Large part of important problems in pure and applied mathematics involves structured matrices. The analysis of theoretical and computational properties of these structures is a fundamental step in the design of efficient algorithms specifically designed for these problems especially in the very frequent cases where the general solution algorithms cannot be applied due to the huge size of the problem.

Matrix structures are inherited by the specific properties of the problem. Common features shared by problems in different areas, like shift invariance, compact support, separability of variables, turn into structures which are common to matrices encountered in different fields, like the Toeplitz structure, the displacement structure (including Hankel, Cauchy, Bezout, Frobenius, Pick, Vandermonde matrices), band matrices, semi-separable matrices, arrowhead matrices etc. In fact, these structures are encountered in many different problems in statistics, probability, queueing theory (including telecommunications, internet, performance analysis, stochastic processes), polynomial computations, computer algebra, Computer Aided Geometric Design, numerical solution of differential and of integral equations, image and signal processing just to quote a few.

Very often the properties of the original problem turn into structures which are not so evident, say, nonlinear structures like displacement structure, semi-separable or quasi-separable matrices, where the matrix under analysis looks like a "general" matrix.

The analysis of matrix structures and the consequent design of specific algorithms for solving problems in numerical linear algebra has received a great attention and has produced important advances in the algorithm design. On one hand, important theoretical advances have been reached in the study of theoretical properties of certain structured matrices, on the other hand the design of specific algorithms has allowed the effective solution of important problems in different applicative fields. Just to quote a few, algorithms for solving matrix equations, or for solving queueing models based on Toeplitz computations or on Wiener-Hopf factorization, recently designed, allowed to reduce the cpu time of a large factor still maintaining very good numerical performances [ALM]. In certain cases, the structured matrix tools have provided fast solution of unstructured problems [D], [DFLZ], [DLZ].

Summing up, the research of this unit are dictated by the following general guidelines and motivations:

1- Structured matrices play a fundamental role in a great part of mathematical models in applications and are almost pervasive.
2- The analysis of structured matrices, besides its theoretical interest, is a fundamental step for the design of specific solution algorithms. The results obtained in this way constitute a general toolbox which is useful in many different contexts and situations and constitutes the structured matrix technology.
3- The structured matrix technology can be used to effectively solve large scale problems of great importance in the applications, which could be hardly solved with general algorithms. At the same time, even for the solution of general unstructured problems, the structured matrix technology may provide important improvements.
4- The design and the analysis of algorithms taylored on the specific structures of the problem is not enough in itself. It must be integrated by a rigorous analysis of robustness and of numerical stability. In certain situations where the ill conditioning is an intrinsic feature of the problem, the use of multi-precision arithmetic or the use of hybrid approaches based on symbolic and numeric tools becomes a required step.


GENERAL RESARCH LINES

The research activity will be carried out according to the following research lines:

-A- To perform the analysis of algebraic, structural, analytic, spectral and computational properties of structured matrices which intervene in theoretical and applied mathematics (including Toeplitz, Hankel Bezout, Cauchy, Frobenius, and displacement generated matrices, banded matrices, semi-separable and quasi-separable matrices, arrowhead matrices, matrices with scalar entries and with blocks, multilevel matrices).
-B- To exploit the properties obtained from this analysis for creating general methodologies for the design and analysis of effective solution algorithms.
-C- To apply the solution algorithms to problems coming from certain applications and from large scale scientific computing.
-D- To implement, in terms of prototype software, the algorithms designed in this research and to perform experimental validation and comparison.
-E- To perform a theoretical or experimental rounding error analysis and robustness analysis in order to select reliable algorithms. For problems endowed of intrinsic ill conditioning, to design new techniques based on the use of multiple precision arithmetic complemented with the strategy of adaptiveness, and integrating numerical and symbolic tools into hybrid algorithms.

SPECIFIC TOPICS:

Specific topics of our research are:

1- Polynomial computations: Analysis of (weak) Wiener-Hopf factorizations where both the factors may have zeros of unit modulus. Polynomial root-finders based on the QR algorithm applied to generalized companion matrices which maintain the semi-separable (or quasi-separable) structure. Polynomial rootfinders based on Newton's and Laguerre's iteration whith Gerschgorin-like inclusion theorems and the concepts of pseudozeros. Using the representation of polynomials in the Bernstein basis together with their structured matrix counterpart in order to design efficient algorithms for the intersection of Bezier curves and of surfaces. Using hybrid techniques for the related computational problems with tools in projective geometry, computer algebra and numerical analysis. Computing eigenvalues of tridiagonal matrices.
2- Matrix equations: the aim is designing algorithms for solving matrix equations based on Wiener-Hopf factorizations. Algorithms for the Riccati equation, algorithms for computing the principal p-th root of a matrix.
3- Structural and computational properties of semi-separable, quasi-separable matrices, Toeplitz matrices and matrices in trigonometric algebras: one of the aim is designing algorithms for the computation of eigenvalues and for polynomial rootfinders (see part 1)
4- Minimization problems: design and analysis of quasi-Newton algorithms for minimization problems where the Hessian matrix is approximated by a structured matrix. Solving total least squares problems.
5- Applications of part 1) to solving Computer Aided Geometric Design problems like intersection of curves and of surfaces parametrically or implicitly assigned; application of part 2) to solving Markov chains encountered in queueing models; application of part 3) to the design of algorithms for polynomial rootfinding; application of part 4) to solving neural networks.

MAIN TASKS

Concerning part 1) we wish to generalize and extend the shift techniques introduced in [HMR], [BM7] for accelerating the convergence of solving certain equations in Markov chains, in order to deal with more general cases where the Wiener-Hopf factors may have zeros of modulus 1.
-- We wish also to investigate more closely on the algorithm designed in [BDG] where the QR iteration applied to an nxn Frobenius matrix is performed in O(n) ops per step. Here the main issues concern numerical stability and robustness. Also the algorithm introduced in [BGP] for the computation of the eigenvalues of a subclass of quasi-separable matrices will be closely investigated. Here the goal is to obtain an efficient polynomial rootfinder and investigate on a possible generalization of the class of matrices for which the QR iteration is closed.
-- Polynomial rootfinders will be analyzed also in terms of simultaneous and independent iterations based on the result of [HSS] trying to replace Newton's iteration with the modified Laguerre iteration and complementing it with some Gerschgorin like inclusion theorems for the approximations provided by the algorithm.
-- Bernstein polynomials and Bezier curves will be investigated as well. Here the main goal is to perform algebraic operations with polynomials represented in the Bernstein basis without changing to the power basis. In fact this transformation is severely ill conditioned. This goal requires that the algebraic operations which are usally performed in the power basis in terms of structured matrices (Bezout and Sylvester matrices) be performed in the Bernstein basis. Bernstein-Bezoutian metrices must be introduced and accurately analyzed. Moreover algorithms for the LU or QR decomposition of such matrices must be designed in such a way that either their numerical implementation is numerically stable or their symbolic implementation does not lead to a huge growth of the intermediate integers.
In this context, also the computation of the zeros of polynomials represented in the Bernstein basis plays an important role and will be investigated. Another related important topic is the application and integration of numerical tools, symbolic tools, computer algebra and projective geometry for solving computer aided geometric design problems [Gianni].

Concerning part 2) we continue our study of [BM], [BGM1],[BGM2],[BG2] on solving matrix equations. An important advance in this field was introduced in [HMR], [BM7] in the context of Markov chains by means of the "shift technique". Here we want to generalize this technique to general matrix equations of polynomial type.
-- We want to apply the cyclic reduction algorithm of G. Golub [BGN] adjusted to Markov chains by [BM3],[BM4] and analyzed in [BGM1],[BGM2] to solving algebraic Riccati equations. Also the computation of the principal p-th root of a matrix will be faced in terms of Wiener-Hopf factorizations, Laurent series inversion and Newton's iteration.
-- Any kind of matrix equations modeling queueing problems will be taken into account.
-- The problem of computing the p-th root of a matrix will be closely investigated by using the techniques of [M4] and of [Ia].

Concerning part 3, we focus our attention on properties of semi-separable and quasi-separable matrices with the aim of extending the algorithms of [BDG] and [BGP] for implementing the QR iteration in linear cost, to more general classes of matrices. This part is related to some polynomial computations of part 1). Concerning Toeplitz matrices and trigonometric algebras the interest is addressed to the analysis of algebraic multigrid methods for solving the related systems.

Concerning part 4), we continue the study of [D], [DFLZ], [DLZ] with the analysis of sufficient conditions which guarantee the convergence of the quasi-Newton minimization techniques where the Hessian matrix is replaced with a matrix in the Hartley algebra [BF] or with some structured matrix. Here the conditions will concern the boundedness of the operators appearing in the minimization process and on the condition number of the structured matrix. Also adaptive methods where the approximation of the Hessian is performed at each step of the iteration will be studied. Algorithms for total least squares problems for structured matrices will be analyzed with applications to blind image deblurring.


Concerning part 5) of applications, we will apply our algorithms for the listed applications. In particular, the results of 1) concerning Bernstein-Bezout matrices, Bernstein polynomials and hybrid algorithms will be applied to solving CAGD problems. We will apply the results of 2) to solving queueing models where the solution of the associated matrix equations is the most expensive task, we will consider some queueing models which are somehow related to the algebraic Riccati equation. Concerning applications to Markov chains and to queueing models, we planned to organize in Pisa in June 2005 the fifth edition of the international conference "Matrix Analytic Methods in Stochastic Models (MAM5)" which will be supported by this project. We will apply the results of 4) to solving neural networks problems which model automatic learning and to blind image restoration.


Part of the research in 3) on semi-separable and quasi-separable matrices and on multigrid and on 4) on image deblurring is performed in collaboration with the Unit of Genova.

[BG2] D.Bini,L.Gemignani, Solving quadratic matrix equations and factoring polynomials: new fixed point iterations based on Schur complements of Toeplitz matrices. To appear in Numer. Lin. Alg. Appl. 2004.
[BGFM] D.Bini, G.Fiorentino,,L.Gemignani,B.Meini, Effective fast algorithms for polynomial spectral factorization , Numer. Algo., 2003.
[BGP] D.Bini,L.Gemignani,V.Pan $QR$-like algorithms for generalized semiseparable matrices T.R. 1470, Dept. of Math., University of Pisa, 2003.
[BGT] D.Bini L.Gemignani F.Tisseur. The Ehrlich-Aberth method for the nonsymmetric tridiagonal eigenvalue problem. TR 428, Univ. of Manchester, England 2003
[BM7] D.Bini,B.Meini, Non-Skip-Free M/G/1-type Markov chains and Laurent matrix power series, Proc. "Numerical Solution of Markov Chains", Urbana, September 2003.
[FGPT] E.Fortuna, P.Gianni, P.Parenti, C.Traverso, Algorithms to compute the topology of orientable real algebraic surfaces, J. Symb. Comp., 36,2003
[HMR] C.He, B.Meini, N.Rhee, A shifted cyclic reduction algorithm for QBDs, SIMAX,23 2001/02.
[HSS] J.Hubbard,D.Schleicher,S.Sutherland, How to find all roots of complex polynomials by Newton's method. Invent. Math. 146 (2001).
[ALM] G.Anastasi, L.Lenzini, B.Meini, Performance evaluation of a worst case model of the MetaRing MAC protocol with global fairness. Performance Evaluation, 29:127--151, 1997.

 

 

     Unity of GENOVA

    Leader FABIO DI BENEDETTO

 

 

     Task


The new program is organized along two directions that were the main themes considered in the project 2002-2003:
A) Theoretical and computational properties of structured matrices;
B) Structured matrix tools in applied problems.
The topics considered in A) are open problems in structured numerical linear algebra which represent natural continuations of investigations performed in the previous project. Special emphasis will be devoted to less investigated structures (such as shifted, semi-separable matrices, new matrix algebras related to the sine transform of type I, "spectral fattening" of circulant-like/trigonometric algebras): the first three structures have received particular attention in the most recent literature and in applications, while the last is a really new proposal trying to combine the idea of circulant-like and sparse approximate inverse preconditioners.
The main applications considered in B) are again functional approximation, inverse problems, PDEs. Here, in addition to the continuation of previous researches, the group will investigate new applied problems and/or some completely new approaches, often based on advanced tools that are recently developed in the context of structured matrices.
The essential difference with respect to the 2002-2003 project concerns the weights of themes A) and B): more emphasis will be given to (real world) applications. A detailed description of the specific points of the program is reported below.

A1) Theoretical issues
a) Asymptotical linear algebra
The powerful results obtained so far will be applied to the convergence analysis of iterative methods. In particular, if matrices show some distribution property for the eigenvalues and a statistical approach to the error analysis is considered, we want to investigate the role of the distribution function: we believe that its behavior is responsible to this statistical behavior as the spectral radius is in the deterministic setting.
The picture concerning the asymptotical spectral analysis (distribution, localization, extremal behavior) of structured matrices of Toeplitz and Locally Toeplitz type (including PDEs with variable coefficients and on involved geometries) is quite complete now: it remains to understand the extremal behavior (spectral conditioning) in the multidimensional setting following the approach of Bottcher and Grudsky.
In a more non-structured framework, due to interest in convergence theory of Krylov methods, we will study how to deduce spectral clustering of the eigenvalues (proper or general) from a spectral clustering of the singular values for which much more tools are avalaible.
b) Relations between structures
The interplay between diagonal-plus-semiseparable matrices and rational Krylov matrices, firstly outlined in [28], is a rational counterpart of the well-known connection between tridiagonal matrices and classical Krylov matrices. We intend to further this study in order to analyze properties of rational variants of the Krylov iteration, and to deepen the theoretical and computational properties of orthogonal rational functions.
A2) Numerical methods
a) Parametric matrices
A sequence of linear systems is called parametric if it can be represented as {Aj xj = bj} where Aj=A+vj Ej and vj is a scalar parameter; this description generalizes the simpler case of shifted matrices. We plan to focus our attention to incomplete updated factorizations for parametric sequences of linear systems such that
i) A nonsymmetric positive definite and Ej=I (shifted case), and
ii) A symmetric (definite and indefinite) and Ej is banded,
cases which are crucial in the related applications.
For instance, shift-and-invert (SI), inverse power (IP) and Jacobi-Davidson (JD) are popular techniques for the computation of some (typically few) eigenvalues and eigenvectors for large and sparse Hermitian problems. We stress that the core of these algorithms require the efficient solution of several large (and usually indefinite) algebraic systems in shifted form. We plan to follow two main paths to solve the underlying problems:
(i) ad-hoc indefinite incomplete updated factorization preconditioners used with a suitable Krylov projection method;
(ii) solution of shifted systems by projection techniques and multilevel incomplete factorization preconditioned algorithms, with reorderings to transform the original problem in a more diagonally dominant one.
Another application involves the solution of large and sparse total least squares (TLS) problems based on Rayleigh quotient iteration. In this approach, the solution of the underlying problem reduces to the solution of a sequence of shifted symmetric positive definite linear systems. In [14] Björck et al. proposed solving these linear systems by the conjugate gradient method preconditioned by the (complete) Choleski factorization of the normal system matrix. We plan to design and analyze iterative algorithms which generate sequences of linear systems whose condition numbers are better, by using techniques based on incomplete updated and updated approximate inverse factorizations.
b) Semiseparable structures
Matrices with generalized semiseparable structures are gaining a great interest as computational tool to solve eigenvalue problems. The recent research in this field has opened new perspectives for what concerns the design of efficient numerical algorithms for special eigenvalue problems, e.g., associated with generalized companion or quasi-unitary Hessenberg matrices.
Moreover, we want to understand the computational relevance of prescribing the diagonal term in the similarity transformation of a symmetric matrix into diagonal-plus-semiseparable form, in order to speed up the computation of its eigenvalues.
c) Multilevel structures
We are interested in the analysis and in the synthesis of fast algorithms for the solution of multilevel structured linear systems (analysis of preconditioners for the conjugate gradient method, of the associated properties of clustering and spectral equivalence, and analysis of multigrid methods with a special emphasis in the choice of the projector/smoother pair). More in detail, we will be interested in the following aspects:
- Preconditioning of multilevel Toeplitz structures, trying to generalize to the multivariate case the optimality results given in [40].
- Investigation of a ''fattening'' of algebras, by taking matrices of the form U T U* where T is a matrix with a given sparsity pattern and for which the solution of an associated linear system has a complexity bounded by the cost of multiplying U or U* by a vector (typically O(n log n) operations). We would like to understand if this new proposal, which sensibly enriches the ordinary algebras where T is just diagonal, is sufficient for finding preconditioners which are superlinear and/or optimal in the multilevel polynomially ill-conditioned case (which is impossible if T is diagonal).
- Extension of the multigrid optimality proof in [2] to the multilevel case and applications to image restoration problems (see B2-a) and to the PDEs (especially in connection with the dense structured systems arising from the Sinc-Galerkin methods applied to elliptic PDEs).
B1) Continuation of previous research topics
a) Approximation by rational functions
We want to pursue our previous research on interpolation and approximation problems with rational functions, particularly with respect to the analysis of structured conditioning issues and optimal pole placements.
b) Inverse problems in imaging
If we use iterative solvers for linear inverse problems, the iteration number represents a regularization parameter; therefore a first goal is the automatic estimate of such parameter. Another direction involves the extension of regularizing families of preconditioners (e.g. filtering optimal and superoptimal, as described in the Scientific Background) to the reconstruction of non-negative images, and their application to problems related to the adaptive optics device equipping the LBT system.
We also plan the analysis and design of new efficient algorithms for linear systems generated by Tikhonov-like regularization. We would restrict our attention to regularization operators representing the identity, the first or the second derivative. The related normal equations are of parametric form, and should be solved for several values of the Tikhonov parameter. The proposed algorithms are new iterative methods using updated incomplete factorizations generalizing those introduced in [3] (which can be applied directly if the least squares are penalized by means of the solution norm).
Finally, a large amount of work is necessary to optimize the image reconstruction algorithms (i.e. OS-EM method for multiple images). For instance, the LBT restoration problem will require the treatment of images with very high-dynamic range (young binary stars or extrasolar planets). In such a case, the problem will be investigated along two different lines. The first one concerns the study of an appropriate functional in order to take into account the particular features of the object. The second one is a sophisticated post-reconstruction denoising in order to preserve both huge gradients and weak structures.
c) PDEs
The convergence properties of the preconditioned GMRES solution will be investigated for linear systems discretizing the convection-diffusion equation with variable diffusion coefficient. The preconditioner is based on the diffusive part of the operator and we plan to consider a finite difference discretization on a d-dimensional rectangular domain. We plan to show the existence of a proper cluster for the eigenvalues of the preconditioned matrix and then the superlinear behavior of iterations. A preliminary analysis has shown that the structure of the cluster is not influenced by the discretization but the number of outliers linearly increases with respect to the viscosity parameter.
The same techniques will be applied to time-dependent convection-diffusion equations arising in a human metabolism model, which has been introduced at the MIMS center of Cleveland.
B2) New techniques
a) Regularization by multigrid
It is well known that for retrieving the real image from the observed one corrupted by blurring and by noise, it is necessary to resort to some regularizing methods. There exist many procedures in this direction as Tikhonov, Riley, iterative methods as CG, CGNE, Landweber etc (see [25]). However often the most precise ones are also terribly slow. We would like to investigate the possibility of using a multigrid technique (based on the algebraic approach in [2,47]) for the regularization problem as well. Our goals are the following:
- obtain the same reconstruction quality as the best regularizing methods but with a great saving in computation time;
- adapt the method to every boundary condition (see the subsequent point b)) because the algebraic method described in [2] is very versatile with respect to the different structures.
b) Boundary condition effects
We concentrate our efforts in the direction of anti-reflective BCs by considering two goals:
- a more precise analysis of the (non-canonical i.e. mixed structures) algebra of matrices appearing with the anti-reflective BCs;
- the study of the role of the noise and more precisely how to use a regularizing method (which is necessary due to the noise and to the ill-posedness of the problem) in connection with the anti-reflective BCs: we would like to modify the classical normal equations and Tikhonov techniques in such a way that the structure of algebra of the resulting system is exploited in the best possible way. It is understood that the main aim is to obtain a better precision with respect to regularizing methods applied with the other BCs.
c) Denoising and edge detection by wavelets
Astronomical images are corrupted by a large amount of noise. Denoising methods are required in many cases, and in particular we will investigate three cases: chopped and nodded images (near-infrared images), measured PSFs used in classical restoration methods, over-iterated reconstruction. The study will be performed using wavelet-based shrinkage methods in order to optimize for any problem the suitable threshold method.
Moreover, it is well known that the wavelet approach can be especially successful in detecting the edges of an image. We would like to study the following research themes:
- finding a relation between the classical wavelet methods for deblurring of images (noisy and blurred) and the multigrid methods (see B2-a) for deblurring of images (noisy and blurred);
- incorporate the more precise anti-reflective BCs (see B2-b)) into a wavelets framework in order to improve the precision of the resulting reconstruction.
B3) New applications
a) Nonlinear inverse problems
The two-level "Quasi-Newton" approach previously implemented for Remote Sensing problems, as described in the Scientific Background, has revealed the possibility of using structured matrix tools in the new applicative context of nonlinear inverse problems. In this project we plan to apply the same technique to a non-linear inverse scattering problem arising in electrical engineering, where the dielectric permittivity has to be determined.
b) Web search engines
Here the structure that dominates the problem is indicated by the key words "sparsity" and "block structure": a further rank-one matrix multiplied by 1-c, 0<c<1, is added for modelistic reasons (see [36]). When c is far from one, the computation of the positive PageRank(c) vector (left stationary vector) can be performed quite fastly by the standard power method since the second large eigenvalue has at most modulus c. However, if c is close to 1 then the technique is too slow also due to the huge size of the WEB matrix. We would like to express the positive vector PageRank(c) as a rational expansion of PageRank(1) in such a way that its computation can be obtained eventually through vector extrapolation methods.
c) Identification problems
Non-classical inverse eigenvalue problems for tridiagonal matrices and other closely related matrix structures (also including diagonal-plus-semiseparable matrices) are recently arising as models in identification problems in civil engineering, where localized damages or constitutive coefficients in (multilayer) rods and frames are to be detected from measures of frequency responses or resonances. We are aimed at solving both theoretical and computational issues arising in this framework, e.g., identifiability results, perturbative analysis, and synthesis of numerical methods to identify a properly structured matrix from spectral data.

All the research themes considered here require numerical validations and the development of dedicated prototype software. Hence the group needs to incorporate additional young people, for which a substantial part of funding will be used: temporary research positions, fellowships and occasional jobs are planned to start during the two years of the project.
The group includes researchers of several different universities and a strong cooperation is needed to obtain the expected results, also with the help of external collaborators: M.Benzi, Atlanta-USA (topics A2-a), D.Calvetti and G.Saidel, Ohio-USA (topic B1-c), G.Golub, Stanford-USA (topics A1-a, A2-a, B1-c), A.Morassi, Udine-Italy (topic B3-c), D.Noutsos and P.Vassalos, Ioannina-Greece (topics A1-a, A2-c), M.Pastorino, DIBE Genova-Italy (topic B3-a), M.Tasche, Rostock-Germany (topic B2-c) and people of other groups in this project (point A2-b).
Therefore the funding will cover also several travel expenses due to individual meetings and participation to conferences.

 

 

     Unity of CAGLIARI

     Leader SEBASTIANO SEATZU

 

 

     Task

As agreed upon with the other units, the Cagliari unit will conduct research on the following arguments:
1) development of innovative methods to solve linear systems with ill-conditioned matrices and linear systems with multi-index matrices;
2) development of preconditioners to solve linear systems with multi-index matrices;
3) generation of families of multi-index test matrices, identifying the asymptotic behaviour of their condition numbers;
4) numerical solution of integral equations with multivariate structured kernels and their application to acoustics and remote sensing;
5) development of prototypical software specialized to structured matrices.

1) As to the first research argument, the Cagliari unit proposes to develop computational methods for ill-conditioned structured linear systems by applying fast algorithms based on displacement structure to Tikhonov regularization methods. In the very frequent cases in which the regularization matrix itself is structured, these methods should allow us to deal with linear systems characterized by various types of displacement structure, such as for instance Toeplitz, Hankel and Cauchy structures, which are of particular interest in the numerical solution of integral equations with structured kernels. We also propose to study the solution of ill-conditioned linear systems by means of preconditioned iterative methods applied in Krylov spaces, such as for instance GMRES, the conjugate gradient method and the Lanczos method. More precisely, we are thinking of using the preconditioning methods mentioned in 2) as preconditioners in various iterative methods and in particular in the GMRES method.
A second part of the research will regard the extension of projection methods in weighted Wiener algebras to the solution of linear systems with structured multi-index matrices. This method, which is well-known when solving linear Toeplitz systems in the Wiener algebra, is easily extended to weighted Wiener algebras, provided one deals with mono-index block Toeplitz systems or with positive definite multi-index block Toeplitz systems with scalar blocks [17],[2]. To the contrary, the extension to the multi-index case is not easy for non positive definite Toeplitz systems or for positive definite block Toeplitz systems. Such an extension would be very useful because it would allow us to characterize the asymptotic behaviour of the error as a function of the decay of the elements of the matrix away from the diagonal and of the decay of the elements of the vector on the right-hand side.
The local unit also proposes to study the solution of linear systems whose matrices have a nonconventional structure stemming from problems in astronomy such as the transfer of polarized light in planetary media where there exists considerable local expertise [25].

2) As is well-known [10,11,39], there exist various methods of optimal computational complexity to solve mono-index linear systems. On the other hand, there exist only partial results on the generation of preconditioners of low computational complexity in the case of multi-index matrices [39]. Research in progress by the Cagliari unit appears able to generate in a uniform way preconditioners for mono- and multi-index matrices. The research is presently evolving towards the optimization of the computational complexity independently of the number of indices. Here we propose to experimentally determine the efficiency of the solution method for an extensive class of linear systems by means of preconditioned iterative methods. The development of a vast class of test matrices, as to be discussed in 3), should allow us to validate the method in an appropriate way. We also plan research on the recursive generation of preconditioners that are to be applied iteratively to e.g. the GMRES method. More explicitly, the basic idea is the following: once the preconditioner has been constructed for the matrix at step k, the preconditioner at the next step is obtained from its predecessor by simply evaluating two additional elements.

3) In [28] a general method has been developed that can be applied in a simple way to generate bi-infinite positive definite matrices. In [35] this method has been improved considerably and generalized to the multi-index case, but only to bi-infinite matrices. Here we propose to adapt the method to generate semi-infinite matrices, aiming at identifying the asymptotic behaviour of the condition numbers of an extensive class of multi-index matrices as a function of their order. Such a result should allow us to have at our disposal a sufficiently large class of matrices for which the asymptotic behaviour of their condition numbers is known. The purpose is to be able to compare the efficiency of the newly proposed methods to those in present use. The Cagliari unit is particularly interested in the evaluation of the efficiency of the preconditioners proposed in connection with the conjugate gradient method and the GMRES method.

4) Because linear systems with structured multi-index matrices represent the discrete analogue of integral equations in several variables with structured kernel, there exists a natural and deep connection between the two research topics that is potentially very profitable but has presently not been valued in an adequate way. This depends fundamentally on the fact that operator theory has so far been developed by functional analysts and structured matrix theory by numerical analysts. Because either type of expertise is present in the local unit, we could in fact develop innovative methods for the numerical solution of multivariate integral equations with structured kernel. From the theoretical point of view the computational method for infinite multi-index systems proposed in [34] should be very useful. An important applied field where there exist significant contributions of the local unit [1],[29], is represented by inverse scattering in acoustics and quantum mechanics. In collaboration with scientists at the Universities of Cagliari and Napels that are experts on remote sensing, the local unit proposes to apply the theoretical and numerical results on integral equations and linear multi-index systems to solve typical problems in this field.

5) In spite of a great interest, both from the scientific and the applied side, in the arguments connected with structured matrices and their treatment, there does not seem to exist a specialized software library which implements the most recent algorithms that have specifically been developed for this class of matrices. Because there exist various algorithms and programs of applied interest developed by the local unit as well as a real expertise in writing and evaluating specialized software (*), we propose to develop a toolbox specialized to structured matrices in MatLab, which is one of the most widely disseminated software instruments for scientific computation and visualization. The characteristics to which we will give particular attention are the following: user friendliness, integration with the Matlab environment and its intrinsic numerical routines, and optimization of computing time and memory requirements. As a first goal we intend to pay attention to certain types of structured matrices and basic algorithms for solving linear systems, both by using preconditioned iterative methods and by using direct methods based on displacement structure. The final goal is to construct software that is flexible and can easily be extended to new computational algorithms, also in the multi-index case.
(*) G. Rodriguez has developed the codes for the numerical algorithms proposed in [35-38] and M. Redivo Zaglia is software editor of the journal Numerical Algorithms.