Numerical computation of an analytic singular value decomposition of a matrix valued function

Summary This paper extends the singular value decomposition to a path of matricesE(t). An analytic singular value decomposition of a path of matricesE(t) is an analytic path of factorizationsE(t)=X(t)S(t)Y(t) ^T whereX(t) andY(t) are orthogonal andS(t) is diagonal. To maintain differentiability the diagonal entries ofS(t) are allowed to be either positive or negative and to appear in any order. This paper investigates existence and uniqueness of analytic SVD's and develops an algorithm for computing them. We show that a real analytic pathE(t) always admits a real analytic SVD, a full-rank, smooth pathE(t) with distinct singular values admits a smooth SVD. We derive a differential equation for the left factor, develop Euler-like and extrapolated Euler-like numerical methods for approximating an analytic SVD and prove that the Euler-like method converges.Partial support received from SFB 343, Diskrete Strukturen in der Mathematik, Universität BielefeldPartial support received from FSP Mathematisierung, Universität BielefeldPartial support received from FSP Mathematisierung, Universität BielefeldPartial support received from National Science Foundation grant CCR-8820882. Some support was also received from the University of Kansas through International Travel Fund 560478 and General Research Allocations # 3758-20-0038 and #3692-20-0038.     

Modifying the generalized singular value decomposition with application with application in direction-of-arrival finding

We consider updating and downdating problems for the generalized singular value decomposition (GSVD) of matrix pairs when new rows are added to one of the matrices or old rows are deleted. Two classes of algorithms are developed, one based on the CS decomposition formulation of the GSVD and the other based on the generalized eigenvalue decomposition formulation. In both cases we show that the updating and downdating problems can be reduced to a rank-one SVD updating problem. We also provide perturbation analysis for the cases when the added or deleted rows are subject to errors. Numerical experiments on direction-of-arrival (DOA) finding with colored noises are carried out to demonstrate the tracking ability of the algorithms. The work of the first author was supported in part by NSF grants CCR-9308399 and CCR-9619452. The work of the second author was supported in part by China State Major Key Project for Basic Researches.     

Schwarz alternating and iterative refinement methods for mixed formulations of elliptic problems,part I: Algorithms and numerical results

Summary We describe sequential and parallel algorithms based on the Schwarz alternating method for the solution of mixed finite element discretizations of elliptic problems using the Raviart-Thomas finite element spaces. These lead to symmetric indefinite linear systems and the algorithms have some similarities with the traditional block Gauss-Seidel or block Jacobi methods with overlapping blocks. The indefiniteness requires special treatment. The sub-blocks used in the algorithm correspond to problems on a coarse grid and some overlapping subdomains and is based on a similar partition used in an algorithm of Dryja and Widlund for standard elliptic problems. If there is sufficient overlap between the subdomains, the algorithm converges with a rate independent of the mesh size, the number of subdomains and discontinuities of the coefficients. Extensions of the above algorithms to the case of local grid refinement is also described. Convergence theory for these algorithms will be presented in a subsequent paper.This work was supported in part by the National Science Foundation under Grant NSF-CCR-8903003, while the author was a graduate student at New York University, and in part by the Army Research Office under Grant DAAL 03-91-G-0150, while the author was a Visiting Assistant Researcher at UCLA     

Wavelet sparse approximate inverse preconditioners

We show how to use wavelet compression ideas to improve the performance of approximate inverse preconditioners. Our main idea is to first transform the inverse of the coefficient matrix into a wavelet basis, before applying standard approximate inverse techniques. In this process, smoothness in the entries ofA ⁻¹ are converted into small wavelet coefficients, thus allowing a more efficient approximate inverse approximation. We shall justify theoretically and numerically that our approach is effective for matrices with smooth inverses. Supported by grants from ONR: ONR-N00014-92-J-1890, and the Army Research Office: DAAL-03-91-C-0047 (Univ. of Tenn. subcontract ORA4466.04 Amendment 1 and 2). The first and the third author also acknowledge support from RIACS/NASA Ames NAS 2-96027 and the Alfred P. Sloan Foundation as Doctoral Dissertation Fellows, respectively. the work was supported by the Natural Sciences and Engineering Research Council of Canada, the Information Technology Research Centre (which is funded by the Province of Ontario), and RIACS/NASA Ames NAS 2-96027.     

High-resolution boosted deconvolution of spectroscopic data

The quality of the analysis of spectrometric data consists in correct identification of the existence of peaks and subsequently in good estimation of their positions. In this paper we present high-resolution deconvolution algorithms. We propose an improvement in the efficiency by introducing a boosting operation into the deconvolution process.     

Contractivity of domain decomposition splitting methods for nonlinear parabolic problems

This work deals with the efficient numerical solution of nonlinear parabolic problems posed on a two-dimensional domain Ω. We consider a suitable decomposition of domain Ω and we construct a subordinate smooth partition of unity that we use to rewrite the original equation. Then, the combination of standard spatial discretizations with certain splitting time integrators gives rise to unconditionally contractive schemes. The efficiency of the resulting algorithms stems from the fact that the calculations required at each internal stage can be performed in parallel.     

Multilevel Schwarz methods

Summary We consider the solution of the algebraic system of equations which result from the discretization of second order elliptic equations. A class of multilevel algorithms are studied using the additive Schwarz framework. We establish that the condition number of the iteration operators are bounded independent of mesh sizes and the number of levels. This is an improvement on Dryja and Widlund's result on a multilevel additive Schwarz algorithm, as well as Bramble, Pasciak and Xu's result on the BPX algorithm. Some multiplicative variants of the multilevel methods are also considered. We establish that the energy norms of the corresponding iteration operators are bounded by a constant less than one, which is independent of the number of levels. For a proper ordering, the iteration operators correspond to the error propagation operators of certain V-cycle multigrid methods, using Gauss-Seidel and damped Jacobi methods as smoothers, respectively.This work was supported in part by the National Science Foundation under Grants NSF-CCR-8903003 at Courant Institute of Mathematical Sciences, New York University and NSF-ASC-8958544 at Department of Computer Science, University of Maryland.     

Schwarz alternating and iterative refinement methods for mixed formulations of elliptic problems,part II: Convergence theory

Summary In this paper we discuss bounds for the convergence rates of several domain decomposition algorithms to solve symmetric, indefinite linear systems arising from mixed finite element discretizations of elliptic problems. The algorithms include Schwarz methods and iterative refinement methods on locally refined grids. The implementation of Schwarz and iterative refinement algorithms have been discussed in part I. A discussion on the stability of mixed discretizations on locally refined grids is included and quantiative estimates for the convergence rates of some iterative refinement algorithms are also derived.Department of Mathematics, University of Wyoming, Laramie, WY 82071-3036. This work was supported in part by the National Science Foundation under Grant NSF-CCR-8903003, while the author was a graduate student at New York University, and in part by NSF Grant ASC 9003002, while the author was a Visiting, Assistant Researcher at UCLA.     

Additive Schwarz algorithms for parabolic convection-diffusion equations

Summary In this paper, we consider the solution of linear systems of algebraic equations that arise from parabolic finite element problems. We introduce three additive Schwarz type domain decomposition methods for general, not necessarily selfadjoint, linear, second order, parabolic partial differential equations and also study the convergence rates of these algorithms. The resulting preconditioned linear system of equations is solved by the generalized minimal residual method. Numerical results are also reported.This work was supported in part by the National Science Foundation under Grant NSF-CCR-8903003 at the Courant Institute, New York University and in part by the National Science Foundation under contract number DCR-8521451 and ECS-8957475 at Yale University     

Tensor properties of multilevel Toeplitz and related matrices

A general proposal is presented for fast algorithms for multilevel structured matrices. It is based on investigation of their tensor properties and develops the idea recently introduced by Kamm and Nagy in the block Toeplitz case. We show that tensor properties of multilevel Toeplitz matrices are related to separation of variables in the corresponding symbol, present analytical tools to study the latter, expose truncation algorithms preserving the structure, and report on some numerical results confirming advantages of the proposal.     

Domain imbedding methods for the Stokes equations

Summary We study direct and iterative domain imbedding methods for the Stokes equations on certain non-rectangular domains in two space dimensions. We analyze a continuous analog of numerical domain imbedding for bounded, smooth domains, and give an example of a simple numerical algorithm suggested by the continuous analysis. This algorithms is applicable for simply connected domains which can be covered by rectangular grids, with uniformly spaced grid lines in at least one coordinate direction. We also discuss a related FFT-based fast solver for Stokes problems with physical boundary conditions on rectangles, and present some numerical results.     

Two new variants of nonlinear inexact Uzawa algorithms for saddle-point problems

Summary. In this paper, we consider some nonlinear inexact Uzawa methods for iteratively solving linear saddle-point problems. By means of a new technique, we first give an essential improvement on the convergence results of Bramble-Paschiak-Vassilev for a known nonlinear inexact Uzawa algorithm. Then we propose two new algorithms, which can be viewed as a combination of the known nonlinear inexact Uzawa method with the classical steepest descent method and conjugate gradient method respectively. The two new algorithms converge under very practical conditions and do not require any apriori estimates on the minimal and maximal eigenvalues of the preconditioned systems involved, including the preconditioned Schur complement. Numerical results of the algorithms applied for the Stokes problem and a purely linear system of algebraic equations are presented to show the efficiency of the algorithms. Received December 8, 1999 / Revised version received September 8, 2001 / Published online March 8, 2002 RID="*" ID="*" The work of this author was partially supported by a grant from The Institute of Mathematical Sciences, CUHK RID="**" ID="**" The work of this author was partially supported by Hong Kong RGC Grants CUHK 4292/00P and CUHK 4244/01P     

Matrix computations and polynomial root-finding with preprocessing

Solution of homogeneous linear systems of equations is a basic operation of matrix computations. The customary algorithms rely on pivoting, orthogonalization and SVD, but we employ randomized preprocessing instead. This enables us to accelerate the solution dramatically, both in terms of the estimated arithmetic cost and the observed CPU time. The approach is effective in the cases of both general and structured input matrices and we extend it and its computational advantages to the solution of nonhomogeneous linear systems of equations, matrix eigen-solving, the solution of polynomial and secular equations, and approximation of a matrix by a nearby matrix that has a smaller rank or a fixed structure (e.g., of the Toeplitz or Hankel type). Our analysis and extensive experiments show the power of the presented algorithms.     

The analysis of multigrid algorithms for cell centered finite difference methods

In this paper, we examine multigrid algorithms for cell centered finite difference approximations of second order elliptic boundary value problems. The cell centered application gives rise to one of the simplest non-variational multigrid algorithms. We shall provide an analysis which guarantees that the W-cycle and variable V-cycle multigrid algorithms converge with a rate of iterative convergence which can be bounded independently of the number of multilevel spaces. In contrast, the natural variational multigrid algorithm converges much more slowly.     

Solution of augmented linear systems using orthogonal factorizations

We consider solvingx+Ay=b andA ^T x=c for givenb, c andm ×n A of rankn. This is called the augmented system formulation (ASF) of two standard optimization problems, which include as special cases the minimum 2-norm of a linear underdetermined system (b=0) and the linear least squares problem (c=0), as well as more general problems. We examine the numerical stability of methods (for the ASF) based on the QR factorization ofA, whether by Householder transformations, Givens rotations, or the modified Gram-Schmidt (MGS) algorithm, and consider methods which useQ andR, or onlyR. We discuss the meaning of stability of algorithms for the ASF in terms of stability of algorithms for the underlying optimization problems.We prove the backward stability of several methods for the ASF which useQ andR, inclusing a new one based on MGS, and also show under what circumstances they may be regarded as strongly stable. We show why previous methods usingQ from MGS were not backward stable, but illustrate that some of these methods may be acceptable-error stable. We point out that the numerical accuracy of methods that do not useQ does not depend to any significant extent on which of of the above three QR factorizations is used. We then show that the standard methods which do not useQ are not backward stable or even acceptable-error stable for the general ASF problem, and discuss how iterative refinement can be used to counteract these deficiencies.Dedicated to Carl-Eric Fröberg on the occasion of his 75th birthdayThis research was partially supported by NSERC of Canada Grant No. A9236.     

ImplicitQR factorization of a product of three matrices

We present a numerical algorithm for computing the implicit QR factorization of a product of three matrices, and we illustrate the technique by applying it to the generalized total least squares and the restricted total least squares problems. We also demonstrate how to take advantage of the block structures of the underlying matrices in order to reduce the computational work.     

Extensions of block-projections methods with relaxation parameters to inconsistent and rank-deficient least-squares problems

There exist many classes of block-projections algorithms for approximating solutions of linear least-squares problems. Generally, these methods generate sequences convergent to the minimal norm least-squares solution only for consistent problems. In the inconsistent case, which usually appears in practice because of some approximations or measurements, these sequences do no longer converge to a least-squares solution or they converge to the minimal norm solution of a “perturbed” problem. In the present paper, we overcome this difficulty by constructing extensions for almost all the above classes of block-projections methods. We prove that the sequences generated with these extensions always converge to a least-squares solution and, with a suitable initial approximation, to the minimal norm solution of the problem. Numerical experiments, described in the last section of the paper, confirm the theoretical results obtained.     

Numerical experiments with two approximate inverse preconditioners

We present the results of numerical experiments aimed at comparing two recently proposed sparse approximate inverse preconditioners from the point of view of robustness, cost, and effectiveness. Results for a standard ILU preconditioner are also included. The numerical experiments were carried out on a Cray C98 vector processor. This work was partially supported by the GA AS CR under grant 2030706 and by the grant GA CR 205/96/0921.     

A domain decomposition algorithm for elliptic problems in three dimensions

Summary Most domain decomposition algorithms have been developed for problems in two dimensions. One reason for this is the difficulty in devising a satisfactory, easy-to-implement, robust method of providing global communication of information for problems in three dimensions. Several methods that work well in two dimension do not perform satisfactorily in three dimensions.A new iterative substructuring algorithm for three dimensions is proposed. It is shown that the condition number of the resulting preconditioned problem is bounded independently of the number of subdomains and that the growth is quadratic in the logarithm of the number of degrees of freedom associated with a subdomain. The condition number is also bounded independently of the jumps in the coefficients of the differential equation between subdomains. The new algorithm also has more potential parallelism than the iterative substructuring methods previously proposed for problems in three dimensions.This work was supported in part by the National Science Foundation under grant NSF-CCR-8903003 and by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.     

Preconditioning of elliptic problems by approximation in the transform domain

Preconditioned conjugate gradient method is applied for solving linear systemsAx=b where the matrixA is the discretization matrix of second-order elliptic operators. In this paper, we consider the construction of the trnasform based preconditioner from the viewpoint of image compression. Given a smooth image, a major portion of the energy is concentrated in the low frequency regions after image transformation. We can view the matrixA as an image and construct the transform based preconditioner by using the low frequency components of the transformed matrix. It is our hope that the smooth coefficients of the given elliptic operator can be approximated well by the low-rank matrix. Numerical results are reported to show the effectiveness of the preconditioning strategy. Some theoretical results about the properties of our proposed preconditioners and the condition number of the preconditioned matrices are discussed.