Generalized Schur Algorithm Research Articles

Ab initio uncertainty quantification in scattering analysis of microscopy.

Estimating parameters from data is a fundamental problem in physics, customarily done by minimizing a loss function between a model and observed statistics. In scattering-based analysis, it is common to work in the reciprocal space. Researchers often employ their domain expertise to select a specific range of wave vectors for analysis, a choice that can vary depending on the specific case. We introduce another paradigm that defines a probabilistic generative model from the beginning of data processing and propagates the uncertainty for parameter estimation, termed the ab initio uncertainty quantification (AIUQ). As an illustrative example, we demonstrate this approach with differential dynamic microscopy (DDM) that extracts dynamical information through minimizing a loss function for the squared differences of the Fourier-transformed intensities, at a selected range of wave vectors. We first show that the conventional way of estimation in DDM is equivalent to fitting a temporal variogram in the reciprocal space using a latent factor model as the generative model. Then we derive the maximum marginal likelihood estimator, which optimally weighs the information at all wave vectors, therefore eliminating the need to select the range of wave vectors. Furthermore, we substantially reduce the computational cost of computing the likelihood function without approximation, by utilizing the generalized Schur algorithm for Toeplitz covariances. Simulated studies of a wide range of dynamical systems validate that the AIUQ method improves estimation accuracy and enables model selection with automated analysis. The utility of AIUQ is also demonstrated by three distinct sets of experiments: first in an isotropic Newtonian fluid, pushing limits of optically dense systems compared to multiple particle tracking; next in a system undergoing a sol-gel transition, automating the determination of gelling points and critical exponent; and lastly, in discerning anisotropic diffusive behavior of colloids in a liquid crystal. These studies demonstrate that the new approach does not require manually selecting the wave vector range and enables automated analysis.

The Generalized Schur Algorithm and Some Applications

The generalized Schur algorithm is a powerful tool allowing to compute classical decompositions of matrices, such as the Q R and L U factorizations. When applied to matrices with particular structures, the generalized Schur algorithm computes these factorizations with a complexity of one order of magnitude less than that of classical algorithms based on Householder or elementary transformations. In this manuscript, we describe the main features of the generalized Schur algorithm. We show that it helps to prove some theoretical properties of the R factor of the Q R factorization of some structured matrices, such as symmetric positive definite Toeplitz and Sylvester matrices, that can hardly be proven using classical linear algebra tools. Moreover, we propose a fast implementation of the generalized Schur algorithm for computing the rank of Sylvester matrices, arising in a number of applications. Finally, we propose a generalized Schur based algorithm for computing the null-space of polynomial matrices.

A Divide and Conquer Algorithm for the Superfast Solution of Toeplitz-like Systems

In this paper a new $O(N \,\log^3N)$ solver for $N\times N$ Toeplitz-like systems, based on a divide and conquer technique, is presented. Similarly to the superfast algorithm MBA for the inversion of a Toeplitz-like matrix [R. R. Bitmead and B. D. O. Anderson, Linear Algebra Appl., 34 (1980), pp. 103--116; M. Morf, Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, 1980, pp. 954--959], it exploits the displacement properties. In order to avoid the well-known numerical instability of the explicit inversion, the new algorithm relies on the triangular factorization and back-substitution formula for the system seen as a $2\times 2$ block system with blocks of half size. This idea is the one used in [M. Stewart, SIAM J. Matrix Anal. Appl., 25 (2003), pp. 669--693] to improve the numerical stability of superfast methods based on the generalized Schur algorithm for positive definite Toeplitz matrices, but the algorithm we propose can be applied also to nonsymmetric Toeplitz-like systems. The stability of the algorithm is examined through numerical experiments.

A fast solver for linear systems with displacement structure

We describe a fast solver for linear systems with reconstructable Cauchy-like structure, which requires O(rn^2) floating point operations and O(rn) memory locations, where n is the size of the matrix and r its displacement rank. The solver is based on the application of the generalized Schur algorithm to a suitable augmented matrix, under some assumptions on the knots of the Cauchy-like matrix. It includes various pivoting strategies, already discussed in the literature, and a new algorithm, which only requires reconstructability. We have developed a software package, written in Matlab and C-MEX, which provides a robust implementation of the above method. Our package also includes solvers for Toeplitz(+Hankel)-like and Vandermonde-like linear systems, as these structures can be reduced to Cauchy-like by fast and stable transforms. Numerical experiments demonstrate the effectiveness of the software.

A fast algorithm for computing the smallest eigenvalue of a symmetric positive‐definite Toeplitz matrix

AbstractRecent progress in signal processing and estimation has generated considerable interest in the problem of computing the smallest eigenvalue of a symmetric positive‐definite (SPD) Toeplitz matrix.An algorithm for computing upper and lower bounds to the smallest eigenvalue of a SPD Toeplitz matrix has been recently derived (Linear Algebra Appl. 2007; DOI: 10.1016/j.laa.2007.05.008). The algorithm relies on the computation of the R factor of the QR factorization of the Toeplitz matrix and the inverse of R. The simultaneous computation of R and R−1 is efficiently accomplished by the generalized Schur algorithm. In this paper, exploiting the properties of the latter algorithm, a numerical method to compute the smallest eigenvalue and the corresponding eigenvector of SPD Toeplitz matrices in an accurate way is proposed. Copyright © 2008 John Wiley & Sons, Ltd.

A Schur-based algorithm for computing bounds to the smallest eigenvalue of a symmetric positive definite Toeplitz matrix

Recent progress in signal processing and estimation has generated c onsiderable interest in the problem of computing the smallest eigenvalue of symmetric positive definite Toeplitz matrices. Several algorithms have been proposed in the literature. They compute the smallest eigenvalue in an iterative fashion, many of them relying on the Levinson–Durbin solution of sequences of Yule–Walker systems. Exploiting the properties of two algorithms recently developed for estimating a lower and an upper bound of the smallest singular value of upper triangular matrices, respectively, an algorithm for computing bounds to the smallest eigenvalue of a symmetric positive definite Toeplitz matrix is derived. The algorithm relies on the computation of the R factor of the QR-factorization of the Toeplitz matrix and the inverse of R. The simultaneous computation of R and R −1 is efficiently accomplished by the generalized Schur algorithm.

A fast algorithm for solving the Sylvester structured total least squares problem

In this paper, we develop a fast structured total least squares (STLS) algorithm for computing an approximate greatest common divisor (GCD) of two univariate polynomials. By exploiting the displacement structure of the Sylvester matrix and applying the generalized Schur algorithm, each single iteration of the proposed algorithm has quadratic computational complexity in the degrees of the given polynomials.

A structured rank-revealing method for Sylvester matrix

We propose a fast algorithm for computing the numeric ranks of Sylvester matrices. Let S denote the Sylvester matrix and H denote the Hankel-like-Sylvester matrix. The algorithm is based on a fast Cholesky factorization of S T S or H T H and relies on a stabilized version of the generalized Schur algorithm for matrices with displacement structure. All computations can be done in O ( r ( n + m ) ) , where n + m and r denote the size and the numerical rank of the Sylvester matrix, respectively.

A generalized isometric Arnoldi algorithm

This paper describes a generalization of the isometric Arnoldi algorithm and shows that it can be interpreted as a structured form of modified Gram–Schmidt. Given an isometry A, the algorithm efficiently orthogonalizes the columns of a sequence of matrices M j for j ⩾ 0 (with M −1 = 0) for which the columns of M j − AM j−1 are in a fixed finite dimensional subspace for each j ⩾ 0. The dimension of the subspace is analogous to displacement rank in the generalized Schur algorithm. The algorithm is described in terms of projections and inner products. This is in contrast to orthogonalization methods based on the generalized Schur algorithm, for which Cholesky factorization is central to the computation. Numerical experiments suggest that, relative to a generalized Schur algorithm, the new algorithm improves the numerical orthogonality of the computed orthonormal sequence.

An Efficient Parallel Algorithm to Solve Block?Toeplitz Systems

In this paper, we present an efficient parallel algorithm to solve Toeplitz–block and block–Toeplitz systems in distributed memory multicomputers. This algorithm parallelizes the Generalized Schur Algorithm to obtain the semi-normal equations. Our parallel implementation reduces the communication cost and optimizes the memory access. The experimental analysis on a cluster of personal computers shows the scalability of the implementation. The algorithm is portable because it is based on standard tools and libraries, such as ScaLAPACK and MPI.

Fast direct solvers for some complex symmetric block Toeplitz linear systems

We consider the solution of a class of complex symmetric block Toeplitz linear systems, arising from integral equation problems. Algorithms that exploit the Toeplitz structure provide considerable savings on the number of arithmetic operations, compared to the classical Cholesky factorization. We propose a generalized Schur algorithm adapted to the complex symmetric case. We detail blocked variants, that perform better by using BLAS 3 primitives. We also propose a solver, based on an augmented matrix approach, that allows a substantial decrease in the use of memory, by avoiding an explicit assembly of the Cholesky factor. All algorithms have been implemented and numerical results are included to illustrate the effectiveness of our approach.

Fast regularized structured total least squares algorithm for solving the basic deconvolution problem

AbstractIn this paper, we develop a fast regularized structured total least squares (RSTLS) algorithm for solving the basic deconvolution problem. The algorithm is based on a particular implementation of the generalized Schur algorithm. We apply the new algorithm on a deconvolution problem arising in a medical application in renography. By means of this example, we show that regularization in the structured total least squares framework yields more accurate results compared to other estimators. Copyright © 2004 John Wiley & Sons, Ltd.

Implementation of the regularized structured total least squares algorithms for blind image deblurring

The structured total least squares (STLS) problem has been introduced to handle problems involving structured matrices corrupted by noise. Often the problem is ill-posed. Recently, regularization has been proposed in the STLS framework to solve ill-posed blind deconvolution problems encountered in image deblurring when both the image and the blurring function have uncertainty. The kernel of the regularized STLS (RSTLS) problem is a least squares problem involving Block–Toeplitz–Toeplitz–Block matrices. In this paper an algorithm is described to solve this problem, based on a particular implementation of the generalized Schur Algorithm (GSA). It is shown that this new implementation improves the computational efficiency of the straightforward implementation of GSA from O( N 2.5) to O( N 2), where N is the number of pixels in the image.

Incremental refinement of a multi-user-detection algorithm (II)

Abstract. Multi-user detection is a technique proposed for mobile radio systems based on the CDMA principle, such as the upcoming UMTS. While offering an elegant solution to problems such as intra-cell interference, it demands very significant computational resources. In this paper, we present a high-level approach for reducing the required resources for performing multi-user detection in a 3GPP TDD multi-user system. This approach is based on a displacement representation of the parameters that describe the transmission system, and a generalized Schur algorithm that works on this representation. The Schur algorithm naturally leads to a highly parallel hardware implementation using CORDIC cells. It is shown that this hardware architecture can also be used to compute the initial displacement representation. It is very beneficial to introduce incremental refinement structures into the solution process, both at the algorithmic level and in the individual cells of the hardware architecture. We detail these approximations and present simulation results that confirm their effectiveness.

A fast algorithm for subspace state-space system identification via exploitation of the displacement structure

Two recent approaches (Van Overschee, De Moor, N4SID, Automatica 30 (1) (1994) 75; Verhaegen, Int. J. Control 58(3) (1993) 555) in subspace identification problems require the computation of the R factor of the QR factorization of a block-Hankel matrix H, which, in general has a huge number of rows. Since the data are perturbed by noise, the involved matrix H is, in general, full rank. It is well known that, from a theoretical point of view, the R factor of the QR factorization of H is equivalent to the Cholesky factor of the correlation matrix H T H, apart from a multiplication by a sign matrix. In Sima (Proceedings Second NICONET Workshop, Paris-Versailles, December 3, 1999, p. 75), a fast Cholesky factorization of the correlation matrix, exploiting the block-Hankel structure of H, is described. In this paper we consider a fast algorithm to compute the R factor based on the generalized Schur algorithm. The proposed algorithm allows to handle the rank–deficient case.

On a class of matrices with low displacement rank

A matrix A such that, for some matrices U and V, the matrix AU−VA or the matrix A−VAU has a rank which is small compared with the order of the matrix is called a matrix with displacement structure. In this paper the authors single out a new class of matrices with displacement structure—namely, finite sections of recursive matrices—which includes the class of finite Hurwitz matrices, widely used in computer graphics. For such matrices it is possible to give an explicit evaluation of the displacement rank and an algebraic description, and hence a stable numerical evaluation, of the corresponding generators. The generalized Schur algorithm can therefore be used for the fast factorization of some classes of such matrices, and exhibits a backward stable behavior.

Schur-Type Methods for Solving Least Squares Problems with Toeplitz Structure

We give an overview of fast algorithms for solving least squares problems with Toeplitz structure, based on generalization of the classical Schur algorithm, and discuss their stability properties. In order to obtain more accurate triangular factors of a Toeplitz matrix as well as accurate solutions for the least squares problems, methods based on corrected seminormal equations (CSNE) can be used. We show that the applicability of the generalized Schur algorithm is considerably enhanced when the algorithm is used in conjunction with CSNE. Several numerical tests are reported, where different variants of the generalized Schur algorithm and CSNE are compared for their accuracy and speed.

Computation of coprime factorizations of rational matrices

We propose numerically reliable state-space algorithms for computing several coprime factorizations of rational matrices: (1) factorizations with factors having poles in a given stability domain; (2) factorizations with proper stable factors; (3) factorizations with inner and J-inner denominators. The new algorithms are based on a recursive generalized Schur algorithm for pole dislocation. They are generally applicable whether the underlying descriptor state-space representation is minimal or not, and whether it is stabilizable/detectable or not. The proposed algorithms are useful in solving various computational problems for both standard and descriptor system representations.

Computation of coprime factorizations of rational matrices

We propose numerically reliable state-space algorithms for computing several coprime factorizations of rational matrices: (1) factorizations with factors having poles in a given stability domain; (2) factorizations with proper stable factors; (3) factorizations with inner and J-inner denominators. The new algorithms are based on a recursive generalized Schur algorithm for pole dislocation. They are generally applicable whether the underlying descriptor state-space representation is minimal or not, and whether it is stabilizable/detectable or not. The proposed algorithms are useful in solving various computational problems for both standard and descriptor system representations.

A Fast Stable Solver for Nonsymmetric Toeplitz and Quasi-Toeplitz Systems of Linear Equations

We derive a stable and fast solver for nonsymmetric linear systems of equations with shift structured coefficient matrices (e.g., Toeplitz, quasi-Toeplitz, and product of two Toeplitz matrices). The algorithm is based on a modified fast QR factorization of the coefficient matrix and relies on a stabilized version of the generalized Schur algorithm for matrices with displacement structure. All computations can be done in O(n2 ) operations, where n is the matrix dimension, and the algorithm is backward stable.