Numerical Linear Algebra

in Signals and Systems

International workshopHotel porto Giardino

Monopoli (Bari)September 11-15, 2006

Patrocinio del Istituto per le Katholieke Patrocinio delPresidente della Applicazioni del Universiteit Presidente dellaGiunta Regionale Calcolo “M.Picone” Leuven Provincia

della Puglia CNR di Bari

Istituto Nazionale di Alta MatematicaGruppo Nazionale per il Calcolo Scientifico

Contents

1 Aims & Topics 4

2 Time Schedule 5

2.1 Short Timetable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Extended Time Schedule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 Contributed Posters 8

4 Sponsors 9

5 Organization 10

6 Special Issue 11

7 Abstracts 12

8 List of participants 56

3

1 Aims & Topics

The cross fertilization between numerical linear algebra and digital signal processing has beenvery fruitful in the last decades. The interaction between them has been growing, leading to manynew algorithms.

Numerical linear algebra tools, such as eigenvalue and singular value decomposition and theirhigher–extension, least squares, total least squares, recursive least squares, regularization, orthog-onality and projections, are the kernels of powerful and numerically robust algorithms.

The goal of this workshop is to bring together researcher in numerical linear algebra anddigital signal processing to discuss new efficient and reliable numerical linear algebra tools forsignal processing applications.

Areas and topics of interest for the workshop include (but are not limited to):

• Singular value and eigenvalue decompositions, including applications.

• Toeplitz, Cauchy, Vandermonde and semiseparable matrices, including special algorithmsand architectures.

• Recursive least squares in digital signal processing.

• Updating and downdating techniques in linear algebra and signal processing.

• Stability and sensitivity analysis of special recursive least squares problems.

Numerical linear algebra in:

• Biomedical signal processing applications.

• Adaptive filters.

• Remote sensing.

• Acoustic echo cancellation.

• Blind signal separation and multiuser detection.

• Multidimensional harmonic retrieval and direction of arrival estimation.

• Applications in wireless communications.

• Applications in pattern analysis and statistical modeling.

• Sensor array processing.

4

2T

ime

Sch

edule

2.1

Short

Tim

eta

ble

Opening and Welcome: Monday from 8.50–9.00.

Monday Tuesday Wednesday Thursday Friday

9.00–9.40 P.C. Hansen E. Tyrtyshnikov P. Dewilde R. Chan Z. Strakos

9.40–10.20 H. Fassbender P. Benner V. Olshevsky L. Reichel H. Voss

10.20–10.50 Coffee break & Poster Session

10.50–11.30 X.-W. Chang I. Markovsky Y. Eidelman O. Menchi D. Fasino

11.30–12.10 J.H. Husøy A. Dax S. Delvaux M. Donatelli V. Sima

16.00–16.40 P. Van Dooren S. Van Huffel Excursion A. Arico N. Mahdavi-Amiri

16.40–17.10 Coffee break & Poster Session

17.10–17.50 K. Hueper D. Sima Excursion S. Kunis Y. Vanberghen

17.50–18.30 J. Duintjer Tebbens L. Dimoccoli & B. Truyen Excursion J. Keiner

5

2.2 Extended Time Schedule

Monday, September 11

Opening and Welcome

• 8.50–9.00: Nicola Mastronardi

Morning Session, Chair: Zdenek Strakos

• 9.00–9.40: Per Christian Hansen, Subspace-Based Noise Reduction via Diagonal and Trian-gular Decompositions.

• 9.40–10.20: Heike Fassbender, A symplectic algorithm for matrix equations arising in theanalysis of stationary Gaussian reciprocal processes.

• 10.50–11.30: Xiao-Wen Chang, Solving Integer Least Squares Problems.

• 11.30–12.10: John Hakon Husøy A unified approach to adaptive filters and their perfor-mances.

Afternoon Session, Chair: Heinrich Voss

• 16.00–16.40:Paul Van Dooren, Optimal correlation between projected matrices.

• 16.40–17.10:Knut Hueper, On the Relationship between FastICA Algorithms and RayleighQuotient Iterations.

• 17.10–17.50:Jurjen Duintjer Tebbens, Implementational Aspects of Linear Discriminant Anal-ysis for Classification Tasks.

Thuesday, September 12

Morning Session, Chair: Paul Van Dooren

• 9.00–9.40: Eugene Tyrtyshnikov, Fast dimensionality reduction and approximate simultane-ous diagonalization of matrices.

• 9.40–10.20: Peter Benner, Model Reduction for Data-Sparse Systems Using HierarchicalMatrices.

• 10.50–11.30:Ivan Markovsky, Structured Low-Rank Approximation in Signal Processing.

• 11.30–12.10:Achiya Dax, Orthogonalization via Deflation :A Minimum Norm Approach forLow Rank Approximations of a Matrix.

Afternoon Session, Chair: Patrick Dewilde

• 16.00–16.40:Sabine Van Huffel, Numerical linear algebra in computational biomedical signalprocessing.

• 16.40–17.10:Diana Sima, Quantification of magnetic resonance spectroscopic signals and sep-arable constrained nonlinear least squares.

• 17.10–17.50: Bart Truyen, A subspace based structure preserving solution method for theElectrical Impedance Tomography problem: Part I.Luca Dimiccoli, A subspace based structure preserving solution method for the ElectricalImpedance Tomography problem: Part II.

6

Wednesday, September 13

Morning Session, Chair: Eugene Tyrtyshnikov

• 9.00–9.40: Patrick De Wilde, Handling Formally Semi-separable Hierarchical System.

• 9.40–10.20: Vadim Olshevsky, to be announced.

• 10.50–11.30: Yuli Eidelman, QR factorization and reduction to Hessenberg form for qua-siseparable matrices.

• 11.30–12.10: Steven Delvaux, Representation and algorithms for unitary-plus-low-rank rankstructured matrices.

Thursday, September 14

Morning Session, Chair: Achiya Dax

• 9.00–9.40: Raymond Chan, Two-Phase Methods for Deblurring Images Corrupted by ImpulsePlus Gaussian Noise.

• 9.40–10.20: Lothar Reichel, Iterative Methods for Inverse Problems.

• 10.50–11.30: Ornella Menchi, Regularizing inverse preconditioners for symmetric band Toeplitzmatrices.

• 11.30–12.10: Marco Donatelli, Regularization by multigrid-type algorithms.

Afternoon Session, Chair: Raymond Chan

• 16.00–16.40: Antonio Arico, The anti-reflective algebra: structural, computational and spec-tral considerations (with application to image deblurring).

• 16.40–17.10: Stefan Kunis, Random sampling of sparse trigonometric polynomials.

• 17.10–17.50: Jens Keiner, Fast evaluation of quadrature formulae on the sphere.

Friday, September 15

Morning Session, Chair: Per Christian Hansen

• 9.00–9.40: Zdenek Strakos, Lanczos tridiagonalization and Golub-Kahan bidiagonalization:Ideas, connections and impact.

• 9.40–10.20: Heinrich Voss, Reducing large gyroscopic eigenproblems by automated multi-levelsubstructuring.

• 10.50–11.30: Dario Fasino, Eigenvalue computations of monodromy operators for linearFDEs with periodic coefficients.

• 11.30–12.10: Vasile Sima, Structure-exploiting algorithms in system identification.

Afternoon Session, Chair: Lothar Reichel

• 16.00–16.40: Nezam Mahdavi-Amiri, Integer ABS class of algorithms for solving linear Dio-phantine equations.

• 16.40–17.10: Yvette Vanberghen, The use of lower semiseparable matrices to solve general-ized eigenvalue problems.

7

3 Contributed Posters

Santa Agreste and Albert Vocaturo, Watermarking Scheme with Multiwavelet Filters.

Rogier Ellenbroek, Distributed numerical linear algebra in control for adaptive optics.

Damiana Lazzaro, Three-Dimensional Dyadic Wavelet Transform for Image Volume Denois-ing.

Ctirad Matonoha, A shifted Steihaug-Toint method for computing a trust-region step.

Serena Papi, A Doubly iterative wavelet algorithm for edge preserving image Reconstruction.

8

4 Sponsors

• Gruppo Nazionale Calcolo Scientifico.http://gruppi.altamatematica.it/gncs/

• INDAM.http://www.altamatematica.it/

• Regione Puglia.http://www.regione.puglia.it/quiregione/start.php

• Katholieke Universiteit Leuven.http://www.kuleuven.be

• The European Association for Signal, Speech and Image Processing.http://www.eurasip.org/content/

• Consiglio Nazionale delle Ricerche.http://www.cnr.it/

• Istituto per le Applicazioni del Calcolo.http://www.iac.rm.cnr.it

• Patrocinio della Provincia di Bari.

9

5 Organization

Organizing committee

Nicola Mastronardi, Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionaledelle Ricerche, Bari, Italy.

Marc Van Barel, Department of Computer Science, Katholieke Universiteit Leuven, Belgium.

Raf Vandebril, Department of Computer Science, Katholieke Universiteit Leuven, Belgium.

Local organizing committee

Diele Fasma, Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionale delleRicerche, Bari, Italy.

Lamura Antonio, Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionaledelle Ricerche, Bari, Italy.

Laudadio Teresa, Istituto di studi sui sistemi intelligenti per l’automazione, Consiglio Nazionaledelle Ricerche, Bari, Italy.

Marangi Carmela, Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionaledelle Ricerche, Bari, Italy.

Nico Giovanni, Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionaledelle Ricerche, Bari, Italy.

Notarnicola Filippo, Istituto per le Applicazioni del Calcolo “M. Picone”, Consiglio Nazionaledelle Ricerche, Bari, Italy.

10

6 Special Issue

A special issue of EURASIP Journal on Applied Signal Processing, will be devoted to NumericalLinear Algebra in Signal Processing Applications. Participants in the workshop can consider tosubmit a manuscript related to the talk presented at the workshop for possible publication in thisspecial issue.

• The manuscripts are refereed, according to the general rules of the Journal.

• EURASIP Journal on Applied Signal Processing.http://www.hindawi.com/journals/asp/

• Numerical Linear Algebra in Signal Processing Applications, with the following call forpapers:http://www.hindawi.com/GetPage.aspx?journal=ASP&page=LASP

11

7 Abstracts

Watermarking Scheme with Multiwavelet Filters

Santa Agreste1

Alberto Vocaturo2

This paper proposes a new watermarking scheme for digital color images, in which a logo [9]watermark is embedded into the multiwavelet coefficients of a Discrete Multiwavelet Transform ofthe image.The diffusion of Internet has increased the phenomenon of the digital piracy, in multimedia objects,like software, image, video, audio and text. A right solution to this problem could be applyinginto the digital object a watermarking to establish the original ownership.In general an watermarking algorithm works on a dual scheme: a) watermark embedding, in whichan inference, called watermark, is added to an original digital object to obtain the watermarkedone; b) watermark detection, in which the existence or not of the watermark is determinated.Our algorithm works with digital color images and, in the first step, it embeddes watermark signalinto Discrete Multiwavelet Transform coefficients.A multiwavelet with multiplicity r is a vector-valued wavelet with r components. The choiceto use multiwavelet functions is not casual for digital color image processing. In fact this kindof images are vector objects with three components. In the case of multiwavelets, unlike thewavelet scalar case, some extra degrees of freedom are allowed, which can be used to constructfunctions with several good properties, combining, for example, orthogonality, symmetry, shortsupport and vanishing moments. Therefore multiwavelet, thanks to their good properties, seem tobe more powerful for color image processing applications than scalar wavelet functions [3],[6],[5].Nevertheless, the application of the multiwavelet decomposition scheme requires an additionalstep with respect to the scalar case. It consists in finding r sequences of initial coefficients, neededby the decomposition scheme, from a given set of data. This step is called prefiltering of data,since it can be seen as the application of a filter to the initial data. Prefiltering is crucial in amultiwavelet setting: a bad filter, which does not take into account the properties of the assignedbasis, can give raise to bad results in applications. The prefiltering step can be avoided if balancedmultiwavelets are used.In this work we describe a watermarking scheme based on decomposition and reconstruction ofthe color digital image by means multiwavelet full rank filters with multiplicity r = 3 [4], that arebalanced multiwavelets.The performance of this multiwavelet-based watermarking scheme is compared with the wavelet-based watermarking algorithm [1].

References

[1] S. Agreste, N. Castorina, S. Giovinazzo, D. Prestipino, L. Puccio,Tutela del diritto di pro-prieta’delle immagini digitali: Implementazione di un algoritmo di Watermark mediante fun-zioni Wavelet, Atti della Accademia Peloritana dei Peric., Classe di Scienze Fis., Mat. eNat. (on line). vol. LXXXI-LXXXII, 2005, pp. 1-15, ISSN: 1825-1242. Identification Number:C1A0401009 ID Code: 261.

[2] M. Barni, F. Bartolini, A. Piva, Improved wavelet-based watermarking through pixel wise-masking, Image Processing, IEEE Transactions, Volume 10, Issue 5, May 2001 pp.783 - 791.

1Department of Information Science, Milan University, Italy.agreste@dsi.unimi.it

2Department of Mathematics, University of Messina, Italy.avocat@dipmat.unime.it

12

[3] M. Cotronei, L.B. Montefusco, and L. Puccio, Multiwavelet analysis and signal processing,IEEE Trans. on Circuits and Systems II: Analog and Digital Signal Processing 45, 970-987(1998).

[4] M.Cotronei, L.Puccio and A. Vocaturo,On the application of full rank filters to color imageprocessings, EUA4X#8, TCN CAE 2005, MASCOT05 Proceedings, Lecce to appear.

[5] M. Cotronei, D. Lazzaro, L.B. Montefusco, and L. Puccio,Image compression through embed-ded multiwavelet transform coding, IEEE Trans. on Image Process vol.9 n.2, 184-189 (2000).

[6] M.B.Martin and A. E.Bell,New image compression techniques using multiwavelets and mul-tiwavelets packets, IEEE Trans. on Image Process vol.10 n.4, 500-510 (2001).

[7] F.Keinert,Wavelets and Multiwavelets, Studies in advanced mathematics vol.42, ChapmanHall/CRC Press, Boca Raton, FL (2003)

[8] I. Pitas,A method for signature casting on digital images, in Proc. of the International Con-ference on Image Processing, v. 3, 1996, pp. 215-218.

[9] G. Xie, H. Shen,A New Fusion Based Blind Logo-Watermarking Algorithm, Japan Societyfor the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (B) under GrantNo.14380139

13

The anti-reflective algebra: structural,computational and spectral considerations (with

application to image deblurring)

Antonio Arico1

The choice of suitable boundary conditions (BC) is very important in denoising problems. Anew kind of BC, the so-called “Anti-Reflective” (AR) boundary conditions, have been recentlyintroduced (in connection with fast deblurring algorithms) both in the case of signals and images;later the AR-BC have been extended to arbitrary dimensions (d = 3, 4, . . . ).The talk will be about some structural, computational and spectral properties of AR-BC: underthe assumption of strong symmetry of the PSF many interesting results can be shown; someexperiments will be presented in the case of image deblurring (d = 2).

References

[1] S. Serra Capizzano 2003 A note on anti-reflective boundary conditions and fast deblurringmodels, SIAM J. Sci. Comput. 25-3 pp. 1307–1325.

[2] A. Arico, M. Donatelli, and S. Serra-Capizzano 2006 The anti-reflective algebra: structuraland computational analysis with application to image deblurring and denoising.

[3] A. Arico, M. Donatelli, and S. Serra-Capizzano 2006 Spectral analysis of the anti-reflectivealgebras and applications.

1Dipartimento di Fisica e Matematica, Universita dell’Insubria, Como, Italy.antonio.arico@unipv.it

14

Model Reduction for Data-Sparse Systems UsingHierarchical Matrices

Peter Benner1

Ulrike Baur2

Model reduction is an ubiquitous tool to facilitate or even enable the simulation, optimiza-tion and control of large-scale dynamical systems. Application areas range from microelectronics(device and circuit simulation) over computational biology, control of mechanical and electricalsystems to flow control and PDE-constrained optimization. In this talk, we will focus on problemsfrom the latter application areas. In particular, we will discuss model reduction techniques basedon system balancing for (optimal) control of parabolic partial differential equations. After spatialsemi-discretization by FEM or BEM methods, large-scale, sparse (in case of FEM) or data-sparse(in case of BEM) linear control systems are obtained. Due to the cubic complexity of standardimplementations of balanced truncation and relatives, it is not possible to apply these methodsdirectly to such systems.

1Department of Mathematics, Chemnitz University of Technology, Chemnitz, Germany.benner@mathematik.tu-chemnitz.de http://www-user.tu-chemnitz.de/∼benner/

2Institut fur Mathematik, TU Berlin,Berlin, Germany.baur@math.tu-berlin.de

15

Two-Phase Methods for Deblurring ImagesCorrupted by Impulse Plus Gaussian Noise 4

Raymond Chan1

Jian-Feng Cai2

Mila Nikolova 3

We consider the problem of restoring blurred images which are corrupted by impulse noise andpossibly Gaussian noise. We solve the problem in two phases. In the first phase, median-typefilters are applied to locate the pixels which are likely to be corrupted by impulse noise. In thesecond phase, we solve an inverse problem from incomplete data by a variational method. Ourscheme can restore blurred images corrupted by 90 % salt-and-pepper noise or 50 % random-valuedimpulse noise.

1Department of Mathematics, The Chinese University of Hong Kong, Hong Kong, China.rchan@math.cuhk.edu.hk http://www.math.cuhk.edu.hk/∼rchan

2Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong.cwho@math.cuhk.edu.hk

3Centre de Mathematiques et de Leus Applications, ENS de Cachan, Cachan Cedex.nikolova@cmla.ens-cachan.fr

4This work was supported by HKRGC Grant CUHK 400503.

16

Solving Integer Least Squares Problems

Xiao-Wen Chang1

Integer least squares (ILS) problems may arise from many applications such as communications,cryptograph, global positioning systems, etc, see, e.g., [1], [6], and [7]. Unlike real LS problems,ILS problems are NP-hard. Solving an ILS problem usually consists of two phases: reduction andsearch. The reduction process transforms the original problem to a new problem to make thesearch process more efficient and we can cast its crucial part as a matrix factorization. The searchprocess is to generate and examine the integer points over a specific region and find the optimalsolution.

In this talk, we first briefly introduce a typical algorithm for ordinary ILS problems andthen present some of our recent work on various ILS problems, which may subject to differentconstraints, see [2], [3], [4] and [5]. Some simulation results will be given to show our algorithmscan be significantly faster than the current ones in the literature.

References

[1] E. Agrell, T. Eriksson, A. Vardy, and K. Zeger, Closest point search in lattices, IEEETrans. Inform. Theory, 48 (2002), pp. 2201–2214.

[2] X.-W. Chang, G. Golub, and Q. Han, Solving elliposoid-constrained integer least squaresproblems. Submitted, 2006.

[3] X.-W. Chang and Q. Han, Solving box-constrained integer least squares problems. Sub-mitted, 2005.

[4] X.-W. Chang and X. Yang, A new fast generalized sphere decoding algorithm for under-determined MIMO systems. In Proceedings of 23rd Queen’s Biennial Symposium on Com-munications, Kingston, Ontario, Canada, May 30 - June 1, 2006.

[5] X.-W. Chang, X. Yang, and T. Zhou, MLAMBDA: A modified LAMBDA method forinteger least-squares estimation, J. of Geodesy, 79 (2005), pp. 525–565.

[6] M. O. Damen, H. El Gamal, and G. Caire, On maximum-likelihood detection and thesearch for the closest lattice point, IEEE Trans. Inform. Theory, 49 (2003), pp. 2389–2402.

[7] P. Teunissen, The least-squares ambiguity decorrelation adjustment: a method for fast GPSambiguity estitmation, J. of Geodesy, 70 (1995), pp. 65–82.

1School of Computer Science, McGill University, Montreal, Canada.chang@cs.mcgill.ca http://www.cs.mcgill.ca/∼chang/index.php

17

Orthogonalization via Deflation :A MinimumNorm Approach for Low Rank Approximations of

a Matrix

Achiya Dax1

In this talk we introduce a new orthogonalization method. Given a real m × n matrix A,the new method constructs an SVD-type decomposition of the form A = UΣV T . The columnsof U and V are orthogonal, or nearly orthogonal, while Σ is a diagonal matrix whose diagonalentries approximate the singular values of A. The method has three versions. A ”left-side” or-thogonalization scheme in which the columns of U constitute an orthonormal basis of Range(A).A ”right-side” orthogonalization scheme in which the columns of V constitute an orthonormalbasis of Range(AT ). In the third version both U and V has orthonormal columns, but the thedecomposition is not exact.

Starting from A1 = A the deflation process generates a sequence of matrices A1, A2, A3, . . . ,by the rule

Ak+1 = Ak − σkukvTk , k = 1, 2, 3, . . . ,

where σk, uk, vk denotes a computed estimate of a singular triplet of Ak. The estimated singularvectors are obtained by a few ”rectangular iterations” for solving the minimum norm problem

minimizeF (u, v) = ‖Ak − uvT ‖2F .

The resulting orthogonal decomposition may substitute the SVD in many applications. Theadvantage of the new method lies in problems with missing data. That is, when some entries ofA are unknown. Standard SVD algorithms are unable to handle such matrices. Yet the minimumnorm approach overcomes this difficulty in an elegant way: The objective function is redefined as

F (u, v) =∑

(aij − uivj)2,

where the sum is restricted to known entries of A. This modification enables us to construct alow-rank approximation of A, or a ”pseudo SVD” of A, in spite of the missing data. Once a pseudoSVD of A is constructed, it can be used for estimating the missing data. Numerical experimentsillustrate the usefulness of the proposed methods.

1Hydrological Service, P.O.B. 36118, Jerusalem 91360, ISRAEL.dax20@water.gov.il

18

Handling Formally Semi-separable HierarchicalSystems

Patrick De Wilde1

Many systems of equations can be put in a form that can be characterized as ’formally semi-separable’, by which is meant that their formal structure looks semi-separable while they are notactually semi-separable in the classical sense because of the non- elementary nature of the compo-nents in the semi-separable representation. An example is a tri-diagonal system of block matrices,in which the constitutive blocks are themselves tri- diagonal matrices (or even tri-diagonal blockmatrices, leading to a deeper hierarchy - the standard case for 3D partial differential equations).The usual way of handling such systems is by finding a new order that is more efficient for systemsolving than the original, e.g. nested bisection or a derivative of it. This strategy may work well insimple cases, but breaks down for larger systems. It also fails to exploit the system structure. Wepropose new strategies that do exploit the formal structure, but have to resort to model reductionsteps in order to keep computational complexity linear in the size of the defining data - in contrastto the class of hierarchical semi-separable systems, the class of formally semi- separable hierarchi-cal systems is not closed under system inversion. We also show for some characteristic examplesthat the model reduction method presented leads to numerically accurate results. In fact, we areable to solve some important cases such as occur in global motion estimation explicitly.

1Faculty of Electrical Engineering, Technische Universiteit Delft, Delft, The Netherlands.dewilde@cas.et.tudelft.nl http://cobalt.et.tudelft.nl/people/dewilde.html

19

Representation and algorithms forunitary-plus-low-rank rank structured matrices

Steven Delvaux1

Marc Van Barel2

In this talk we consider unitary/Givens representations for unitary-plus-low-rank rank struc-tured matrices. We design several matrix algorithms using these representations, e.g., the Hessen-berg reduction, implicit QR-algorithm, . . .

1Department of Computer Science, Katholieke Universiteit Leuven, Belgium.Steven.Delvaux@cs.kuleuven.be http://www.cs.kuleuven.be/∼steven

2Department of Computer Science, Katholieke Universiteit Leuven, Belgium.Marc.VanBarel@cs.kuleuven.be http://www.cs.kuleuven.be/∼marc

20

A subspace based structure preserving solutionmethod for the Electrical Impedance Tomography

problem: Part II

Luca Dimiccoli1

Bart Truyen2

Jan Cornelis3

We present a subspace based, structure preserving solution method for the problem of ElectricalImpedance Tomography, where the conductivity inside a simply connected 2-dimensional domainis sought from noisy and incomplete boundary data. Unlike conventional output-least squaresalgorithms that can be regarded as minimizing a certain error norm, solutions are recovered hereas the minimizers of a closely related residual norm problem. An iterative solution scheme is shownto lead to a sequence of sparse matrix subproblems, with conditioning far more favourable thantypically observed in output-least squares. We find that these sparse subproblems demonstrate aparticular form of displacement structure that can be further elaborated to finally arrive upon anefficient computational implementation. In the second part of this contribution, we investigate inmore detail the different sparsity patterns deriving from the problem formulation, and illustratehow these connect with the original concept of displacement structure.

1ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUBldimicco@etro.vub.ac.be http://www.etro.vub.ac.be/Members/Personal Page.asp?PM ID=1463

2ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUBbatruyen@etro.vub.ac.be http://www.etro.vub.ac.be/Personal/btruyen/

3ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUB.

21

Regularization by multigrid-type algorithms

Marco Donatelli1

Stefano Serra-Capizzano2

We consider the de-blurring problem of noisy and blurred signals/images in the case of space in-variant point spread functions. The use of appropriate boundary conditions leads to linear systemswith structured coefficient matrices related to space invariant operators like Toeplitz, circulants,trigonometric matrix algebras etc. We combine an algebraic multigrid (which is designed ad hocfor structured matrices) with the low-pass projectors typical of the classical geometrical multigridemployed in a PDEs context. Thus, using an appropriate smoother, we obtain an iterative reg-ularizing method (see [2, 1]) based on: projection in a subspace where it is easier to distinguishbetween the signal and the noise and then application of an iterative regularizing method in theprojected subspace. Therefore any iterative regularizing method like conjugate gradient (CG),conjugate gradient for normal equation (CGNE), Landweber etc., can be used as smoother in ourmultigrid algorithm. The projector is chosen in order to maintain the same algebraic structure ateach recursion level and having a low-pass filter property, which is very useful in order to reducethe noise effects. In this way, we obtain a better restored image with a flatter restoration errorcurve and also in less time than the auxiliary method used as smoother. As direct consequencethe choice of the exact iteration where to stop is less critical than in other regularizing iterativemethods. Furthermore, we can choose multigrid procedures which are extremely more efficientthan classical techniques without losing accuracy in the restored image.

References

[1] M. Donatelli. Image Deconvolution and Multigrid Methods, PhD Thesis in Applied andComputational Mathematics, University of Milano, March 2006.

[2] M. Donatelli and S. Serra Capizzano, On the regularizing power of multigrid-typealgorithms, SIAM J. Sci. Comput., 27–6 (2006), pp. 2053–2076.

1Dipartimento di Fisica e Matematica, Universita dell’Insubria, Como, Italy.marco.donatelli@uninsubria.it

2Dipartimento di Fisica e Matematica, Universita dell’Insubria - Sede di Como, Italy.s.serracapizzano@uninsubria.it

22

Implementational Aspects of Linear DiscriminantAnalysis for Classification Tasks

Jurjen Duintjer Tebbens1

Pavel Schlesinger2

In applications such as pattern analysis, face recognition, text document classification or proteinfold prediction, one often faces the problem to classify with a number of characterizing variableslargely exceeding the number of given samples. When applying one of the most popular classifi-cation methods, Fisher’s Linear Discriminant Analysis (FLDA) [2], this leads to large generalizedeigenproblems with singular matrices.

The talk addresses a natural modification of FLDA to circumvent the singularity of the matrixpencil [1, 3] and focusses on its efficient implementation. It considers, among others, tools fromPrincipal Component Analysis as well as new techniques to optimize implementation [1], appro-priate choices of numerical methods for the individual steps of the algorithm and exploitation ofpossible sparsity.

This work was supported by the Program Information Society under project 1ET400300415(first author) and the MSMT CR Project LC536 (second author).

References

[1] J. Duintjer Tebbens and P. Schlesinger: Improving Implementation of Linear DiscriminantAnalysis for the Small Sample Size Problem. submitted to Computational Statistics and DataAnalysis (special issue), in March 2006.

[2] R.-A. Fisher: The use of multiple measurements in taxonomic problems. Ann. Eugenics, vol.7, 178–188, 1936.

[3] J. Yang and J.-Y. Yang: Why can LDA be performed in PCA transformed space? PatternRecognition, vol. 36(2), 563–566, 2003.

1Institute of Computer Science, Academy of Sciences of the Czech Republic.tebbens@cs.cas.cz

2Institute of Formal and Applied Linguistics, Charles University, Prague.schlesinger@ufal.ms.mff.cuni.cz

23

QR factorization and reduction to Hessenbergform for quasiseparable matrices

Yuli Eidelman1

Luca Gemignani2

Israel Gohberg3

We discuss a fast algorithm for reduction of a quasiseparable matrix to a Hessenberg form. Theimportant part of the method is the factorization of a given quasiseparable matrix as a product ofa unitary block lower triangular matrix and a block upper triangular matrix. This factorizationis also an essential part of the QR factorization algorithm for quasiseparable matrices.

1Department of Mathematics, Tel-Aviv University, Israel.eideyu@post.tau.ac.il

2Departimento di Matematica, Universita di Pisa, Italy.gemignan@dm.unipi.it

3Department of Mathematics, Tel-Aviv University, Israel.gohberg@post.tau.ac.il

24

Distributed numerical linear algebra in controlfor adaptive optics

Rogier Ellenbroek1

Michel Verhaegen2

Niek Doelman3

Adaptive Optics systems used in modern ground-based telescopes use an optical corrector such asa Deformable Mirror (DM) to reduce wavefront aberrations due to atmospheric turbulence. Thewavefront aberrations can be measured and are used in a feedback control system. However, thecontrol algorithms consist of a number of steps that are computationally demanding and scale uprapidly with the diameter of the telescope’s primary mirror. For AO systems of future large tele-scopes, this computational complexity rises far beyond the capabilities of single processor systems,which raises the need for efficiently parallelizable control algorithms. Based on the structure ofthe problem, a distributed approach towards controller design is being taken [1]. It is distributedin the sense that each actuator (or set of neighboring actuators) has its own local controller thatcan communicate only directly with its neighbors.One particularly demanding step in the control algorithm, is the so called wavefront reconstructionstep, in which a 2D surface has to be reconstructed from locally measured spatial gradients. Thisyields a large but sparse system of linear equations, that can be efficiently solved using iterativemethods. Numerical complexity that grows linearly with the number of actuators can be achieved[2] using a preconditioned conjugate gradient method, but this method does not fit into the pro-posed distributed processor topology. This is possible for a successive over-relaxation method [1],but this does not have the order-N efficiency.A possible solution is to exploit the spatio-temporal structure in the wavefront aberrations. Insimulations this shows promise, but it represents a possibly restrictive assumption on the allowedaberrations. A distributed reconstruction algorithm is now sought that does not have this restric-tiveness, while it is scalable in the sense that for a fixed accuracy, the required computationalpower of each distributed processor is independent of the number of actuators of the DM.

References

[1] Rogier Ellenbroek, Michel Verhaegen, Roger Hamelinck, Niek Doelman, Maarten Steinbuch,and Nick Rosielle. Distributed control in adaptive optics - deformable mirror and turbulencemodeling. In Proceedings of the SPIE conference on astronomical telescopes and instrumenta-tion. SPIE, May 2006.

[2] Luc Gilles. Order-n sparse minimum-variance open-loop reconstructor for extreme adaptiveoptics. Optics Letters, 28(20):1927–1929, October 2003.

[3] Roger Hamelinck, Nick Rosielle, Maarten Steinbuch, Rogier Ellenbroek, Michel Verhaegen,and Niek Doelman. Actuator tests for a large deformable membrane mirror. In Proceedingsof the SPIE conference on astronomical telescopes and instrumentation - advances in adaptiveoptics. SPIE, May 2006.

1Delft University of Technology, Delft, The Netherlands.r.m.l.ellenbroek@dcsc.tudelft.nl http://www.dcsc.tudelft.nl/∼rellenbroek

2Delft University of Technology, Delft, The Netherlands.M.Verhaegen@dcsc.tudelft.nl http://www.dcsc.tudelft.nl/∼mverhaegen/

3TNO Science and Industry, Opto-mechanical Instrumentation, Stieltjesweg 1, 2628 CK Delft, The Netherlands.niek.doelman@tno.nl

25

Eigenvalue computations of monodromyoperators for linear FDEs with periodic

coefficients

Dario Fasino1

Functional-delayed differential equations are the most natural way to model a variety of dy-namical systems whose status depends also on its past values.

Here, we consider the problem of assessing asymptotic stability of dynamical systems governedby a delayed differential equation,

x′(t) = F (t, xt) 0 ≤ t < +∞x(t) = f(t) −τ ≤ t ≤ 0,

where the right-hand side F (t, xt) is τ -periodic in t and linear in xt(s) = x(t + s), −τ ≤ s ≤ 0.Any such equation admits a monodromy operator S such that, for any integer k, one has

x(kτ + t) = (Skf)(t), −τ ≤ t ≤ 0.

Thus, it is the spectral decomposition of this operator (in particular, its eigenvalues) that primarilymatters in understanding the long-term behavior of those systems.

The purpose of this talk is to lay down an operator-theoretic framework for the numericalcomputation of eigenvalues and eigenfunctions of S, reducing the problem to standard numericallinear algebra computations. We also discuss the way structured matrix technology enters intothe picture.

1Dipartimento di Matematica e Informatica Universit degli Studi di Udine, Italy.fasino@dimi.uniud.it http://www.dimi.uniud.it/fasino/

26

A symplectic algorithm for matrix equationsarising in the analysis of stationary Gaussian

reciprocal processes

Heike Fassbender1

Peter Benner2

The nonlinear matrix equation

X = f(X) with f(X) = Q+ LX−1LT , (1)

where Q = QT ∈ Rn×n is positive definite and L ∈ R

n×n is nonsingular, arises in the analysis ofstationary Gaussian reciprocal processes over a finite interval; see, e.g. [2]. The problem is to finda positive definite symmetric solution X of (1).

In [1] it is noted that if X solves (1), then it also obeys the iterated equation

X = f(f(X)) = Q+ F (R−1 +X−1)−1FT

with F = LL−T and R = LTQ−1L = RT positive definite. Using the Sherman-Morrison-Woodbury formula to derive an expression for (R−1 +X−1)−1, we obtain

DR(X) = Q+ FXF T − FX(X +R)−1XFT −X,

a discrete-time algebraic Riccati equation. Because (F, I) is observable and (F,Q1

2 ) is controllable,a unique positive definite solution X? exists. This unique solution coincides with that solution of(1) one is interested in.

In this talk, we propose to compute an approximate solution of the DARE by the (butterfly)SZ algorithm applied to the corresponding symplectic pencil where zero and infinity eigenvaluesare removed using an iterative deflation strategy. This algorithm is a fast, reliable and structure-preserving algorithm for computing the stable deflating subspace of the symplectic matrix pencilassociated with the DARE. We compare our algorithm to existing iterative methods to solve (1).

References

[1] A. Ferrante and B.B. Levy. Hermitian Solutions of the Equation X = Q+NX−1N?. LinearAlgebra and its Applications 247, pp. 359-373, 1996.

[2] B.C. Levy, R. Frezza and A.J. Kerner. Modeling and Estimation of Discrete-Time GaussianReciprocal Processes. IEEE Trans. Automat. Control AC-90, pp. 1013-1023, 1990.

1Institut Computational Mathematics, TU Braunschweig, Germany.h.fassbender@tu-bs.de http://www-public.tu-bs.de:8080/∼hfassben/index en.html

2Technische Universitat Chemnitz, Fakultat fur Mathematik, Chemnitz, Germany.benner@mathematik.tu-chemnitz.de

27

Subspace-Based Noise Reduction via Diagonaland Triangular Decompositions

Per Christian Hansen1

Søren Holdt Jensen2

Rank reduction is a general principle for finding the right trade-off between model bias andmodel variance when reconstructing signals from noisy data. – L. L. Scharf and D. W. Tufts, 1987.

The signal subspace approach has proved itself useful for noise reduction, speech enhancementand many other applications in signal processing. The area has grown dramatically over thelast 20 years, along with advances in efficient computational algorithms for matrix computationg,especially singular value decompositions and rank-revealing decompositions.

We survey the underlying mathematics and present the techniques and algorithms in a commonframework. In addition to methods based on diagonal decompositions (eigenvalues and singularvalues), we survey the use of rank- revealing triangular decompositions. We also discuss alterna-tives to the classical least-squares formulation, and we show how signals with general (non-white)noise is treated by prewhitening.

Finally we show how all the above algorithms can be interpreted in terms of FIR filters definedfrom the decompositions involved, and we introduce a new analysis tool called “canonical filters”that allows us to compare the performance of different subspace algorithms.

1Informatics and Mathematical Modelling, Technical University of Denmark.pch@imm.dtu.dk http://www.imm.dtu.dk/∼pch

2Information and Signals Division, Department of Communication Technology, Aalborg University, FredrikBajers Vej 7A-3, DK-9220 Aalborg, Denmark.

28

On the Relationship between FastICAAlgorithms and Rayleigh Quotient Iterations

Knut Hueper1

Hao Shen2

In this talk we present recent results in the analysis of FastICA algorithms. FastICA is apopular algorithm for linear Independent Component Analysis (ICA), which is a powerful tool forsolving the problem of Blind Source Separation (BSS) in statistical signal processing [1]. Our anal-ysis suggests a close connection between FastICA algorithms and the Rayleigh Quotient Iteration(RQI) method, which is well known to the numerical linear algebra community.

Originally, RQI was developed as a scalar-shifted version of a simpler algorithm, namely theinverse iteration method [2]. It ensures a local cubic convergence to one eigenvector of a realsymmetric matrix. To analyse RQI rigorously, a proper setting is to do the analysis on the realprojective space RP

m−1 rather than on the unit sphere Sm−1. On RPm−1, RQI is free of a sign

flipping behavior, which the RQI on Sm−1 does always suffer from. More important is thatsmoothness of RQI iterations will hold on RP

m−1, but not on Sm−1. As an aside, RQI can bealso developed as a special case of a Newton type method on Sm−1, specifically the so-calledNewton-Rayleigh quotient method [3].

Analogously, our work shows that FastICA is in fact a scalar-shifted version of a simpler ICAalgorithm developed in [4]. Although the original FastICA was commonly considered as a selfmapon Sm−1, RP

m−1 is actually a more proper setting for analysing FastICA than Sm−1 [5], for thesame reasons as RQI (free of the sign flipping phenomena and smoothness of iterations). Notsurprisingly, the original FastICA is also shown as a special case of a Newton type method, the so-called Approximate Newton Parametric ICA algorithm (ANPICA) [6]. By using techniques similarto the analysis of RQI [7], FastICA is proven to be locally quadratically convergent to a fixed pointcorresponding to correct source separations. Moreover, some generalisations of ANPICA in theframework of a scalar shift strategy are studied in this talk as well.

Finally, due to the discovered similarity between FastICA and RQI, we parallelise ANPICA tothe case where all units converge simultaneously, i.e., we consider some kind of parallel algorithmsliving on the orthogonal group [8]. Several parallel algorithms are proposed using the similarparallelising strategies as for RQI, see [7, 9, 10, 11]. The local quadratic convergence properties ofthese parallel algorithms are discussed as well.

References

[1] P. Comon, “Independent component analysis, a new concept?” Signal Processing, vol. 36,no. 3, pp. 287–314, 1994.

[2] B. N. Parlett, The Symmetric Eigenvalue Problem, vol. 20, SIAM Classics In AppliedMathematics Series, Philadelphia, 1998.

[3] S. T. Smith, “Optimization techniques on Riemannian manifolds,” in Hamiltonian andGradient Flows, Algorithms and Control, A. Bloch, Ed., Fields Institute Communications,pp. 113–136. American Math. Soc., Providence, 1994.

[4] P. Regalia and E. Kofidis, “Monotonic convergence of fixed-point algorithms for ICA,” IEEETransactions on Neural Networks, vol. 14, no. 4, pp. 943–949, 2003.

[5] H. Shen and K. Huper, “Local convergence analysis of FastICA,” in Lecture Notes in Com-puter Science, Proceedings of the 6th International Conference on Independent Component

1National ICT Australia Ltd.knut.hueper@nicta.com.au http://users.rsise.anu.edu.au/∼hueper/

2National ICT Australia Ltd. and RSISE, The Australian National University.hao.shen@rsise.anu.edu.au

29

Analysis and Blind Source Separation (ICA 2006), J. Rosca, D. Erdogmus, J. C. Principe,and S. Haykin, Eds., Berlin/Heidelberg, 2006, vol. 3889, pp. 893–900, Springer-Verlag.

[6] H. Shen, K. Huper, and A.-K. Seghouane, “Geometric optimisation and FastICA algorithms,”To appear in MTNS 2006, Kyoto, Japan, July 24-28, 2006.

[7] K. Huper, A Calculus Approach to Matrix Eigenvalue Algorithms, Habilitation Dissertation,Department of Mathematics, University of Wurzburg, Germany, July 2002.

[8] K. Huper, H. Shen, and A.-K. Seghouane, “Local convergence properties of FastICA andsome generalisations,” in Proceedings of the 31st IEEE International Conference on Acoustics,Speech, and Signal Processing (ICASSP 2006), Toulouse, France, 2006, pp. V1009–V1012.

[9] P.-A. Absil, R. Mahony, R. Sepulchre, and P. Van Dooren, “A Grassmann–Rayleigh quotientiteration for computing invariant subspaces,” SIAM Review, vol. 44, no. 1, pp. 57–73, 2002.

[10] M. Nikpour, Reduced Rank Signal Processing, Ph.D. thesis, Department of Electric andElectronic Engineering, University of Melbourne, Australia, July 2002.

[11] M. Nikpour, K. Huper, and J.H. Manton, “Generalizations of the Rayleigh quotient iterationfor the iterative refinement of the eigenvectors of real symmetric matrices,” in Proceedings ofthe 30th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP2005), Philadelphia, PA, USA, 2005, pp. V1041–V1044.

30

A unified approach to adaptive filters and theirperformances

John Hakon Husøy1

All adaptive filter algorithms share the same common goals, the most important being therapid convergence to a good approximation of the solution of the Wiener-Hopf equation in astationary environment. In spite of the commonality in goals, the theory of adaptive filters ischaracterized by a multitude of algorithms whose derivations rely on a large number ideas thatare often perceived as somewhat unrelated. Approaching the subject of adaptive filtering throughthe elementary theory of stationary iterative linear equation solvers, we present a unified theoryfor adaptive filter algorithms and their performances.

More specifically, our starting point is the Wiener-Hopf equation Rxxht = rxd, where ht isthe M × 1 vector of filter coefficients constituting what we refer to as the true Wiener solution,Rxx is the autocorrelation matrix of the filter input signal and rxd is a cross correlation vector.An iterative solution scheme is found using a preconditioned – C being the preconditioner –Richardson iteration:

h(k + 1) = h(k) + µC[rxd −Rxxh(k)]. (2)

Substituting estimates based on observed signals up to time kN – denote them Rxx(kN), rxd(kN),and C(kN) – for the quantities involved in Eq. 2, and accepting the premise that the mostgeneral estimates for the autocorrelation matrix and the cross correlation vector are Rxx(kN) =X(kN)FFTXT (kN), and rxd(kN) = X(kN)FFT d(kN), where F is a matrix whose columns areanalysis filter unit pulse responses of a K channel orthogonal perfect reconstruction criticallysampled filter bank system and X(kN) is a matrix derived from the input signal of the adaptivefilter, we obtain one of two forms of what we call a generic adaptive filter:

h(k + 1) = h(k) + µC(kN)XF (kN)eF (kN), (3)

where XF (kN) = X(kN)F, dF (kN) = FT d(kN), eF (kN) = FT e(kN), e(kN) = d(kN)−XT (kN)h(k).Based on this, we first show that all major adaptive filter algorithms can be derived from one oftwo forms of our generic adaptive filter through the specification of F, the dimensions of the ma-trix/vector quantities involved in forming Rxx(kN), rxd(kN), and the selection of one out of threeintuitively plausible choices for the structure of C(kN). Finally, following the approach of [1], wederive compact and explicit general performance results that are applicable to all major adaptivefilter families.

References

[1] H.-C. Shin and A. H. Sayed, “Mean-square performance of a family of affine projection algo-rithms,” IEEE Trans. Signal Processing, vol. 52, pp. 90–102, Jan. 2004.

1University of Stavanger.john.h.husoy@uis.no

31

Fast evaluation of quadrature formulae on thesphere

Jens Keiner1

Daniel Potts2

Recently, a fast approximate algorithm for the evaluation of expansions in terms of standardL2(

S2)

-orthonormal spherical harmonics at arbitrary nodes on the sphere S2 was proposed in [3].

The aim of this paper is to develop a fast algorithm for the adjoint problem which on the contrarycan be used to compute expansion coefficients from sampled data by means of quadrature rules.

We give a formulation in matrix-vector notation and an explicit factorisation of the sphericalFourier matrix based on the latter algorithm. Starting from this, we obtain the correspondingfactorisation of the adjoint spherical Fourier matrix and are able to describe the algorithm for theadjoint transformation which can be employed to evaluate quadrature rules for arbitrary nodesand weights. We provide results of numerical tests showing the stability of the obtained algorithmusing as examples classical Gauß-Legendre and Clenshaw-Curtis quadrature rules as well as theHEALPix pixelation scheme in [2] and an equidistribution from [1].

References

[1] W. Freeden, T. Gervens, and M. Schreiner. Constructive Approximation on the Sphere. OxfordUniversity Press, Oxford, 1998.

[2] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, andM. Bartelmann. Healpix: A framework for high-resolution discretization and fast analysis ofdata distributed on the sphere. The Astrophysical Journal, 622:759?–771, 2005.

[3] S. Kunis and D. Potts. Fast spherical Fourier algorithms. J. Comput. Appl. Math., 161:75 –98, 2003.

1Institute of Mathematics, University of Lubeck, Germany.keiner@math.uni-luebeck.de http://www.math.uni-luebeck.de/keiner

2Chemnitz University of Technology, Department of Mathematics, Chemnitz, Germany.potts@mathematik.tu-chemnitz.de,

32

Random sampling of sparse trigonometricpolynomials

Stefan Kunis1

Holger Rauhut2

Recently, the surprising fact that it is possible to recover functions having only few non-zero coef-ficients with respect to some basis from vastly incomplete information has gained much attention.Such functions are commonly called sparse or compressible and they naturally appear in a widerange of applications. We study sparse trigonometric polynomials

f (x) =∑

k∈IN

fke−2πikx, IN := −

N

2, . . . ,

N

2− 1,

with non-zero Fourier coefficients fk ∈ C only on a set Ω ⊂ IN with size |Ω| ¿ N . However,a priori nothing is known about Ω apart from a maximum size. Our aim is to sample f at Mrandomly chosen nodes xj ∈ [− 1

2, 12] and try to reconstruct f from these samples.

Due to a recent result, the greedy algorithm Orthogonal Matching Pursuit succeeds in this taskwith high probability by doing |Ω| iterations if the number of samples fulfils M = O(|Ω| logN).The proposed scheme adds in each iteration one ’likely’ index k ∈ IN to its active set Ω′ andprojects the right hand of given samples (f(xj))j=0,...,M−1 onto the columns of

AX,Ω′ :=(

e−2πikxj)

j=0,...,M ;k∈Ω′

In this talk, we focus on the computational complexity of Orthogonal Matching Pursuit withinthe present setting. We use the nonequispaced FFT and particular updating techniques for theprojection to obtain a fast total scheme. We illustrate our observations with numerical experimentsthat indeed show that the proposed method works well in practice.

1Department of Mathematics, Chemnitz Universitiy of Technology, Germany.kunis@mathematik.tu-chemnitz.de http://www.tu-chemnitz.de/∼skunis

2Faculty of Mathematics, University of Vienna, Austria.holger.rauhut@univie.ac.at http://homepage.univie.ac.at/holger.rauhut

33

Three-Dimensional Dyadic Wavelet Transform forImage Volume Denoising

Damiana Lazzaro1

Montefusco, L.B., Guerrini, C.

Reducing noise in images is a preliminary step in many image processing applications. In thiswork a method for edge-preserving denoising of three dimensional image data is presented. Itis based on the 3D generalization of the dyadic wavelet transform and overcomes the problemspresented by the extension to higher dimensions of the classical orthogonal wavelet transforms.In fact, a straightforward extension of these separable transforms to multidimensional data givesvery poor results, due to the lack of translation invariance that produces undesired artifacts.On the contrary, dyadic wavelets are truly multidimensional translation invariant bases, and thecoefficients of the dyadic expansion of an image f at different scales contain informations of thegradient of f smoothed according to the scale. By integrating these geometric informations in thesmoothing process we obtain a 3D denoising algorithm capable to preserve as much as possible ofthe signal features while reducing the noise to a sufficiently low level.

References

[1] S. Mallat and S. Zhong, Characterization of signals form multiscale edges, IEEE Transactionson Pattern Analysis and Machine Intelligence, vol. 14(7), 710–732, 1992.

[2] P. Mrazek, J.Weickert, and G. Steidl, Correspondences between wavelet shrinkage and nonlin-ear diffusion, L.D. Griffin and M. Lillholm (Eds.): Scale-Space 2003, LNCS 2695, pp. 101–116,2003. Springer-Verlag Berlin Heidelberg 2003.

[3] L.B.Montefusco, D. Lazzaro, Edge-Preserving Wavelet Thresholding for Image Denoising, toappear on Journal of Computationa and Applied Mathematics.

1Department of Mathematics, University of Bologna, Italy.damiana.lazzaro@unibo.it

1Department of Mathematics, University of Bologna, Italy.damiana.lazzaro@unibo.it

2Department of Mathematics, University of Bologna, Italy.laura.montefusco@unibo.it

3Department of Mathematics, University of Bologna, Italy.carla.guerrini@unibo.it

34

Integer ABS class of algorithms for solving linearDiophantine equations

Nezam Mahdavi-Amiri1

In the past two decades, a class of methods for solving linear systems of equations, first intro-duced by Abaffy, Broyden, and Spedicato (ABS) in 1984, has subsequently been extended to solvenonlinear systems of equations and optimization problems. These methods are a generalizationof the rank-one reducing process proposed by Egervary in 1955, and can also be explained by arank-one reducing formula given by Wedderburn in 1934. In 2001, we presented an ABS-typeclass of algorithms for solving linear Diophantine systems of equations giving complete character-izations of the integer solutions. Later in 2003, we extended the algorithms to solve the scaled(or preconditioned) Diophantine systems. Recently, we have shown how to efficiently solve a rankone modified Diophantine system using the information obtained from solving the original systemby an integer ABS algorithm. The updating strategy can be quite useful in implementing integerprogramming algorithms based on active set strategies. We describe the integer ABS approach ascited here.

References

[1] J. Abaffy, C.G. Broyden and E. Spedicato. A class of direct methods for linear systems.Numerische Mathematik, 45:361–376, 1984.

[2] J. Abaffy and E. Spedicato. ABS Projection Algorithms: Mathematical Techniques for Linearand Nonlinear Equations. Ellis Horwood, Chicheste, 1989.

[3] E. Egervary. Auflosung eines homogenen Diophantischen, gleichungsystems mit hilfe vonprojectormatrizen. Publ. Math. Debercen, 4:481-483, 1955.

[4] M.T. Chu, R.E. Funderlick and G.H. Golub. A rank-one reduction formula and its applicationto matrix factorizations. SIAM Review, 37(4):512–530, 1995.

[5] H. Esmaeili, N. Mahdavi-Amiri and E. Spedicato. A class of algorithms for Diophantine linearsystems. Numerische Mathematik, 90:101–115, 2001.

[6] E. Spedicato, E. Bodon, A. Del Popolo and N. Mahdavi-Amiri. ABS mathods and ABSPACKfor linear systems and optimization: A review. 4OR, 1:51–66, 2003.

[7] K. Amini and N. Mahdavi-Amiri. Solving rank one perturbed linear Diophantine systems bythe ABS methods. Optimization Methods and Software, to appear.

1Sharif University of Technology.nezamm@sharif.edu

35

Structured Low-Rank Approximation in SignalProcessing

Ivan Markovsky1

We review a representative sample of signal processing problems—linear prediction, harmonicretrieval, deconvolution, and image deblurring—and relate them to a generic structured low-rankapproximation problem. On one hand, these links show a conceptual unification of seeminglydifferent problems and put into evidence the central role of low-rank approximation in signal pro-cessing. On another hand, they allow to solve different problems with a single algorithm andtherefore use in practice the same piece of numerical software. We describe two algorithms forsolving the generic structured low-rank approximation problem. The first one is based on the vari-able projections method and the second one is based on the alternating projections method. Thealgorithm using the variable projections method is implemented and its performance is illustratedon simulation examples.

1Department of Electrical Engineering, Katholieke Universiteit Leuven, Belgium.Ivan.Markovsky@esat.kuleuven.be http://www.esat.kuleuven.be/sista-cosic-docarch/person.php?view=0&persid=34

36

A shifted Steihaug-Toint method for computing atrust-region step

Ctirad Matonoha1

Ladislav Luksan2

Jan Vcek

Introduction

Basic optimization methods for minimization of function F : Rn → R can be realized in variousways which differ in direction determination and step-size selection. Line-search and trust-regionrealizations are most popular. Trust-region methods can be advantageously used when the Hessianmatrix of the objective function (or its approximation) is indefinite, ill-conditioned or singular.This situation often arises in connection with the Newton method for general objective function(indefiniteness) or with the Gauss-Newton method for nonlinear least-squares (near-singularity).

Trust-region methods

Trust-region methods generate points xi ∈ Rn, i ∈ N , in such a way that x1 is arbitrary and

xi+1 = xi + αidi, i ∈ N , (4)

where di ∈ Rn are direction vectors and αi > 0 are step-sizes.A crucial part is a direction determination. There are various commonly known methods

for computing direction vectors satisfying certain conditions which we now mention briefly. Tosimplify the notation, we omit index i.

The most sophisticated method is based on a computation of the optimal locally constrainedstep. In this case, vector d ∈ Rn is obtained by solving subproblem

min Q(d) =1

2dTBd+ gT d subject to ‖d‖ ≤ ∆, (5)

where function Q(d) locally approximates difference F (xi + d) − F (xi). Necessary and sufficientconditions for this solution are

‖d‖ ≤ ∆, (B + λI)d = −g, B + λI º 0, λ ≥ 0, λ(∆− ‖d‖) = 0, (6)

where λ is a Lagrange multiplier. The More-Sorensen (MS) method [8] is based on solving non-linear equation 1/‖d(λ)‖ = 1/∆ with (B + λI)d(λ) + g = 0 by the Newton method using thesparse Choleski decomposition of B + λI. This method is very robust but requires 2-3 Choleskidecompositions per iteration.

Simpler methods are based on minimization of Q(d) on the two-dimensional subspace contain-ing Cauchy step dC = −(gT g/gTBg)g and Newton step dN = −B−1g. The most popular is thedog-leg (DL) method [2],[9], where d = dN if ‖dN‖ ≤ ∆ and d = (∆/‖dC‖)dC if ‖dC‖ ≥ ∆. Inthe remaining case, d is a convex combination of dC and dN such that ‖d‖ = ∆. This methodrequires only one Choleski decomposition per iteration.

If B is not sufficiently sparse, then the sparse Choleski decomposition of B is expensive. Inthis case, iterative methods based on conjugate gradients are more suitable. Steihaug [11] andToint [12] proposed a method based on the fact that Q(dk+1) < Q(dk) and ‖dk+1‖ > ‖dk‖ hold inthe subsequent CG iterations if CG coefficients are positive. We either obtain an unconstrainedsolution with a sufficient precision or stop on the trust-region boundary if a negative curvature

1Institute of Computer Science, Academy of Sciences of the Czech Republic, Praguematonoha@cs.cas.cz

2Institute of Computer Science, Academy of Sciences of the Czech Republic, Czech.luksan@cs.cas.cz http://www.cs.cas.cz/∼luksan/

37

is indicated or the trust-region is left. When suitable preconditioning is used, then this method(PST) is very efficient in practice. Note that ‖dk+1‖C > ‖dk‖C (where ‖dk‖

2C = dTkCdk) holds

instead of ‖dk+1‖ > ‖dk‖ if preconditioner C (symmetric and positive definite) is used. Thus thesolution on the trust-region boundary obtained by the preconditioned CG method can be fartherfrom the optimal locally constrained step than the solution obtained without preconditioning (seeFigure 1). This insufficiency is usually compensated by the rapid convergence of the preconditionedCG method.

||d|| < ∆x

||d||C < ∆

Figure 1: Preconditioned CG method.

The solution on the trust-region boundary obtained by the Steihaug-Toint method can berather far from the optimal solution. This insufficiency can be overcame by using the Lanczosprocess [3]. Initially, the conjugate gradient algorithm is used as in the Steihaug-Toint method.At the same time, the Lanczos tridiagonal matrix is constructed from the CG coefficients. If anegative curvature is indicated or the trust-region is left, we turn to the Lanczos process. In thiscase, d = Zd, where d is obtained by solving subproblem

min1

2dTT d+ ‖g‖eT1 d subject to ‖d‖ ≤ ∆. (7)

Here T = ZTBZ (with ZTZ = I) is the Lanczos tridiagonal matrix and e1 is the first columnof the unit matrix. This method cannot be successfully preconditioned, since preconditioningchanges trust-region constraint ‖d‖ ≤ ∆ to ‖d‖C ≤ ∆, where C changes in each major iterationand can be ill-conditioned.

Therefore, we apply the Steihaug-Toint method to subproblem

min Q(d) = Qλ(d) =1

2dT (B + λI)d+ gT d subject to ‖d‖ ≤ ∆. (8)

Number λ ≥ 0, which approximates λ in (6), is found by solving a small-size subproblem (7) withtridiagonal matrix T obtained by using a small number of the Lanczos steps. This method [6],[7]like method [3], combines good properties of the More-Sorensen and the Steihaug-Toint methods.Moreover, it can be successfully preconditioned. The point on the trust-region boundary obtainedby this method is usually closer to the optimal solution in comparison with the point obtained bythe original Steihaug-Toint method.

A shifted Steihaug-Toint method

A (preconditioned) shifted Steihaug-Toint method (PSST) differs from the standard one byusing shifted subproblem (8), where number λ approximates λ in (6). Number λ has to be chosen

38

in such a way that λ = 0 if ‖d‖ < ∆, where d is a solution of (5), which is true if 0 ≤ λ ≤ λ.If we denote Kk = spang,Bg, . . . , Bk−1g the Krylov subspace of dimension k, then (under

some assumptions) we can prove the following assertions. Let

dk(λi) = arg mind∈Kk

Qλi(d), where Qλ(d) =

1

2dT (B + λI)d+ gT d.

Thenλi ≤ λj ⇔ ‖dk(λi)‖ ≥ ‖dk(λj)‖.

Moreover, if

dj = arg mind∈Kj

Q(d) subject to ‖d‖ ≤ ∆, where Q(d) =1

2dTBd+ gT d

with corresponding Lagrange multipliers λj , j ∈ 1, . . . , n, then for 1 ≤ k ≤ l ≤ n we have

λk ≤ λl.

Let’s return to subproblem (8). If we set λ = λk for some k ≤ n, then 0 ≤ λ = λk ≤ λn = λ.As a consequence of this inequality, one has that λ = 0 implies λ = 0, so that ‖d‖ < ∆ impliesλ = 0. Thus the shifted Steihaug-Toint method reduces to the standard one in this case. At thesame time, if B is positive definite and λ > 0, then one has ∆ ≤ ‖(B + λI)−1g‖ < ‖B−1g‖. Thusthe unconstrained minimizer of (8) is closer to the trust-region boundary than the unconstrainedminimizer of (5) and we can expect that d(λ) is closer to the optimal locally constrained step thand. Finally, if λ > 0, then matrix B + λI is better conditioned than B and we can expect that theshifted Steihaug-Toint method will converge more rapidly than the original one.

The shifted Steihaug-Toint method consists of the three major steps. First, we carry out k ¿ nsteps of the unpreconditioned Lanczos method [3] to obtain tridiagonal matrix T ≡ Tk = ZT

k BZk,where Zk ∈ R

n×k is the matrix whose columns form an orthonormal basis in Kk. Then we solvesubproblem (7) using the More-Sorensen method [8] to obtain Lagrange multiplier λ. Finally, weapply the (preconditioned) Steihaug-Toint method [11],[12] to subproblem (8) to obtain directionvector d = d(λ).

Computational experiments

A numerical comparison of methods for computing direction vectors mentioned in subsection2 implies several conclusions [6],[7]. If problems do not have sparse Hessian matrices, then directmethods MS and DL can be much worse than iterative methods PST and PSST. On the otherhand, direct methods can be more efficient for ill-conditioned but reasonably sparse problems.Comparing PST and PSST, we can see that PSST is usually slightly worse than PST, measuredby the computational time, since it uses additional operations for determining the Lanczos matrixT and computing parameter λ. Nevertheless, if the problems are difficult, then PSST is muchbetter than PST. Thus the total computational time can be lower for PSST.

• MS - the method of More and Sorensen [8] for computing the optimal locally constrainedstep.

• DL - the dog-leg strategy of Powell [9] or Dennis and Mei [2].

• MDL - the multiple dog-leg strategy (k = 5) mentioned in [11].

• ST - the basic (unpreconditioned) Steihaug [11] and Toint [12] method.

• GLRT - the method of Gould, Lucidi, Roma and Toint [3] which combines CG method withthe Lanczos process to give a good approximation of the optimal locally constrained step.

• PST - Preconditioned Steihaug-Toint method. The incomplete Choleski preconditioner isused.

39

• PSST - Preconditioned shifted Steihaug-Toint method (k = 5). The incomplete Choleskipreconditioner is used.

These algorithms were used for solving trust-region subproblems arising in a realization ofa discrete Newton’s method. They were tested by using a set of 22 sparse least-squares testproblems with 1000 and 5000 variables (subroutine TEST14 [5], which can be found on the pagewww.cs.cas.cz/~luksan/test.html). The results are given in Table 1, where NIT is the totalnumber of iterations, NFV is the total number of function evaluations, NFG is the total number ofgradient evaluations, NCG is the total number of CG iterations and Time is the total computationaltime (in seconds).

N Method NIT NFV NFG NCG Time

1000 MS 1918 1955 8797 - 4.65DL 2515 2716 11859 - 4.42MDL 2292 2456 10673 12203 4.61ST 3329 3784 16456 53573 8.20GLRT 3107 3444 15306 55632 8.53PST 2631 2823 13019 910 5.14PSST 1999 2046 9201 1161 4.25

5000 MS 8391 8566 35824 - 122.44DL 9657 10133 42425 - 115.77MDL 8938 9276 39032 47236 122.84ST 16894 19163 83933 358111 364:42GLRT 14679 16383 71483 366695 401.45PST 10600 11271 50365 3767 145.42PSST 8347 8454 35939 4329 108.87

Table 1: Comparison of methods using TEST14.

For a better comparison of methods PST, PSST, DL and MS, we have performed additionaltests with problems from the widely used CUTE collection [1] which can be found in [6],[7].

Acknowledgement: This work was supported by the Grant Agency CAS, project IAA1030405,the Grant Agency CR, project 201/06/P397, and the Institutional Research Plan AV0Z10300504.

References

[1] I.Bongartz, A.R.Conn, N.Gould, P.L.Toint: CUTE: constrained and unconstrained testingenvironment. ACM Transactions on Mathematical Software 21 (1995), pp. 123-160.

[2] J.E.Dennis, H.H.W.Mei: An unconstrained optimization algorithm which uses function andgradient values. Report No. TR 75-246, 1975.

[3] N.I.M.Gould, S.Lucidi, M.Roma, P.L.Toint: Solving the trust-region subproblem using theLanczos method. Report No. RAL-TR-97-028, 1997.

[4] L.Luksan, M.Tuma, J.Vlcek, N.Ramesova, M.Siska, J.Hartman, C.Matonoha: UFO 2004 –Interactive System for Universal Functional Optimization. Technical report no. V923-04, ICSAS CR, 2004, 227 p.

[5] L.Luksan, J.Vlcek: Sparse and partially separable test problems for unconstrained and equalityconstrained optimization. Report No. V767-98, Institute of Computer Science AS CS, 1998.

[6] L.Luksan, C.Matonoha, J.Vlcek: A shifted Steihaug-Toint method for computing a trust-region step. Technical report No. V914-04, ICS AS CR, 2004, 9 p.

40

[7] L.Luksan, C.Matonoha, J.Vlcek: A shifted Steihaug-Toint method for computing a trust-region step. To appear in BIT.

[8] J.J.More, D.C.Sorensen: Computing a trust region step. SIAM Journal on Scientific andStatistical Computations 4 (1983), pp. 553-572.

[9] M.J.D.Powell: A new algorithm for unconstrained optimization. In: Nonlinear Programming,J.B.Rosen, O.L.Mangasarian, K.Ritter, eds., Academic Press, London, 1970.

[10] M.J.D.Powell: On the global convergence of trust region algorithms for unconstrained opti-mization. Mathematical Programming 29 (1984), pp. 297-303.

[11] T.Steihaug: The conjugate gradient method and trust regions in large-scale optimization.SIAM Journal on Numerical Analysis 20 (1983), pp. 626-637.

[12] P.L.Toint: Towards an efficient sparsity exploiting Newton method for minimization. In:Sparse Matrices and Their Uses, I.S.Duff, ed., Academic Press, London, 1981, pp. 57-88.

41

Regularizing inverse preconditioners forsymmetric band Toeplitz matrices

Ornella Menchi1

Paola Favati2

Grazia Lotti3

The image restoration is a widely studied discrete ill-posed problem. Among the many regu-larization methods used for treating the problem efficiently the iterative methods have been shownto be effective. We consider here the case of a blurring function defined by space invariant andband limited PSF, which is modelled by a linear system having a band block Toeplitz structurewith band Toeplitz blocks. In order to reduce the number of iterations required to obtain accept-able reconstruction, an Inverse Toeplitz preconditioner for problems with Toeplitz structure hasbeen proposed by Hanke and Nagy [1]. The cost per iteration of their method is of O(n2 log n)operations, where n2 is the pixels number of the 2D image, regardless of the band structure of theoriginal matrix. We propose inverse preconditioners with band Toeplitz structure, which lowerthe cost to O(n2) and show experimentally the same speed of convergence and reconstructionefficiency of the Inverse Toeplitz preconditioner.

References

[1] M. Hanke and J. Nagy, Inverse Toeplitz preconditioner for ill-posed problems. Linear AlgebraAppl., 284:137–156, 1994.

1Dipartimento di Informatica, Universita di Pisa, Italy.menchi@di.unipi.it

2IIT - CNR Via G. Moruzzi 1, 56124 Pisa, Italy.paola.favati@iit.cnr.it

3Department of Mathematics, University of Parma, Parco Area delle Scienze 53/A, 43100 Parma, Italy.grazia.lotti@unipr.it

42

To be announced

Vadim Olshevsky1

1Department of Mathematics, University of Connecticut, USAolshevsky@math.uconn.edu http://www.math.uconn.edu/∼olshevsky/

43

A Doubly iterative wavelet algorithm for edgepreserving image Reconstruction

Serena Papi1

Bacchelli Silvia2

Many image processing problems are ill posed and must be regularized. This is the case ofimage reconstruction in computerized tomography. Tomographic image reconstruction refers to abroad class of problems in which the goal is to recover a two-dimensional object of interest froma set of measurements that may be described as noisy line integrals, or approximations thereof.The degradation model is represented by the linear equation

g = Kf + e

where Kf represents the Radon Operator acting on the two-dimensional image function f , namely

Kf =

∫ +∞

−∞

∫ +∞

−∞

f(x, y) δ(x cosθ + y sinθ − t) dx dy

and e is the independent and identically distributed term of white noise with mean zero and vari-ance σ2. The aim of this work is to consider a variational approach in a deterministic waveletsetting, consisting in the minimization of a functional that incorporates both a discrepancy term,ensuring that the estimated f is a good approximation of the original image, and a term repre-senting a priori knowledge one may have about the solution. This latter term enforces a roughnesspenalty on the estimate, and, in order to obtain an edge-preserving reconstruction, it is definedas a sum of potentials which are functions of image variations.In this paper we approach the variational problem by means of an iterative algorithm, similarin the spirit to that proposed in [1], which amounts to a Landweber iteration with a non linearshrinkage applied at each iteration step. Differently from [1], we work in the wavelet packet domain,which is particularly suited to drastically reduce the computational complexity of the algorithm,while exploiting the wavelet packet capability of decorrelating the noise from the significant signalfeatures, also in the highest frequency subbands. Moreover, at each iteration step, the amountof shrinkage depends on the value of the wavelet coefficients that are arguments of the potentialfunction We get so a non linear equation for each wavelet coefficient that we solve by using aniterative procedure. This approach hence produces a doubly iterative algorithm which conjugatethe effectiveness of the wavelet packet bases in image processing with the edge preserving propertiesof the used potential functions.

References

[1] I.Daubechies, M. Defrise, C. De Mol An iterative thresholding algorithm for linear inverseproblems with a sparsity constraint, Communication on Pure and Applied Mathematics,57,11, pp.1413-1457.

[2] P. Charbonnier, L. Blanc-Feraud, G. Aubert, M. Barlaud, Deterministic Edge-PreservingRegularization in Computed Imaging, IEEE Transaction on Image Processing, vol. 6, 2, 1997,pp.298-311.

[3] L.B.Montefusco, D. Lazzaro, Edge-Preserving Wavelet Thresholding for Image Denoising,preprint.

1Department of Mathematics, University of Bolognaserena.papi@unibo.it

2Department of Mathematics, University of Bologna, Italy.silvia@csr.unibo.it

44

Iterative Methods for Inverse Problems

Lothar Reichel1

Inverse problem are concerned with the determination of the cause of an observed or desiredeffect. Many inverse problems can be formulated as a Fredholm integral equation of the first kind,whose solution is an ill-posed problem. This talk discusses numerical methods for the solution oflarge-scale ill-posed inverse problems by iterative methods.

1Department of Mathematical Sciences, Kent State University Kent, USA.reichel@math.kent.edu http://www.math.kent.edu/∼reichel/

45

Structure-exploiting algorithms in systemidentification

Vasile Sima1

The availability of powerful system identification tools is very important in practice, since mod-ern control techniques are heavily dependent on suitable dynamical models. Efficient, structure-exploiting algorithms for input/output (I/O) data processing in subspace-based multivariable sys-tem identification are presented. These algorithms use a matrix H, built from concatenated block-Hankel matrices, using part of the available I/O data. The first step in a subspace identificationprocedure is to perform a “data compression” by computing an upper triangular factor R of a QRfactorization of the matrix H. Since many theoretical results concerning, e.g., consistency and nor-mality of the estimates of the system matrices hold asymptotically, the algorithms should be ableto deal with large amounts of data, if available. The standard QR factorization algorithm could betoo costly, since H can have very large dimensions. Efficient data processing is possible by usingthe problem structure. Two such techniques are based on fast Cholesky and fast QR factorizationalgorithms, which exploit the special structure of the matrix H. The algorithms are implementedin the system identification toolbox for discrete-time systems, SLIDENT, incorporated in the For-tran 77 Subroutine Library In COntrol Theory (SLICOT). Besides drivers and computationalroutines, this toolbox provides MATLAB or Scilab interfaces, implementing several algorithmicapproaches. Extensive numerical testing and comparisons with similar MATLAB tools show thatSLIDENT is reliable, efficient, and powerful enough to solve industrial identification problems.

References

[1] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Green-baum, S. Hammarling, A. McKenney and D. Sorensen. LAPACK Users’ Guide: Third Edi-tion, 1999. SIAM, Philadelphia.

[2] P. Benner, V. Mehrmann, V. Sima, S. Van Huffel and A. Varga. SLICOT — A subroutinelibrary in systems and control theory. In: Applied and Computational Control, Signals, andCircuits (B. N. Datta, Ed.), Vol. 1, Chap. 10, pp. 499–539, 1999. Birkhauser, Boston.

[3] C. Gomez (Ed.). Engineering and Scientific Computing with Scilab, 1999. Birkhauser, Boston.

[4] W. Larimore. System identification, reduced order filtering and modeling via canonical variateanalysis. In: Proc. of the American Control Conference, San Francisco, CA, USA, pp. 445–451, 1983.

[5] Ljung, L. System Identification Toolbox For Use with MATLAB. User’s Guide, Version 5.The MathWorks, Inc. 3 Apple Hill Drive, Natick, MA, 2000.

[6] N. Mastronardi, D. Kressner, V. Sima, P. Van Dooren and S. Van Huffel. A fast algorithmfor subspace state-space system identification via exploitation of the displacement structure.J. Comput. Appl. Math., 132(1), 71–81, 2001.

[7] The MathWorks, Inc. Using MATLAB. Version 5. 3 Apple Hill Drive, Natick, MA., 1999.

[8] V. Sima. Cholesky or QR factorization for data compression in subspace-based identification?In: Proc. of the Second NICONET Workshop on “Numerical Control Software: SLICOT, aUseful Tool in Industry”, December 3, 1999, INRIA Rocquencourt, France, pp. 75–80, 1999.

1National Institute for Research & Development in Informatics.vsima@ici.ro http://www.ici.ro/ici/organizare/directory/vsima.html

46

[9] V. Sima. SLICOT linear systems identification toolbox. SLICOT Working Note 2000-4,Katholieke Universiteit Leuven (ESAT/SISTA), Leuven, Belgium, 2000. Available fromhttp://www.slicot.org.

[10] V. Sima, and S. Van Huffel. Performance investigation of SLICOT system identificationtoolbox. In: Proc. of the European Control Conference, ECC 2001, 4–7 September, 2001,Seminario de Vilar, Porto, Portugal, pp. 3586–3591, 2001.

[11] V. Sima, D. M. Sima and S. Van Huffel. High-performance numerical algorithms and softwarefor subspace-based linear multivariable system identification. J. Comput. Appl. Math, 2004.

[12] V. Sima, D. M. Sima and S. Van Huffel. SLICOT system identification software and appli-cations. In: Proc. of the 2002 IEEE International Conference on Control Applications andIEEE International Symposium on Computer Aided Control System Design, CCA/CACSD2002, September 18–20, 2002, Scottish Exhibition and Conference Centre, Glasgow, Scotland,U.K., pp. 45–50, 2002. Omnipress.

[13] P. Van Overschee, and B. De Moor. N4SID: Two subspace algorithms for the identificationof combined deterministic-stochastic systems. Automatica, 30(1), 75–93, 1994.

[14] P. Van Overschee, and B. De Moor. Subspace Identification for Linear Systems : Theory– Implementation – Applications. Kluwer Academic Publishers, Boston/London/Dordrecht,1996.

[15] M. Verhaegen. Subspace model identification. Part 3: Analysis of the ordinary output-errorstate-space model identification algorithm. Int. J. Control, 58(3), 555–586, 1993.

[16] M. Verhaegen, and P. Dewilde. Subspace model identification. Part 1: The output-error state-space model identification class of algorithms. Int. J. Control, 56(5), 1187–1210, 1992.

47

Quantification of magnetic resonancespectroscopic signals and separable constrained

nonlinear least squares

Diana Sima1

Sabine Van Huffel2

A magnetic resonance spectroscopic (MRS) time-domain signal can capture, in a non-invasiveway, relevant biochemical information from the human body. The quantification of all substancesthat are present in a tissue (such as a small region in the brain), is important for medical diag-nosis. An MRS signal can be mathematically modeled as a linear combination of some nonlinearfunctions. Fitting this model to an in vivo measured signal leads to a separable nonlinear leastsquares fitting problem, with linear bound constraints on some variables. The Variable Projection(VARPRO) [1] technique can be applied to this problem, but needs to be adapted in several re-spects. If only the nonlinear variables are subject to constraints, then the Levenberg-Marquardtminimization algorithm that is classically used by the VARPRO method should be replaced with aversion that can incorporate those constraints. If some of the linear variables are also constrained,then they cannot be projected out in closed-form expression, as is the case of the classical VARPROtechnique. We show how quadratic programming problems can be solved instead, and we providedetails on efficient function and approximate Jacobian evaluations for the inequality constrainedVARPRO method.

References

[1] G. H. Golub and V. Pereyra. “The differentiation of pseudo-inverses and nonlinear leastsquares problems whose variables separate,” SIAM Journal on Numerical Analysis, 10:413–432, 1973.

1Department of Electrical Engineering, Katholieke Universiteit Leuven, Belgium.diana.sima@esat.kuleuven.be http://www.esat.kuleuven.ac.be/sista-cosic-docarch/person.php?view=1&persid=159

2Department of Electrical Engineering, Katholieke Universiteit Leuven, Belgium.Sabine.VanHuffel@esat.kuleuven.behttp://www.esat.kuleuven.ac.be/sista/members/vanhuffel.html

48

Lanczos tridiagonalization and Golub-Kahanbidiagonalization: Ideas, connections and impact

Zdenek Strakos1

The Lanczos algorithm was introduced in the middle of the last century in the context ofcomputing eigenvalues and solving linear algebraic systems. It can be viewed in many differentways — for example as a (partial) reduction of the given real symmetric matrix to tridiagonalform. The closely related Golub-Kahan bidiagonalization was introduced in 1965 in the context ofcalculating the singular value decomposition, and applies to a general rectangular matrix. Bothalgorithms have influenced the theory and practice of scientific computing immensely. We willnot survey the multitude of related algorithms and their implementations here. Instead, thispresentation will concentrate on several mathematical ideas which link together both differentperiods, and different disciplines, and will concentrate on a model reduction aspect.

The list of authors who contributed significantly to the development and understanding of linearalgebraic and/or eigenvalue solvers based on the Lanczos tridiagonalization and the Golub-Kahanbidiagonalization would be very long. Various successful applications, descriptions of convergencebased on the relationship to Gauss quadrature, and rounding error analyses which have led to ourunderstanding of the behavior in finite precision arithmetic, certainly belong among the importantachievements of numerical linear algebra in the second half of the 20th century. Related algorithmshave been widely used in Statistics, Engineering, and the Sciences; with a very important role inComputational Physics and Computational Quantum Chemistry. From more recent developmentsit is worth mentioning, e.g., the bidiagonalization-based methods for solving ill-posed problems,model-reduction applications and the core problem formulation in errors-in-variables modelling.

In our contribution we recall some basic historical and interdisciplinary links. We close bygiving examples of seemingly unrelated results. These demonstrate how a classical analytic viewcan be useful in approaching very recently formulated open problems, and vice versa. These openproblems are motivated by computational issues — such as rates of convergence, and effects ofrounding errors.

This work has been supported by the National Program of Research under project 1ET400300415.

1Institute of Computer Science, Academy of Sciences of the Czech Republic, Praha.strakos@cs.cas.cz http://www.cs.cas.cz/∼strakos/

49

A subspace based structure preserving solutionmethod for the Electrical Impedance Tomography

problem: Part I

Bart Truyen1

Luca Dimiccoli2

Jan Cornelisfootnotemark[3]

We present a subspace based, structure preserving solution method for the problem of ElectricalImpedance Tomography, where the conductivity inside a simply connected 2-dimensional domainis sought from noisy and incomplete boundary data. Unlike conventional output-least squaresalgorithms that can be regarded as minimizing a certain error norm, solutions are recovered hereas the minimizers of a closely related residual norm problem. An iterative solution scheme is shownto lead to a sequence of sparse matrix subproblems, with conditioning far more favourable thantypically observed in output-least squares. We find that these sparse subproblems demonstrate aparticular form of displacement structure that can be further elaborated to finally arrive upon anefficient computational implementation. In the first part of this contribution, we introduce thestructured problem formulation, outline the algorithmic approach taken, and summarize some ofits numerical properties.

1ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUBbatruyen@etro.vub.ac.be http://www.etro.vub.ac.be/Personal/btruyen/

2ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUB.ldimicco@etro.vub.ac.be http://www.etro.vub.ac.be/Members/Personal Page.asp?PM ID=1463

3ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUB.

50

Fast dimensionality reduction and approximatesimultaneous diagonalization of matrices

Eugene Tyrtyshnikov1

Ivan Oseledets2

Dmitry Savostianov3

In some applications we are interested to work with huge amounts of data for which standardmethods of numerical analysis cannot be applied. For instance, consider a 3-dimensional arrayof size of several thousand in every dimension. Is it possible to work with this size array anyefficiently? In the talk, we are going to show that the answer is ”yes” - however, under certainassumptions on the data.

The main assumption is that the data admits a sufficiently accurate low-tensor-rank approxima-tion. Then, we present a special cross-approximation technique using a relatively small amount ofthe entries of the given array and providing the so-called Tucker decomposition with the complexitylinear in one-dimensional size of the array. The technique generalizes the low-rank approximationtechnique for matrices to 3-dimensional arrays. The Tucker core can be further approximated bya trilinear (canonical) decomposition. Development of good algorithms solving the latter problemis still topical. To this end, we present a new algorithm that seems to be competitive with thebest known algorithms.

1Institute of Numerical Mathematics, Russian Academy of Sciencestee@inm.ras.ru http://bach.inm.ras.ru/chair.htm

2Institute of Numerical Mathematics, Russian Academy of Sciences.3Institute of Numerical Mathematics, Russian Academy of Sciences.

51

Optimal correlation between projected matrices

Paul Van Dooren1

Catherine Fraikin2

Yurii Nesterov2

We consider the problem of finding the optimal correlation between two projected matrices U ∗AUand V ∗BV . The square matrices A and B may be of different dimensions, but the isometriesU and V have a common column dimension k. The correlation is measured by the real functionc(U, V ) = < Tr(U∗AUV ∗B∗V ), which we maximize of the isometries U ∗U = V ∗V = Ik.

This problem can be viewed as an extension of the generalized numerical range of two ma-trices, which are now allowed to be of different dimension. We discuss several properties of thisoptimization problem, characterize its extremal points and propose an algorithm converging tosuch an extremal point.

1Department of Mathematical Engineering, Catholic University of Louvain, Belgium.vdooren@csam.ucl.ac.be http://www.inma.ucl.ac.be/∼vdooren/

2Universite catholique de Louvain, Department of Mathematical Engineering, Louvain-la-Neuve, Belgium.

52

Numerical linear algebra in computationalbiomedical signal processing

Sabine Van Huffel1

In biomedical signal processing, the aim is to extract clinically, biochemically or pharmaceu-tically relevant information (e.g metabolite concentrations in the brain) in terms of parametersout of low-quality measurements in order to enable an improved medical diagnosis. Typically,biomedical data are affected by large measurement errors, largely due to the non-invasive natureof the measurement process or the severe constraints to keep the input signal as low as possible forsafety and bio-ethical reasons. Accurate and automated quantification of this information requiresan ingenious combination of the following 4 issues:

• an adequate pretreatment of the data,

• the design of an appropriate model and model validation,

• a fast and numerically robust model parameter quantification method,

• an extensive evaluation and performance study, using in-vivo and patient data, up to theembedding of the advanced tools into user-friendly user interfaces to be used by clinicians.

The underlying computational signal processing problems can be solved by making use oflinear algebra, signal processing, linear algebra, system theory and optimisation. In particular,it is shown how computational linear algebra kernels, such as SVD, PCA, Least Squares,TLS,CCA, ICA, ..., can be used as building blocks of higher-level signal processing algorithms. Inaddition, the application of these algorithms will be briefly illustrated in a variety of case studies,including Magnetic Resonance Spectroscopic Imaging, and epileptic seizure detection. In this talkwe present a new method for solving the Toeplitz system

1Department of Electrical Engineering, Katholieke Universiteit Leuven, Belgium.Sabine.VanHuffel@esat.kuleuven.behttp://www.esat.kuleuven.ac.be/sista/members/vanhuffel.html

53

The use of lower semiseparable matrices to solvegeneralized eigenvalue problems

Yvette Vanberghen1

Raf Vandebril2

Marc Van Barel3

Lower semiseparable matrices have already been used as an alternative to Hessenberg matricesfor computing the eigenvalues.

In this talk we will discuss the use of lower semiseparable matrices when solving generalizedeigenvalue problems. Given a matrix pair (A,B), the matrix A is classically reduced to Hessenbergform by unitary transformations, while the matrix B is reduced to upper triangular form. As analternative, we will reduce the matrix A into lower semiseparable form. Furthermore we will designa QZ-step that preserves the (lower semiseparable, upper triangular) structure. The algorithmswill be illustrated with some numerical examples.

1Department of Computer Science, Katholieke Universiteit Leuven, Belgiumyvette.vanberghen@cs.kuleuven.be

2Department of Computer Science, Katholieke Universiteit Leuven, Belgium.raf.vandebril@cs.kuleuven.be http://www.cs.kuleuven.be/∼raf/

3Department of Computer Science, Katholieke Universiteit Leuven, Belgium.marc.vanbarel@cs.kuleuven.be http://www.cs.kuleuven.ac.be/∼marc/

54

Reducing large gyroscopic eigenproblems byautomated multi-level substructuring

Heinrich Voss1

Simulating numerically the sound radiation of a rolling tire requires the solution of a verylarge and sparse gyroscopic eigenvalue problem. Taking advantage of the automated multi–levelsubstructuring (AMLS) method it can be projected to a much smaller gyroscopic problem, thesolution of which however is still quite costly since the eigenmodes are non–real and complex arith-metic is necessary. This paper discusses the application of AMLS to huge gyroscopic problems andthe numerical solution of the AMLS reduction. A numerical example demonstrates the efficiencyof AMLS. This is joint work with Kolja Elssel.

1Section of Mathematics, Technical University of Hamburg-Harburg, Germany.voss@tu-harburg.de http://www.tu-harburg.de/mat/hp/voss/

55

8 List of participants

Adesokan, Bolaji, African institute for mathematical sciences.www.aims.za.ac/∼bjaysokanbj@yahoo.com

Agreste, Santa, Department of Information Science, Milan University, Italy.agreste@dsi.unimi.it

Arico, Antonio, Dipartimento di Fisica e Matematica, Universita dell’Insubria, Como, Italy.antonio.arico@unipv.it

Benner, Peter, Department of Mathematics, Chemnitz University of Technology, Chemnitz,Germany.http://www-user.tu-chemnitz.de/∼benner/benner@mathematik.tu-chemnitz.de

Benzi, Michele, Department of Mathematics and Computer Science, Emory University,Atlanta, USA.http://www.mathcs.emory.edu/∼benzi/benzi@mathcs.emory.edu

Bini, Dario A., Dipartimento di Matematica, Universita degli Studi di Pisa, Italy.http://www.dm.unipi.it/∼bini/bini@dm.unipi.it

Chan, Raymond, Department of Mathematics, The Chinese University of Hong Kong, HongKong, China.http://www.math.cuhk.edu.hk/∼rchanrchan@math.cuhk.edu.hk

Chang, Xiao-Wen, School of Computer Science, McGill University, Montreal, Canada.http://www.cs.mcgill.ca/∼chang/index.phpchang@cs.mcgill.ca

Dax, Achiya, Hydrological Service, P.O.B. 36118, Jerusalem 91360, ISRAEL.dax20@water.gov.il

De Wilde, Patrick, Faculty of Electrical Engineering, Technische Universiteit Delft, Delft,The Netherlands.http://cobalt.et.tudelft.nl/people/dewilde.htmldewilde@cas.et.tudelft.nl

Delvaux, Steven, Department of Computer Science, Katholieke Universiteit Leuven, Belgium.http://www.cs.kuleuven.be/∼stevenSteven.Delvaux@cs.kuleuven.be

Diele, Fasma, Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche,Bari, Italy.http://www.ba.cnr.it/∼irmafd03/Welcome.htmlf.diele@ba.iac.cnr.it

Dimiccoli, Luca, ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUB.http://www.etro.vub.ac.be/Members/Personal Page.asp?PM ID=1463ldimicco@etro.vub.ac.be

Donatelli, Marco, Dipartimento di Fisica e Matematica, Universita dell’Insubria, Como,Italy.marco.donatelli@uninsubria.it

56

Duintjer Tebbens, Jurjen, Institute of Computer Science, Academy of Sciences of the CzechRepublic.tebbens@cs.cas.cz

Eidelman, Yuli, Department of Mathematics, Tel-Aviv University, Israel.eideyu@post.tau.ac.il

Ellenbroek, Rogier, Delft University of Technology, Delft, The Netherlands.http://www.dcsc.tudelft.nl/∼rellenbroekr.m.l.ellenbroek@dcsc.tudelft.nl

Fasino, Dario, Dipartimento di Matematica e Informatica Universit degli Studi di Udine,Italy.http://www.dimi.uniud.it/fasino/fasino@dimi.uniud.it

Fassbender, Heike, Institut Computational Mathematics, TU Braunschweig, Germany.http://www-public.tu-bs.de:8080/∼hfassben/index en.htmlh.fassbender@tu-bs.de

Hansen, Per Christian, Informatics and Mathematical Modelling, Technical University ofDenmark.http://www.imm.dtu.dk/∼pchpch@imm.dtu.dk

Hnetynkova, Iveta, Dep. of Numerical Mathematics, Faculty of Mathematics and Physics,Charles University, Prague.iveta.hnetynkova@gmail.com

Hueper, Knut, National ICT Australia Ltd.http://users.rsise.anu.edu.au/∼hueper/knut.hueper@nicta.com.au

Husøy, John Hakon, University of Stavanger.john.h.husoy@uis.no

Keiner, Jens, Institute of Mathematics, University of Lubeck, Germany.http://www.math.uni-luebeck.de/keinerkeiner@math.uni-luebeck.de

Kunis, Stefan, Department of Mathematics, Chemnitz Universitiy of Technology, Germany.http://www.tu-chemnitz.de/∼skuniskunis@mathematik.tu-chemnitz.de

Lamura, Antonio, Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche,Bari, Italy.http://www.ba.cnr.it/∼irmaal21/Welcome.htmla.lamura@ba.iac.cnr.it

Laudadio, Teresa, Istituto di studi sui sistemi intelligenti per l’automazione, Consiglio Nazionaledelle Ricerche, Bari, Italy.http://www.ba.cnr.it/∼irman21/Welcome Teresa.htmllaudadio@ba.issia.cnr.it

Lazzaro, Damiana, Department of Mathematics, University of Bologna, Italy.damiana.lazzaro@unibo.it

Mahdavi-Amiri, Nezam, Sharif University of Technology.nezamm@sharif.edu

57

Marangi, Carmela, Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche,Bari, Italy.c.marangi@ba.iac.cnr.it

Markovsky, Ivan, Department of Electrical Engineering, Katholieke Universiteit Leuven, Bel-gium.http://www.esat.kuleuven.be/sista-cosic-docarch/person.php?view=0&persid=34Ivan.Markovsky@esat.kuleuven.be

Mastronardi, Nicola, Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche,Bari, Italy.http://www.ba.cnr.it/∼irmanm21/Welcome.htmln.mastronardi@ba.iac.cnr.it

Matonoha, Ctirad, Institute of Computer Science, Academy of Sciences of the Czech Repub-lic, Prague.matonoha@cs.cas.cz

Menchi, Ornella, Dipartimento di Informatica, Universita di Pisa, Italy.menchi@di.unipi.it

Nico, Giovanni, Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche,Bari, Italy.g.nico@ba.iac.cnr.it

Notarnicola, Filippo, Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delleRicerche, Bari, Italy.f.notarnicola@ba.iac.cnr.it

Olshevsky, Vadim, Department of Mathematics, University of Connecticut, USA.http://www.math.uconn.edu/∼olshevsky/olshevsky@math.uconn.edu

Papi, Serena, Department of Mathematics, University of Bologna.serena.papi@unibo.it

Pavel, Jiranek, Department of Modelling of Processes, Technical University of Liberec, CzechRepublic.http://flow.kmo.tul.cz/∼paya/pavel.jiranek@tul.cz

Plesinger, Martin, Dept. of Modelling of Processes, Technical University of Liberec, CzechRepublic.martin.plesinger@seznam.cz

Poloni, Federico, student at: Scuola Normale Superiore, Pisa.f.poloni@sns.it

Reichel, Lothar, Department of Mathematical Sciences, Kent State University Kent, USA.http://www.math.kent.edu/∼reichel/reichel@math.kent.edu

Rodriguez, Giuseppe, Dipartimento di Matematica, Universit di Cagliari, Italy.http://bugs.unica.it/∼gppe/rodriguez@unica.it

Sayed, Ali H., Electrical Engineering Department, University of California, Los Angeles,USA.http://www.ee.ucla.edu/faculty/bios/sayed.htmsayed@ee.ucla.edu

58

Serra-Capizzano, Stefano, Dipartimento di Fisica e Matematica, Universit dell’Insubria,Como, Italy.http://www3.uninsubria.it/pls/uninsubria/uninsubria docenti.h preview?id doc=P000411s.serracapizzano@uninsubria.it

Sima, Diana, Department of Electrical Engineering, Katholieke Universiteit Leuven, Bel-gium.http://www.esat.kuleuven.ac.be/sista-cosic-docarch/person.php?view=1&persid=159diana.sima@esat.kuleuven.be

Sima, Vasile, National Institute for Research & Development in Informatics.http://www.ici.ro/ici/organizare/directory/vsima.htmlvsima@ici.ro

Strakos, Zdenek, Institute of Computer Science, Academy of Sciences of the Czech Republic,Praha.http://www.cs.cas.cz/∼strakos/strakos@cs.cas.cz

Truyen, Bart, ETRO-IRIS Research Group, Vrije Universiteit Brussel - VUB.http://www.etro.vub.ac.be/Personal/btruyen/batruyen@etro.vub.ac.be

Tyrtyshnikov, Eugene, Institute of Numerical Mathematics, Russian Academy of Sciences.http://bach.inm.ras.ru/chair.htmtee@inm.ras.ru

Van Barel, Marc, Department of Computer Science, Katholieke Universiteit Leuven, Bel-gium.http://www.cs.kuleuven.be/∼marcmarc.vanbarel@cs.kuleuven.be

Van Dooren, Paul, Department of Mathematical Engineering, Catholic University of Louvain,Belgium.http://www.inma.ucl.ac.be/∼vdooren/vdooren@csam.ucl.ac.be

Van Huffel, Sabine, Department of Electrical Engineering, Katholieke Universiteit Leuven,Belgium.http://www.esat.kuleuven.ac.be/sista/members/vanhuffel.htmlSabine.VanHuffel@esat.kuleuven.be

Vanberghen, Yvette, Department of Computer Science, Katholieke Universiteit Leuven, Bel-gium.yvette.vanberghen@cs.kuleuven.be

Vandebril, Raf, Department of Computer Science, Katholieke Universiteit Leuven, Belgium.http://www.cs.kuleuven.be/∼rafraf.vandebril@cs.kuleuven.be

Vocaturo, Alberto, Department of Matematics, University of Messina, Italy.avocat@dipmat.unime.it

Voss, Heinrich, Section of Mathematics, Technical University of Hamburg-Harburg, Ger-many.http://www.tu-harburg.de/mat/hp/voss/voss@tu-harburg.de

59

List of Speakers & Poster Presenters

Achiya Dax, 18Albert Vocaturo, 8Antonio Arico, 14

Bart Truyen, 50

Ctirad Matonoha, 8, 37

Damiana Lazzaro, 8, 34Dario Fasino, 26Diana Sima, 48

Eugene Tyrtyshnikov, 51

Heike Fassbender, 27Heinrich Voss, 55

Ivan Markovsky, 36

Jens Keiner, 32John Hakon Husøy, 31Jurjen Duintjer Tebbens, 23

Knut Hueper, 29

Lothar Reichel, 45Luca Dimiccoli, 21

Marco Donatelli, 22

Nezam Mahdavi-Amiri, 35

Ornella Menchi, 42

Patrick De Wilde, 19Paul Van Dooren, 52Per Christian Hansen, 28Peter Benner, 15

Raymond Chan, 16Rogier Ellenbroek, 8, 25

Sabine Van Huffel, 53Santa Agreste, 8, 12Serena Papi, 8, 44Stefan Kunis, 33Steven Delvaux, 20

Vadim Olshevsky, 43Vasile Sima, 46

Xiao-Wen Chang, 17

Yuli Eidelman, 24Yvette Vanberghen, 54

Zdenek Strakos, 49

60