UTV Expansion Pack: Special-purpose rank-revealing algorithms

This collection of Matlab 7.0 software supplements and complements the package UTV Tools from 1999, and includes implementations of special-purpose rank-revealing algorithms developed since the publication of the original package. We provide algorithms for computing and modifying symmetric rank-revealing VSV decompositions, we expand the algorithms for the ULLV decomposition of a matrix pair to handle interference-type problems with a rank-deficient covariance matrix, and we provide a robust and reliable Lanczos algorithm which – despite its simplicity – is able to capture all the dominant singular values of a sparse or structured matrix. These new algorithms have applications in signal processing, optimization and LSI information retrieval. AMS subject classification 65F25     

Accuracy of TSVD solutions computed from rank-revealing decompositions

Summary. Rank-revealing decompositions are favorable alternatives to the singular value decomposition (SVD) because they are faster to compute and easier to update. Although they do not yield all the information that the SVD does, they yield enough information to solve various problems because they provide accurate bases for the relevant subspaces. In this paper we consider rank-revealing decompositions in computing estimates of the truncated SVD (TSVD) solution to an overdetermined system of linear equations , where is numerically rank deficient. We derive analytical bounds which show how the accuracy of the solution is intimately connected to the quality of the subspaces. Received July 12, 1993 / Revised version received November 14, 1994     

A block algorithm for computing rank-revealing QR factorizations

We present a block algorithm for computing rank-revealing QR factorizations (RRQR factorizations) of rank-deficient matrices. The algorithm is a block generalization of the RRQR-algorithm of Foster and Chan. While the unblocked algorithm reveals the rank by peeling off small singular values one by one, our algorithm identifies groups of small singular values. In our block algorithm, we use incremental condition estimation to compute approximations to the nullvectors of the matrix. By applying another (in essence also rank-revealing) orthogonal factorization to the nullspace matrix thus created, we can then generate triangular blocks with small norm in the lower right part ofR. This scheme is applied in an iterative fashion until the rank has been revealed in the (updated) QR factorization. We show that the algorithm produces the correct solution, under very weak assumptions for the orthogonal factorization used for the nullspace matrix. We then discuss issues concerning an efficient implementation of the algorithm and present some numerical experiments. Our experiments show that the block algorithm is reliable and successfully captures several small singular values, in particular in the initial block steps. Our experiments confirm the reliability of our algorithm and show that the block algorithm greatly reduces the number of triangular solves and increases the computational granularity of the RRQR computation.This work was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, US Department of Energy, under Contract W-31-109-Eng-38. The second author was also sponsored by a travel grant from the Knud Højgaards Fond.This work was partially completed while the author was visiting the IBM Scientific Center in Heidelberg, Germany.     

An efficient rank detection procedure for modifying the ULV decomposition

The ULV decomposition (ULVD) is an important member of a class of rank-revealing two-sided orthogonal decompositions used to approximate the singular value decomposition (SVD). The problem of adding and deleting rows from the ULVD (called updating and downdating, respectively) is considered. The ULVD can be updated and downdated much faster than the SVD, hence its utility. When updating or downdating the ULVD, it is necessary to compute its numerical rank. In this paper, we propose an efficient algorithm which almost always maintains rank-revealing structure of the decomposition after an update or downdate without standard condition estimation. Moreover, we can monitor the accuracy of the information provided by the ULVD as compared to the SVD by tracking exact Frobenius norms of the two small blocks of the lower triangular factor in the decomposition.     

REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems

The package REGULARIZATION TOOLS consists of 54 Matlab routines for analysis and solution of discrete ill-posed problems, i.e., systems of linear equations whose coefficient matrix has the properties that its condition number is very large, and its singular values decay gradually to zero. Such problems typically arise in connection with discretization of Fredholm integral equations of the first kind, and similar ill-posed problems. Some form of regularization is always required in order to compute a stabilized solution to discrete ill-posed problems. The purpose of REGULARIZATION TOOLS is to provide the user with easy-to-use routines, based on numerical robust and efficient algorithms, for doing experiments with regularization of discrete ill-posed problems. By means of this package, the user can experiment with different regularization strategies, compare them, and draw conclusions from these experiments that would otherwise require a major programming effert. For discrete ill-posed problems, which are indeed difficult to treat numerically, such an approach is certainly superior to a single black-box routine. This paper describes the underlying theory gives an overview of the package; a complete manual is also available.This work was supported by grants from Augustinus Fonden, Knud Højgaards Fond, and Civ. Ing. Frants Allings Legat.     

Backward errors for eigenvalue and singular value decompositions

Summary. We present bounds on the backward errors for the symmetric eigenvalue decomposition and the singular value decomposition in the two-norm and in the Frobenius norm. Through different orthogonal decompositions of the computed eigenvectors we can define different symmetric backward errors for the eigenvalue decomposition. When the computed eigenvectors have a small residual and are close to orthonormal then all backward errors tend to be small. Consequently it does not matter how exactly a backward error is defined and how exactly residual and deviation from orthogonality are measured. Analogous results hold for the singular vectors. We indicate the effect of our error bounds on implementations for eigenvector and singular vector computation. In a more general context we prove that the distance of an appropriately scaled matrix to its orthogonal QR factor is not much larger than its distance to the closest orthogonal matrix. Received July 19, 1993     

A multiwave approximate Riemann solver for ideal MHD based on relaxation. I: theoretical framework

We present a relaxation system for ideal magnetohydrodynamics (MHD) that is an extension of the Suliciu relaxation system for the Euler equations of gas dynamics. From it one can derive approximate Riemann solvers with three or seven waves, that generalize the HLLC solver for gas dynamics. Under some subcharacteristic conditions, the solvers satisfy discrete entropy inequalities, and preserve positivity of density and internal energy. The subcharacteristic conditions are nonlinear constraints on the relaxation parameters relating them to the initial states and the intermediate states of the approximate Riemann solver itself. The 7-wave version of the solver is able to resolve exactly all material and Alfven isolated contact discontinuities. Practical considerations and numerical results will be provided in another paper.     

The QRPS algorithm: A generalization of the QR algorithm for the singular values decomposition of rectangular matrices

A generalization of the QR algorithm proposed by Francis [2] for square matrices is introduced for the singular values decomposition of arbitrary rectangular matrices. Geometrically the algorithm means the subsequent orthogonalization of the image of orthonormal bases produced in the course of the iteration. Our purpose is to show how to get a series of lower triangular matrices by alternate orthogonal-upper triangular decompositions in different dimensions and to prove the convergence of this series.     

A framework for polynomial preconditioners based on fast transforms I: Theory

Optimal and superoptimal approximations of a complex square matrix by polynomials in a normal basis matrix are considered. If the unitary transform associated with the eigenvectors of the basis matrix is computable using a fast algorithm, the approximations may be utilized for constructing preconditioners. Theorems describing how the parameters of the approximations could be efficiently computed are given, and for special cases earlier results by other authors are recovered. Also, optimal and superoptimal approximations for block matrices are determined, and the same type of theorems as for the point case are proved. This research was supported by the Swedish National Board for Industrial and Technical Development (NUTEK) and by the U.S. National Science Foundation under grant ASC-8958544.     

Analogues of Cayley graphs for topological groups

We define for a compactly generated totally disconnected locally compact group a graph, called a rough Cayley graph, that is a quasi-isometry invariant of the group. This graph carries information about the group structure in an analogous way to the ordinary Cayley graph for a finitely generated group. With this construction the machinery of geometric group theory can be applied to topological groups. This is illustrated by a study of groups where the rough Cayley graph has more than one end and a study of groups where the rough Cayley graph has polynomial growth. Supported by project J2245 of the Austrian Science Fund (FWF) and be an IEF Marie Curie Fellowship of the Commission of the European Union.     

Bounds for the least squares distance using scaled total least squares

Summary. The standard approaches to solving overdetermined linear systems construct minimal corrections to the data to make the corrected system compatible. In ordinary least squares (LS) the correction is restricted to the right hand side c, while in scaled total least squares (STLS) [14,12] corrections to both c and B are allowed, and their relative sizes are determined by a real positive parameter . As , the STLS solution approaches the LS solution. Our paper [12] analyzed fundamentals of the STLS problem. This paper presents a theoretical analysis of the relationship between the sizes of the LS and STLS corrections (called the LS and STLS distances) in terms of . We give new upper and lower bounds on the LS distance in terms of the STLS distance, compare these to existing bounds, and examine the tightness of the new bounds. This work can be applied to the analysis of iterative methods which minimize the residual norm, and the generalized minimum residual method (GMRES) [15] is used here to illustrate our theory. Received July 20, 2000 / Revised version received February 28, 2001 / Published online July 25, 2001     

Continuous crystal and Duistermaat-Heckman measure for Coxeter groups

We introduce a notion of continuous crystal analogous, for general Coxeter groups, to the combinatorial crystals introduced by Kashiwara in representation theory of Lie algebras. We explore their main properties in the case of finite Coxeter groups, where we use a generalization of the Littelmann path model to show the existence of the crystals. We introduce a remarkable measure, analogous to the Duistermaat-Heckman measure, which we interpret in terms of Brownian motion. We also show that the Littelmann path operators can be derived from simple considerations on Sturm-Liouville equations.     

Bounded-from-below solutions of the Hamilton-Jacobi equation for optimal control problems with exit times: vanishing lagrangians,eikonal equations,and shape-from-shading

We study the Hamilton-Jacobi equation for undiscounted exit time control problems with general nonnegative Lagrangians using the dynamic programming approach. We prove theorems characterizing the value function as the unique bounded-from-below viscosity solution of the Hamilton-Jacobi equation that is null on the target. The result applies to problems with the property that all trajectories satisfying a certain integral condition must stay in a bounded set. We allow problems for which the Lagrangian is not uniformly bounded below by positive constants, in which the hypotheses of the known uniqueness results for Hamilton-Jacobi equations are not satisfied. We apply our theorems to eikonal equations from geometric optics, shape-from-shading equations from image processing, and variants of the Fuller Problem.     

Iterative methods for overflow queueing models I

Summary Markovian queueing networks having overflow capacity are discussed. The Kolmogorov balance equations result in a linear homogeneous system, where the right null-vector is the steady-state probability distribution for the network. Preconditioned conjugate gradient methods are employed to find the null-vector. The preconditioner is a singular matrix which can be handled by separation of variables. The resulting preconditioned system is nonsingular. Numerical results show that the number of iterations required for convergence is roughly constant independent of the queue sizes. Analytic results are given to explain this fast convergence.     

Iterative methods for overflow queuing models II

Summary Preconditioned conjugate gradient methods are employed to find the steady-state probability distribution of Markovian queuing networks that have overflow capacity. Different singular preconditioners that can be handled by separation of variables are discussed. The resulting preconditioned systems are nonsingular. Numerical results show that the number of iterations required for convergence grows very slowly with the queue sizes.This research was supported in part by the National Science Foundation grant DCR-8405506 and DCR-8602563     

Fast Fourier-Galerkin methods for solving singular boundary integral equations: Numerical integration and precondition

We develop a fast fully discrete Fourier-Galerkin method for solving a class of singular boundary integral equations. We prove that the number of multiplications used in generating the compressed matrix is O(nlog³n), and the solution of the proposed method preserves the optimal convergence order O(n^−t), where n is the order of the Fourier basis functions used in the method and t denotes the degree of regularity of the exact solution. Moreover, we propose a preconditioning which ensures the numerical stability when solving the preconditioned linear system. Numerical examples are presented to confirm the theoretical estimates and to demonstrate the approximation accuracy and computational efficiency of the proposed algorithm.     

Convergence analysis of an accelerated expectation‐maximization algorithm for ill‐posed integral equations

The maximum‐likelihood expectation‐maximization (EM) algorithm has attracted considerable interest in single‐photon emission computed tomography, because it produces superior images in addition to be being flexible, simple, and allowing a physical interpretation. However, it often needs a large number of calculations because of the algorithm's slow rate of convergence. Therefore, there is a large body of literature concerning the EM algorithm's acceleration. One of the accelerated means is increasing an overrelaxation parameter, whereas we have not found any analysis in this method that would provide an immediate answer to the questions of the convergence. In this paper, our main focus is on the continuous version of an accelerated EM algorithm based on Lewitt and Muehllenner. We extend their conclusions to the infinite‐dimensional space and interpret and analyze the convergence of the accelerated EM algorithm. We also obtain some new properties of the modified algorithm. Copyright © 2015 John Wiley & Sons, Ltd.     

An augmented galerkin procedure for the boundary integral method applied to two-dimensional screen and crack problems

Here we present a new solution procedure for Helm-holtz and Laplacian Dirichlet screen and crack problems in IR² via boundary integral equations of the first kind having as an unknown the jump of the normal derivative across the screen or a crack curve T. Under the assumption of local finite energy we show the equivalence of the integral equations and the original boundary value problem. Via the method of local Mellin transform in [5]-[lo] and the calculus of pseudodifferential operators we derive existence, uniqueness and regularity results for the solution of our boundary integral equations together with its explicit behaviour near the screen or crack tips.With our integral equations we set up a Galerkin scheme on T and obtain high quasi-optimal convergence rates by using special singular elements besides regular splines as test and trial functions.     

Sample path Large Deviations and optimal importance sampling for stochastic volatility models

Sample path Large Deviation Principles (LDP) of the Freidlin–Wentzell type are derived for a class of diffusions, which govern the price dynamics in common stochastic volatility models from Mathematical Finance. LDP are obtained by relaxing the non-degeneracy requirement on the diffusion matrix in the standard theory of Freidlin and Wentzell. As an application, a sample path LDP is proved for the price process in the Heston stochastic volatility model.     

Principal minors, Part I: A method for computing all the principal minors of a matrix

An order O(2ⁿ) algorithm for computing all the principal minors of an arbitrary n × n complex matrix is motivated and presented, offering an improvement by a factor of n³ over direct computation. The algorithm uses recursive Schur complementation and submatrix extraction, storing the answer in a binary order. An implementation of the algorithm in MATLAB® is also given and practical considerations are discussed and treated accordingly.