Second Workshop – Rennes, February 14-15, 2002

Session : Correlation Analysis and Independent Component Analysis
Nicholas Higham. Computing the Nearest Correlation Matrix - A Problem from Finance.  Abstract  Slides
Eric Moreau. Independent Component Analysis (Sources Separation) and Diagonalization.  Abstract. Slides
Mohamed Hanafi. Global optimality of the successive MAXBET algorithm. Abstract Slides

Session : Triangular and QR factorisation
Per Christian Hansen. Computing Symmetric Rank-Revealing Decompositions via Triangular Factorization.  Abstract.  Slides
Patrick R. Amestoy, Symmetric minimum priority orderings for sparse unsymmetric factorization.  Abstract. Slides
Cristian Gatu. Algorithms for computing the best subset regression models using the QR decomposition.  Abstract. Slides
Petko Yanev. Computing the QR decomposition of a set of matrices with common columns.  Abstract.  Slides

 Session : Regression models
Rasmus Bro. Jack-knife for estimation of standard errors and outlier detection in E.G. PARAFAC models. Abstract. Slides
Giorgio Tomasi. Fitting the PARAFAC model.  Abstract. Slides
Paolo Foschi. A comparative study of algorithms for solving seemingly unrelated regressions models.   Abstract. Slides

SVD and Least-squares computations
Lars Elde’n, PCR and PLS vs. TSVD and LSQR. AbstractSlides References
George Moustakides. SVD methods for CDMA communications. Abstract. Slides
Yann-Hervé De Roeck. Regularization of large least square problems in geophysics. Abstract.  Slides
Jocelyne Erhel. Some inverse problems in electrocardiography. Abstract.Slides
Andy Nisbet. Hardware Acceleration of Applications Using FPGAs. AbstractSlides


Computing the Nearest Correlation Matrix - A Problem from Finance

Nicholas Higham
University of Manchester
Department of Mathematics
Manchester, M13 9PL, England

Given a symmetric matrix what is the nearest correlation matrix, that is, the nearest symmetric positive semidefinite matrix with unit diagonal? This problem arises in the finance industry, where the correlations are between stocks. For distance measured in a weighted Frobenius norm we characterize the solution using convex analysis. We show how the modified alternating projections method can be used to compute the solution. In the finance application the original matrix has many zero or negative eigenvalues; we show that the nearest correlation matrix has correspondingly many zero eigenvalues and that this fact can be exploited in the computation.

Independent Component Analysis (Sources Separation) and  Diagonalization

Eric Moreau

SIS-SD, ISITV, Univ. Toulon
Av. G. Pompidou, BP 56
F-83162 La Valette du Var Cedex, France

In recent years, the independent component analysis problem and more specifically the sources separation problem have received a lot of attentions in the signal processing community because of their potential practical interest. Numerous solutions have been proposed. Among them one can find matrices (or tensors) joint-diagonalization based algorithms.
We first review the sources separation problem. Then we present the principle of many matrices joint-diagonalization algorithms in an unifying point of view. We also show some links of joint-diagonalization criteria with independence measure called contrast functions. We present some extensions to the joint-diagonalization of tensors. Finally some perspectives are given together with references.

Global optimality of the successive MAXBET algorithm

Mohamed Hanafi (1)  and  Jos M. F. Ten Berge (2)

(1) SMAD, ENITIAA, Domaine de la Géraudière, BP 82225, 44322 Nantes Cedex 3,
(2) Department of Psychology, University of Groningen, Grote Kruisstraat 2/1, 9712 TS Groningen, Netherlands,

The Maxbet method is an alternative to the method of generalized canonical correlation analysis and of generalized Procrustes analysis. Contrary to these methods, it does not maximize the inner products (covariances) between linear composites, but also takes their sums of squares (variances) into account. It is well-known that the Maxbet algorithm, which has been proven to converge monotonically, may converge to local maxima. The present paper discusses an eigenvalue criterion which is sufficient, but not necessary for global optimality. However, in two special cases, the eigenvalue criterion is shown to be necessary and sufficient for global optimality. The first case is when there are only two data sets involved; the second case is when the inner products between all variables involved are positive, regardless of the number of data sets.

Computing Symmetric Rank-Revealing Decompositions via Triangular Factorization

Per Christian Hansen

Informatics and Mathematical Modelling
Building 321, Technical University of Denmark
DK-2800 Lyngby, Denmark

We present a family of algorithms for computing symmetric rank-revealing VSV decompositions, based on triangular factorization of the matrix. The VSV decomposition consists of a middle symmetric matrix that reveals the numerical rank in having three blocks with small norm, plus an orthogonal matrix whose columns span approximations to the numerical range and null space. We show that for semi-definite matrices the VSV decomposition should be computed via the ULV decomposition, while for indefinite matrices it must be computed via a URV-like decomposition that involves hyperbolic rotations. Truncated VSV solutions are obtained by neglecting the small-norm blocks in the middle matrix. We compare these solutions with truncated SVD solutions, and give perturbation bounds for VSV solutions. Keywords: rank-revealing decompositions, matrix approximation, symmetric matrices, perturbation bounds

Symmetric minimum priority orderings for sparse unsymmetric factorization

P.R. Amestoy (1), X.S. Li, and E.G. Ng (2)

(1) LIMA-IRIT, Toulouse
(2) Lawrence Berkeley Lab.

We introduce minimum priority ordering algorithms forsparse LU factorization where we assume that pivots can be chosen on the diagonal. All approaches based on a separate preordering phase traditionally involve the structure of A'+A to compute the sparsity ordering. But forming A'+A loses some information of A itself. Our ordering algorithms are based solely on the (unsymmetric) structure of A. Its efficient implementation requires us to extend the (symmetric) quotient graph model to the bipartite quotient graph. We will present the experimental results obtained on both multifrontal and supernodal approaches  (MUMPS from Lyon-Toulouse and SuperLU from Berkeley) to show that with the new orderings we can better preserve the sparsity and exploit the asymmetry of the matrix.

Algorithms for computing the best subset regression models using the QR decomposition

Cristian Gatu  and  Erricos John Kontoghiorghes

Institut d'informatique, Univ. of Neuchatel, Switzerland

Efficient parallel algorithms for computing all possible subset regression models are proposed. The algorithms are based on the dropping columns method that generates a regression tree. The properties of the tree are exploited in order to provide an efficient load balancing which results in no inter-processor communication. Theoretical measures of complexity suggest linear speed-up. The parallel algorithms are extended to deal with the general linear and seemingly unrelated regression models.  Experimental results on a shared memory machine are presented and analyzed.
A branch and bound algorithm which obtains the best models without investigating the whole regression tree that generates all possible subset models is proposed. The algorithm has the QR decomposition as a main computational tool and outperforms an existing branch and bound strategy based on inverse matrices. Computational results demonstrate the efficiency of the new algorithm.

Computing the QR decomposition of a set of matrices with common columns

Petko Yanev, Paolo Foschi and Erricos Kontoghiorghes

Institut d'informatique, Univ. of Neuchatel, Switzerland

The QR decomposition of a set of matrices which have common columns is investigated.  The triangular factors of the QR decompositions are represented as nodes of a directed graph, where edges have weights. An edge between two nodes exists if and only if the columns of one of the matrices is a subset of the other.  The weight of an edge denotes the computational complexity of deriving the triangular factor of the destination matrix from the triangular factor of the source node. Within this context the problem becomes equivalent to construct the graph and finding the minimum cost for visiting all the nodes.  An algorithm which computes the QR decompositions by deriving the minimum spanning arborescence of the graph is proposed. Theoretical measures of complexity and numerical results from the implementation of this and alternative algorithms are derived and analysed.

Jack-knife for estimation of standard errors and outlier detection in E.G. PARAFAC models

Rasmus Bro and Jordi Riu
LMT, MLI, KVL, Rolighedsvej 30, DK-1958 Frederiksberg C
Denmark Phone +45 3528 3296 - Fax   +45 3528 3505 .

In the last years multi-way analysis has become increasingly important because it has proved to be a valuable tool e.g. in interpreting data provided by current instrumental methods. PARAFAC is one of the most widely used multi-way models. But despite its usefulness in some real examples, up to date there is no available tool in the literature to estimate the standard errors associated with the parameter estimates. In this study we apply the so-called jack-knife technique to PARAFAC in order to find the associated standard errors to the parameter estimates from the PARAFAC model. The jack-knife technique is also shown to be useful for detecting outliers. A real example corresponding to fluorescence data (emission/excitation landscapes) is used to show the applicability of the method.

Fitting the PARAFAC model

Giorgio Tomasi

Rolighedsvej 30
1958 Frederiksberg
Tel: +45 35 28 35 01
Mail: /

A multitude of different algorithms have been developed to fit a trilinear model to a three-way array. The aim of this study is to compare some of the available methods and to point out their limits and advantages. The algorithms considered in this paper are: GRAM-DTLD, PARAFAC-ALS, ASD, SWATLD, dGN (also referred to as Levenberg-Marquadt) and PMF3. These algorithms can be classified according to different criteria. One such classification is based on the "general" method they use to calculate the solution: GRAM is the only direct method, PARAFAC-ALS, ASD and SWATLD are based on the alternating least squares idea while dGN and PMF3 are derived from the Gauss-Newton method for the least squares problems. The algorithms and the underlying theories are exposed in general terms and two accessorial ethods, i.e. line search and compression, are illustrated when applied to the multiway case. In order to compare the algorithms, 720 sets of artificial data were generated varying features such as level and type of noise, collinearity of the factors and rank, and two PARAFAC models were fitted on each data set using the different methods: the first fitting the correct number of factors F and the second with F+1 components (the objective being to assess the sensibility of the different approaches to the over-factoring problem, i.e. when the number of extracted components exceeds the rank of the array). The algorithms have also been tested on two different real data sets constituted by fluorescence measurements again by extracting both the right and an exceeding number of factors. The evaluations are based on the number of iterations necessary to reach convergence, the number of operations per iteration, the time consumption, the quality of the solution and the amount of resources required for the calculations (mainly memory requirement) even though the last might be considered secondary.

PCR and PLS vs. TSVD and LSQR

Lars Elde’n

Department of Mathematics,
Linköpings University
SE-581 83 Linköpings

There are several connections between certain methods for regression analysis and numerical algorithms for solving ill-conditioned least squares problems. In fact, principal component regression, PCR, is the same as truncated singular value expansion, TSVD, and partial least squares, PLS, is the same as Lanczos bidiagonalization (one implementation of which is called LSQR). By considering so called filter factors, which relate PLS to PCR, we demonstrate how PLS is capable of quickly reducing the least squares residual at the same time as an approximation of the data matrix is obtained. We will also show that similar residual reduction can be obtained using principal components taken in an order determined by the right hand side.

SVD methods for CDMA communications

George Moustakides

Campus de Beaulieu
35042 RENNES Cedex - FRANCE
email :

Code Division Multiple Access (CDMA) is a technology applied to a number of important communications applications as mobile telephony, wireless networks and personal communciations. We present the theoretical modeling of this popular communication system, focusing on its main characteristic which consists in the simultaneous usage of the same communication channel by several users. This important advantage however, tends also to be its main handicap since it is the principal source of performance degradation. Indeed, it is clear that, to each user all other users play the role of interference. We therefore present popular signal processing algorithms that are used to combat multiuser interference. In particular we focus on adaptive methods, that rely on subspace tracking, which have a significantly better performance than their conventional least squares counterparts.

Regularization of large least square problems in geophysics

Yann-Hervé De Roeck, Sandie Le Conte

Acoustic and seismic team,
Ifremer, Brest, France.
BP 70
F-29280 Plouzane
Email :
+33 (0)2 98 22 41 43

The processing of reflection seismic data in geophysics requires some very CPU-intensive computations. Prestack depth migration can namely be viewed, in the framework of waveform inversion, as a very large least square problem to be solved with much care about regularization issues. Although iterative methods are still preferred for their matrix-free implementations, direct methods might be tailored for this purpose. The application of both types of methods has been reviewed in our recent studies. A special attention is also paid to the deconvolution algorithm, as it provides a clear view about the difficult task of setting a threshold for the regularization parameters that should reach an optimal value depending on the noise level in the data.

Some inverse problems in electrocardiography.

Jocelyne Erhel and Edouard Canot

Campus de Beaulieu
35042 RENNES Cedex - FRANCE
email :,

The electrical activity generated by the heart propagates through the chest. Direct problems in electrocardiography aims at computing the potential at the surface of the torso, using a model of the electrical cardiac activity. For inverse problems, we consider a model where the cardiac activity is represented by the potential at the surface of the heart. The solution is then unique but not continuous. After a discretisation using either boundary element methods or finite element methods, we are faced to a discrete ill-posed problem. We formulate it as a general least-squares problem, which has to be regularized. We consider also time regularization, when a sequence of problems has to be solved.

Hardware Acceleration of Applications Using FPGAs

Andy Nisbet

ORI.LG.19, Dept. of Computer Science,
Trinity College, Dublin 2, Ireland
Phone: +353-(01)-608-3682      FAX: +353-(01)-677-2204

High-level programming languages based on C, such as HandelC/SystemC can be used to program FPGA reconfigurable logic devices. In this talk we will present an overview of HandelC and discuss how it can be used to accelerate the execution of matrix and statistic applications. Single-precision floating-point arithmetic can be handled as fixed-point arithmetic, or by using the logarithmic arithmetic unit developed under the European Logarithmic Microprocessor, High Speed Logarithmic Arithmetic EU Long Term Research Project 33544. Current work at Trinity College has started into automatic conversion of C applications into HandelC/SystemC and  the acceleration of lattice QCD and image segmentation applications using FPGAs.

Tree-based methods and generalized additive multi-mixture models for data mining.

 Claudio Conversano and Roberta Siciliano

Dept. of Mathematics & Statistics
University of Naples Federico II

In the framework of supervised learning, a methodology based on a combination of classification/prediction procedures is presented as a data mining tool enabling to extract information from complex and huge data sets for statistical purposes. This approach consists in a three stages data-driven procedure: tree-partitioning, modeling and model fusion. Main tools are represented by the tree production rules and nonlinear regression models. Particularly, in the first stage a recursive partitioning procedure is applied in order to find a suitable partitioning of the raw data into a relatively small number of internal homogeneous subgroups on which standard statistical modelling can be applied. Then, in each subgroup the dependence relation is described on the basis of the Generalized Additive-Multi-Mixture Models. This methodology, introduced by Conversano et al. (see reference), derives from a suitable combinations of different classifiers for statistical learning purposes. It results in a regression-like model defined upon a multiple combination of mixtures of models associated to the predictors, and has proven to be effective and feasible to capture the non-linearity structure of the data. At the end, in order to combine the results of the two previous stages, a final decision rules is formulated enabling to classify, or to predict, new instances. This methodology can be fruitfully applied for many data mining applications. Particularly, we focalize our attention to the analysis financial data, with particular emphasis on the identification of a productive asset allocation strategy in the framework of the management of financial portfolios.

Main References

Conversano, C. (2002). Regression Tree Criteria for Smoothing and Variable Selection in Generalized Additive Models: Modeling MIB30 Index Returns, Technical Report, Dept. of Mathematics and Statistics.

Conversano, C. (2002). Bagged Mixtures of Classifiers using Model Scoring Criteria, Technical Report, Dept. of Mathematics and Statistics.

Conversano, C., Mola, F. and Siciliano, R. (2001). Partitioning Algorithms and Combined Model Integration for Data Mining, Computational Statistics, 16: 323-339.

Conversano, C., Siciliano, R. and Mola, F. (2002). Generalized Additive Multi-Mixture Model for Data Mining, Computational Statistics and Data Analysis, 38/4:487-500.

A comparative study of algorithms for solving seemingly unrelated regressions models

Paolo Foschi (1), Dave Belsley (2) and Erricos J. Kontoghiorghes (1)

1. Institut d'informatique, Univ. of Neuchatel, Switzerland
2. Department of Economics, Boston College, USA

The computational efficiency of various algorithms for solving Seemingly Unrelated Regressions (SUR) models is investigated.  Some of the algorithms adapt known methods; others are new. The first transforms the SUR to an ordinary linear model and uses the QR decomposition to solve it. Three others employ the generalized QR decomposition to solve the SUR model formulated as a generalized linear least squares problem.  Strategies to exploit the structure of the matrices involved are developed.  The algorithms are reconsidered for the solution of the SUR model after it has been transformed to one of smaller dimensions.