Search

A framework, PhaseLift, was recently proposed to solve the phase retrieval problem. In this framework, the problem is solved by optimizing a cost function over the set of complex Hermitian positive semidefinite matrices. This approach to phase retrieval motivates a more general consideration of optimizing cost functions on semidefinite Hermitian matrices where the desired minimizers are known to have low rank. This paper considers an approach based on an alternative cost function defined on a union of appropriate manifolds. It is related to the original cost function in a manner that preserves the ability to find a global minimizer and is significantly more efficient computationally. A rank-based optimality condition for stationary points is given and optimization algorithms based on state-of-the-art Riemannian optimization and dynamically reducing rank are proposed. Empirical evaluations are performed using the PhaseLift problem. The new approach is shown to be an effective method of phase retrieval with computational efficiency increased substantially compared to a state-of-the-art algorithm, the Wirtinger flow algorithm. A preliminary version of this paper can be found in [W. Huang, K. A. Gallivan, and X. Zhang, Procedia Comput. Sci., 80 (2016), pp. 1125-1134]. [ABSTRACT FROM AUTHOR]

A systematic mathematical framework for the study of numerical algorithms would allow comparisons, facilitate conjugacy arguments, and enable the discovery of improved, accelerated, data-driven algorithms. Over the course of the past century, the Koopman operator has provided a mathematical framework for the study of dynamical systems which facilitates conjugacy arguments and can provide efficient reduced descriptions. More recently, numerical approximations of the operator have enabled the analysis of a large number of deterministic and stochastic dynamical systems in a completely data-driven, essentially equation-free pipeline. Discrete- or continuous-time numerical algorithms (integrators, nonlinear equation solvers, optimization algorithms are themselves dynamical systems. In this paper, we use this insight to leverage the Koopman operator framework in the data-driven study of such algorithms and discuss benefits for analysis and acceleration of numerical computation. For algorithms acting on high-dimensional spaces by quickly contracting them toward low-dimensional manifolds, we demonstrate how basis functions adapted to the data help to construct efficient reduced representations of the operator. Our illustrative examples include the gradient descent and Nesterov optimization algorithms as well as the Newton-Raphson algorithm. [ABSTRACT FROM AUTHOR]

We introduce a new algorithm of proper generalized decomposition (PGD) for parametric symmetric elliptic partial differential equations. For any given dimension, we prove the existence of an optimal subspace of at most that dimension which realizes the best approximation---in the mean parametric norm associated to the elliptic operator---of the error between the exact solution and the Galerkin solution calculated on the subspace. This is analogous to the best approximation property of the proper orthogonal decomposition (POD) subspaces, except that in our case the norm is parameter-dependent. We apply a deflation technique to build a series of approximating solutions on finite-dimensional optimal subspaces, directly in the online step, and we prove that the partial sums converge to the continuous solution in the mean parametric elliptic norm. We show that the standard PGD for the considered parametric problem is strongly related to the deflation algorithm introduced in this paper. This opens the possibility of computing the PGD expansion by directly solving the optimization problems that yield the optimal subspaces. [ABSTRACT FROM AUTHOR]

The Euler genus of a graph is a fundamental and well-studied parameter in graph theory and topology. Computing it has been shown to be NP-hard by Thomassen [J. Algorithms, 10 (1989), pp. 568{576; J. Combin. Theory, Ser. B, 57 (1993), pp. 196{206], and it is known to be fixed-parameter tractable. However, the approximability of the Euler genus is wide open. While the existence of an O(1)-approximation is not ruled out, only an O( p n)-approximation [J. Chen, S. P. Kanchi, and A. Kanevsky, Inform. Process. Lett., 61 (1997), pp. 317{322] is known even in bounded-degree graphs. In this paper we give a polynomial-time algorithm which, given a boundeddegree graph of Euler genus g, computes a drawing in a surface of Euler genus gO(1) · logO(1) n. Combined with the upper bound from [J. Chen, S. P. Kanchi, and A. Kanevsky, Inform. Process. Lett., 61 (1997), pp. 317{322], our result also implies a O(n1/2-α)-approximation for some constant α > 0. Using our algorithm for approximating the Euler genus as a subroutine, we obtain, in a uniform fashion, algorithms with approximation ratios of the form OPTO(1) · logO(1) n for several related problems on bounded-degree graphs. These include the problems of orientable genus, crossing number, and planar edge and vertex deletion. Our algorithm and proof of correctness for the crossing number problem are simpler compared to the long and difficult proof in the recent breakthrough by Chuzhoy [Proceedings of the ACM Symposium on Theory of Computing, 2011, pp. 303{312], while essentially obtaining a qualitatively similar result. For planar edge and vertex deletion problems our results are the first to obtain a bound of the form poly(OPT; log n). We also highlight some further applications of our results in the design of algorithms for graphs with small genus. Many such algorithms require that a drawing of the graph is given as part of the input. Our results imply that in several interesting cases, we can implement such algorithms even when the drawing is unknown. [ABSTRACT FROM AUTHOR]

Let F be a compact set of a Banach space X. This paper analyzes the "generalized empirical interpolation method," which, given a function f ∊ F, builds an interpolant Jn[f] in an n dimensional subspace Xn ⊂ X with the knowledge of n outputs (σi(f))ni=1, where σi ∊ X' and X' is the dual space of X. The space Xn is built with a greedy algorithm that is adapted to F in the sense that it is generated by elements of F itself. The algorithm also selects the linear functionals (σi)ni=1 from a dictionary Σ ⊂ X'. In this paper, we study the interpolation error maxf∊F ‖f - Jn[f]‖X by comparing it with the best possible performance on an n dimensional space, i.e., the Kolmogorov n-width of F in X, dn(F,X). For polynomial or exponential decay rates of dn(F,X), we prove that the interpolation error has the same behavior modulo the norm of the interpolation operator. Sharper results are obtained in the case where X is a Hilbert space. [ABSTRACT FROM AUTHOR]

An earlier paper proved the convergence of a variable stepsize Bregman operator splitting algorithm (BOSVS) for minimizing ø(Bu) + H(u), where H and ø are convex functions, and ø is possibly nonsmooth. The algorithm was shown to be relatively efficient when applied to partially parallel magnetic resonance image reconstruction problems. In this paper, the convergence rate of BOSVS is analyzed. When H(u) = ‖Au - f ‖2, where A is a matrix, it is shown that for an ergodic approximation uk obtained by averaging k BOSVS iterates, the error in the objective value ø(Buk)+H(uk) is O(1/k). When the optimization problem has a unique solution u*, we obtain the estimate ‖uk-u*‖ = O(1/ √ k). The theoretical analysis is compared to observed convergence rates for partially parallel magnetic resonance image reconstruction problems where A is a large dense ill-conditioned matrix. [ABSTRACT FROM AUTHOR]

Threshold dynamics and its extensions have proven useful in computing interfacial motions in applications as diverse as materials science and machine learning. Certain desirable properties of the algorithm, such as unconditional monotonicity in two-phase ows and gradient stability more generally, hinge on positivity properties of the convolution kernel and its Fourier transform. Recent developments in the analysis of this class of algorithms indicate that sometimes, as in the case of certain anisotropic curvature ows arising in materials science, these properties of the convolution kernel cannot be expected. Other applications, such as machine learning, would benefit from as great a level of exibility in choosing the convolution kernel as possible. In this paper, we establish certain desirable properties of threshold dynamics, such as gamma convergence of its associated energy, for a substantially wider class of kernels than has been hitherto possible. We also present variants of the algorithm that extend some of these properties to even wider classes of convolution kernels. [ABSTRACT FROM AUTHOR]

The roundoff-error-free (REF) LU and Cholesky factorizations, combined with the REF substitution algorithms, allow rational systems of linear equations to be solved exactly and efficiently by working entirely in integer arithmetic. These REF computational tools share two key properties: their constituent divisions are exact, and their matrix entries' bit-lengths (i.e., required number of storage bits) are bounded polynomially. However, as this paper explains, updating the REF factorizations via the common delete-insert-reduce approach turns out to be inefficient in terms of operand bit-length growth and increased computational effort. In fact, we prove that this inefficiency can eliminate the computational savings expected of the factorization update process. Hence, the current work develops REF update algorithms that differ significantly from their traditional counterparts. To achieve this, it introduces the frame matrix methodology for representing and applying the REF factorizations, and then it devises new REF operations and shortcuts for transitioning between "nearby" frame matrices. The presented push-and-swap update approach utilizes these transitions to preserve the special structure of the REF factorizations and to avoid additional growth in their matrix entries' bit-lengths; this update process requires the expected O(n²) operations. The featured updates are column addition, deletion, and replacement with respect to the REF LU factorization. We also prove that analogous row updates can be performed via the column updates and discuss special considerations for updating the REF Cholesky factorization. An accompanying set of computational tests attests that for fully dense random instances, the push-and-swap column replacement update for the REF LU factorization is on average an order of magnitude faster than the corresponding delete-insert-reduce update for the exact rational arithmetic LU factorization. Altogether, the REF update algorithms represent a promising preliminary step toward enhancing the basic solution validation subroutines employed by state-of-the-art mixed-precision methodologies for exact linear and mixed-integer programming. [ABSTRACT FROM AUTHOR]

We seek discrete approximations to solutions u : Ω → ℝ of semilinear elliptic PDE of the form Δu + fs (u) = 0, where fs is a one-parameter family of nonlinear functions and Ω is a domain in ℝd. The main achievement of this paper is the approximation of solutions to the PDE on the cube Ω = (0,π)³ ⊆ ℝ³. There are 323 possible isotropy subgroups of functions on the cube, which fall into 99 conjugacy classes. The bifurcations with symmetry in this problem are quite interesting, including many with 3-dimensional critical eigenspaces. Our automated symmetry analysis is necessary with so many isotropy subgroups and bifurcations among them, and it allows our code to follow one branch in each equivalence class that is created at a bifurcation point. Our most complicated result is the complete analysis of a degenerate bifurcation with a 6-dimensional critical eigenspace. This article extends our work in [Internat. J. Bifur. Chaos Appl. Sci. Engrg., 19 (2009), pp. 2531-2556], wherein we combined symmetry analysis with modified implementations of the gradient Newton-Galerkin algorithm (GNGA [J. M. Neuberger and J. W. Swift, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 11 (2001), pp. 801-820]) to automatically generate bifurcation diagrams and solution graphics for small, discrete problems with large symmetry groups. The code described in the current paper is efficiently implemented in parallel, allowing us to investigate a relatively fine-mesh discretization of the cube. We use the methodology and corresponding library presented in our previous paper [Int. J. Parallel Program., 40 (2012), pp. 443-464]. [ABSTRACT FROM AUTHOR]

Solving feasibility problems is a central task in mathematics and the applied sciences. One particularly successful method is the Douglas--Rachford algorithm. In this paper, we provide many new conditions sufficient for finite convergence. Numerou s examples illustrate our results. [ABSTRACT FROM AUTHOR]

The development of strong linearizations preserving whatever structure a matrix polynomial might possess has been a very active area of research in the last years, since such linearizations are the starting point of numerical algorithms for computing eigenvalues of structured matrix polynomials with the properties imposed by the considered structure. In this context, Hermitian matrix polynomials are one of the most important classes of matrix polynomials arising in applications and their real eigenvalues are of great interest. The sign characteristic is a set of signs attached to these real eigenvalues which is crucial for determining the behavior of systems described by Hermitian matrix polynomials and, therefore, it is desirable to develop linearizations that preserve the sign characteristic of these polynomials, but, at present, only one such linearization is known. In this paper, we present a complete characterization of all the Hermitian strong linearizations that preserve the sign characteristic of a given Hermitian matrix polynomial and identify several families of such linearizations that can be constructed very easily from the coefficients of the polynomial. [ABSTRACT FROM AUTHOR]

Fast algorithms for matrix multiplication, namely those that perform asymptotically fewer scalar operations than the classical algorithm, have been considered primarily of theoretical interest. Apart from Strassen's original algorithm, few fast algorithms have been efficiently im- plemented or used in practical applications. However, there exist many practical alternatives to Strassen's algorithm with varying performance and numerical properties. Fast algorithms are known to be numerically stable, but because their error bounds are slightly weaker than the classical al- gorithm, they are not used even in cases where they provide a performance benefit. We argue in this paper that the numerical sacrifice of fast algorithms, particularly for the typical use cases of practical algorithms, is not prohibitive, and we explore ways to improve the accuracy both theoret- ically and empirically. The numerical accuracy of fast matrix multiplication depends on properties of the algorithm and of the input matrices, and we consider both contributions independently. We generalize and tighten previous error analyses of fast algorithms and compare their properties. We discuss algorithmic techniques for improving the error guarantees from two perspectives: manipulat- ing the algorithms, and reducing input anomalies by various forms of diagonal scaling. Finally, we benchmark performance and demonstrate our improved numerical accuracy. [ABSTRACT FROM AUTHOR]

The purpose of this paper is to present a theoretic and numerical study of utilizing squeezing and phase shift in coherent feedback control of linear quantum optical systems. A quadrature representation with built-in phase shifters is proposed for such systems. Fundamental structural characterizations of linear quantum optical systems are derived in terms of the new quadrature representation. These results reveal considerable insights into the issue of the physical realizability of such quantum systems. The problem of coherent quantum linear quadratic Gaussian (LQG) feedback control studied in H. I. Nurdin, M. R. James, and I. R. Petersen, Automatica, IFAC, 45 (2009), pp. 1837-1846; G. Zhang and M. R. James, IEEE Trans. Automat. Control, 56 (2011), pp. 1535-1550 is reinvestigated in depth. First, the optimization methods in these papers are extended to a multistep optimization algorithm which utilizes ideal squeezers. Second, a two-stage optimization approach is proposed on the basis of controller parametrization. Numerical studies show that closed-loop systems designed via the second approach may offer LQG control performance even better than that when the closed-loop systems are in the vacuum state. When ideal squeezers in a closed-loop system are replaced by (more realistic) degenerate parametric amplifiers, a sufficient condition is derived for the asymptotic stability of the resultant new closed-loop system; the issue of performance convergence is also discussed in the LQG control setting. [ABSTRACT FROM AUTHOR]

Randomized sampling has recently been proven a highly efficient technique for computing approximate factorizations of matrices that have low numerical rank. This paper describes an extension of such techniques to a wider class of matrices that are not themselves rank-deficient but have off-diagonal blocks that are; specifically, the class of so-called hierarchically semiseparable (HSS) matrices. HSS matrices arise frequently in numerical analysis and signal processing, particularly in the construction of fast methods for solving differential and integral equations numerically. The HSS structure admits algebraic operations (matrix-vector multiplications, matrix factorizations, matrix inversion, etc.) to be performed very rapidly, but only once the HSS representation of the matrix has been constructed. How to rapidly compute this representation in the first place is much less well understood. The present paper demonstrates that if an N x N matrix can be applied to a vector in O(N) time, and if individual entries of the matrix can be computed rapidly, then provided that an HSS representation of the matrix exists, it can be constructed in O(N k²) operations, where k is an upper bound for the numerical rank of the off-diagonal blocks. The point is that when legacy codes (based on, e.g., the fast multipole method) can be used for the fast matrix-vector multiply, the proposed algorithm can be used to obtain the HSS representation of the matrix, and then well-established techniques for HSS matrices can be used to invert or factor the matrix. [ABSTRACT FROM AUTHOR]

In this paper we derive a framework for multiscale approximation of elliptic problems on standard and mixed form. The method presented is based on a splitting into coarse and fine scales together with a systematic technique for approximation of the fine scale part, based on the solution of decoupled localized subgrid problems. The fine scale approximation is then used to modify the coarse scale equations. A key feature of the method is that symmetry in the bilinear form is preserved in the discrete system. Other key features are a posteriori error bounds and adaptive algorithms based on these bounds. The adaptive algorithms are used for automatic tuning of the method parameters. In the last part of the paper we present numerical examples where we apply the framework to a problem in oil reservoir simulation. [ABSTRACT FROM AUTHOR]

This paper introduces a new sweeping preconditioner for the iterative solution of the variable coefficient Helmholtz equation in two and three dimensions. The algorithms follow the general structure of constructing an approximate LDLI factorization by eliminating the unknowns layer by layer starting from an absorbing layer or boundary condition. The central idea of this paper is to approximate the Schur complement matrices of the factorization using moving perfectly matched layers (PMLs) introduced in the interior of the domain. Applying each Schur complement matrix is equivalent to solving a quasi-1D problem with a banded LU factorization in the 2D case and to solving a quasi-2D problem with B multifrontal method in the 3D case. The resulting preconditioner has linear application cost, and the preconditioned iterative solver converges in a number of iterations that is essentially independent of the number of unknowns or the frequency. Numerical results are presented in both two and three dimensions to demonstrate the efficiency of this new preconditioner. [ABSTRACT FROM AUTHOR]

Let A ∈ ℂN x N be a nonsingular complex matrix and b and c be complex vectors of length N. The goal of this paper is to investigate approaches for efficient approximations of the bilinear form c*A-1b. Equivalently, we wish to approximate the scalar value c*x, where x solves the linear system Ax = b. Here the matrix A can be very large or its elements can be too costly to compute so that A is not explicitly available and it is used only in the form of the matrix-vector product. Therefore a direct method is not an option. For A Hermitian positive definite, b*A-1b can be efficiently approximated as a by-product of the conjugate-gradient iterations, which is mathematically equivalent to the matching moment approximations computed via the Gauss-Christoffel quadrature. In this paper we propose a new method using the biconjugate gradient iterations which is applicable to the general complex case. The proposed approach will be compared with existing ones using analytic arguments and numerical experiments. [ABSTRACT FROM AUTHOR]

This paper concerns numerical methods for computing complex geometrical optics (CGO) solutions to the conductivity equation ∇ ⋅ σ∇u (⋅, k) = 0 in ℝ² for piecewise smooth conductivities σ, where k is a complex parameter. The key is to solve an R-linear singular integral equation defined in the unit disk. Recently, Astala et al. [Appl. Comput. Harmon. Anat., 29 (2010), pp. 2-17] proposed a complicated method for numerical computation of CGO solutions by solving a periodic version of the ℝ-linear integral equation in a rectangle containing the unit disk. In this paper, based on the fast algorithms in [P. Daripa and D. Mashat, Numer. Algorithms, 18 (1998), pp. 133-157] for singular integral transforms, we propose a simpler numerical method which solves the ℝ-linear integral equation in the unit disk directly. For the resulting ℝ-linear operator equation, a minimal residual iterative method is proposed. Numerical examples illustrate the accuracy and efficiency of the new method. [ABSTRACT FROM AUTHOR]

As a solution to the frequently reported problem of the high frequencies reconstructed sometimes many iterations after the low frequencies in tomography, we proposed in a recent paper [V. Israel-Jost, Ph. Choquet, and A. Constantinesco, Int. J. Biomed. Imaging, 2006, paper 34043.] a general adaptation that can apply to many algorithms: the incomplete backprojection. We stressed the fact that this type of frequency adapted (FA) algorithm works only when the linear problem of tomography to be solved involves deblurring, i.e., when the point spread functions (PSFs) associated with each voxel are large. Thus, depending on the desired level of detail, a parameter is set to achieve the reconstruction within a small number of iterations. The aim of this paper is to give a mathematical analysis of both the acceleration obtained and the convergence of an FA algorithm derived from the simultaneous algebraic reconstruction technique (SART) algorithm of Andersen and Kak, the so-called FA-SART algorithm. [ABSTRACT FROM AUTHOR]

This paper is concerned with the development of well-balanced high order numerical schemes for systems of balance laws with a linear flux function, whose coefficients may be variable. First, well-balanced first order numerical schemes are obtained based on the use of exact solvers of Riemann problems that include both the flux and the source terms. Godunov's methods so obtained are extended to higher order schemes by using a technique of reconstruction of states. The main contribution of this paper is to introduce a reconstruction technique that preserves the well-balanced property of Godunov's methods. Some numerical experiments are presented to verify in practice the properties of the developed numerical schemes. [ABSTRACT FROM AUTHOR]

In this paper we consider the numerical solution of discrete boundary integral equations on polyhedral surfaces in three dimensions. When the solution contains typical edge singularities, highly stretched meshes are preferred to uniform meshes, since they reduce the number of degrees of freedom needed to obtain a fixed accuracy. The classical panel-clustering method can still be applied in the presence of such highly stretched meshes. However, we will show that the savings in computation time and storage become suboptimal because the nearfield matrix arising in the panel-clustering algorithm is no longer as sparse as it is in the case of uniform meshes. Hence, a natural question arises as to whether a new enhanced panel-clustering algorithm can be designed which performs efficiently even in the presence of highly stretched meshes. The main result of this paper is to formulate such an enhanced version of the panel-clustering algorithm. The key features of the algorithm are (i) the employment of partial analytic integration in the direction of stretching, yielding a new kernel function on a one-dimensional manifold where the influence of high aspect ratios in the stretched elements is removed, and (ii) the introduction of a generalized admissibility condition with respect to the partially integrated kernel, which ensures that certain stretched clusters which are inadmissible in the classical sense now become admissible. In the context of a model problem, we prove that our algorithm yields an accurate (up to the discretization error) matrix-vector multiplication which requires ...(NlogκN) operations, where N is the number of degrees of freedom and k is small and independent of the aspect ratio. The generalized admissibility condition can be viewed as an addition to the classical method which may be useful in general when stretched meshes are present. We also have performed a numerical experiment which shows that the sparsity of the nearfield matrix for the enhanced panel-clustering method is not negatively affected by stretched elements, and the method will perform optimally. [ABSTRACT FROM AUTHOR]

In this paper, we show how to perform the Hessenberg reduction of a rank structured matrix under unitary similarity operations in an efficient way, using the Givens-weight representation. This reduction can be used as a first step for eigenvalue computation. We also show how the algorithm can be modified to compute the bidiagonal reduction of a rank structured matrix, this latter method being a preprocessing step for computing the singular values of the matrix. For the main cases of interest, the algorithms we describe in this paper are of complexity O((ar + bs)n2), where n is the matrix size, r is some measure for the average rank index of the rank structure, s is some measure for the bandwidth of the unstructured matrix part around the main diagonal, and a, b ? R are certain weighting parameters. Numerical experiments demonstrate the stability of this approach. [ABSTRACT FROM AUTHOR]

In this paper, an efficient and stable numerical algorithm for computing multiple saddle points with symmetries is developed by modifying the local minimax method established in [Y. Li and J. Zhou, SIAM J. Sci. Comput. 23 (2001), pp. 840--865; Y. Li and J. Zhou, SIAM J. Sci. Comput., 24 (2002), pp. 840--865]. First an invariant space is defined in a more general sense and a principle of invariant criticality is proved for the generalization. Then the orthogonal projection to the invariant space is used to preserve the invariance and to reduce computational error across iterations. Simple averaging formulas are used for the orthogonal projections. Numerical computations of examples with various symmetries, of which some can and others cannot be characterized by a compact group of linear isomorphisms, are carried out to confirm the theory and to illustrate applications. The mathematical features of various problems demonstrated in these examples fall into two categories: nodal solutions of saddle-point type with large Morse indices and nonradial positive solutions via symmetry breaking in radially symmetric elliptic problems. The new numerical algorithm generates these rather unstable solutions in an efficient and stable way. The existence of many unstable solutions and their behavior found in this paper remain to be investigated. [ABSTRACT FROM AUTHOR]

Local mesh-refining algorithms known from adaptive finite element methods are adopted for locally conservative and monotone finite volume discretizations of boundary value problems for steady-state convection-diffusion-reaction equations. The paper establishes residual-type explicit error estimators and averaging techniques for a posteriori finite volume error control with and without upwind in global H1- and L2-norms. Reliability and efficiency are verified theoretically and confirmed empirically with experimental support for the superiority of the suggested adaptive mesh-refining algorithms over uniform mesh refining. A discussion of adaptive computations in the simulation of contaminant concentration in a nonhomogeneous water reservoir concludes the paper. [ABSTRACT FROM AUTHOR]

In this paper we use a displacement structure approach to design a class of new preconditioners for the conjugate gradient method applied to the solution of large Toeplitz linear equations. Explicit formulas are suggested for the G.Strang-type and for the T.Chan-type preconditioners belonging to any of eight classes of matrices diagonalized by the corresponding discrete cosine or sine transforms. Under the standard Wiener class assumption the clustering property is established for all of these preconditioners, guaranteeing rapid convergence of the preconditioned con- jugate gradient method. All the computations related to the new preconditioners can be done in real arithmetic, and to fully exploit this advantageous property one has to suggest a fast real-arithmetic algorithm for multiplication of a Toeplitz matrix by a vector. It turns out that the obtained formulas for the Strang-type preconditioners allow a number of representations for Toeplitz matrices leading to a wide variety of real-arithmetic multiplication algorithms based on any of eight discrete cosine or sine transforms. Recently, transformations of Toeplitz matrices to Vandermonde-like or Cauchy-like matrices have been found to be useful in developing accurate direct methods for Toeplitz linear equations. In this paper we suggest further extending the range of the transformation approach by exploring it for iterative methods; this technique allowed us to reduce the complexity of each iteration of the preconditioned conjugate gradient method. The results of this paper were announced in [T. Kailath and V. Olshevsky, Calcolo, 33 (1996), pp. 191-208]. [ABSTRACT FROM AUTHOR]

In (SIAM J. Matrix Anal. Appl., 23 (2002), pp. 706–727, we developed numerical algorithms for computing sparse low-rank approximations of matrices, and we also provided a detailed error analysis of the proposed algorithms together with some numerical experiments. The low-rank approximations are constructed in a certain factored form with the degree of sparsity of the factors controlled by some user-specified parameters. In this paper, we cast the sparse low-rank approximation problem in the framework of penalized optimization problems. We discuss various approximation schemes for the penalized optimization problem which are more amenable to numerical computations. We also include some analysis to show the relations between the original optimization problem and the reduced one. We then develop a globally convergent discrete Newton-like iterative method for solving the approximate penalized optimization problems. We also compare the reconstruction errors of the sparse low-rank approximations computed by our new methods with those obtained using the methods in the earlier paper and several other existing methods for computing sparse low-rank approximations. Numerical examples show that the penalized methods are more robust and produce approximations with factors which have fewer columns and are sparser. [ABSTRACT FROM AUTHOR]

This paper presents numerical results for the asynchronous version of the fast adaptive composite-grid algorithm (AFACx). These results confirm the level-independent convergence bounds established theoretically in a companion paper. These numerical results include the case of AFACx applied to first-order system least-squares finite element discretizations of the stationary Stokes equations on curvilinear adaptive mesh refinement grids. [ABSTRACT FROM AUTHOR]

This paper discusses convergence behavior of the nonlinear inexact Uzawa algorithm for solving saddle-point problems presented in a recent paper of Bramble, Pasciak, and Vassilev (SIAM J. Numer. Anal., 34 (1997), pp. 1072­1092). It is shown that this algorithm converges under a condition weaker than that stated in their paper. Key words: saddle-point problems; nonlinear inexact Uzawa algorithm; convergence behavior [ABSTRACT FROM AUTHOR]

In an earlier paper, we developed a new formulation of the shallow-water equations on the sphere, based on the semi-Lagrangian advection of potential vorticity. An important feature of our two-time-level model is that forward extrapolation of the nonlinear terms is avoided. In this paper, we develop a multigrid algorithm for solving the resulting system of three coupled nonlinear equations. We include a discussion of basic multigrid methods, along with some more advanced techniques required for the development of the multigrid scheme for this system. Results are presented showing that usual optimal multigrid efficiency is obtained. [ABSTRACT FROM AUTHOR]

Dynamical systems can often be decomposed into loosely coupled subsystems. The system of ordinary differential equations (ODEs) modelling such a problem can then be partitioned corresponding to the subsystems, and the loose couplings can be exploited by special integration methods to solve the problem using a parallel computer or just solve the problem more efficiently than by standard methods. This paper presents accuracy analysis of methods for the numerical integration of stiff partitioned systems of ODEs. The discretization formulas are based on the implicit Euler formula and the second order implicit backward differentiation formula (BDF2). Each subsystem of the partitioned problem is discretized independently, and the couplings to the other subsystems are based on solution values from previous time steps. Applied this way, the discretization formulas are called decoupled. The stability properties of the decoupled implicit Euler formula are well understood. This paper presents error bounds and asymptotic error expansions to be used in controlling step size, relaxation between subsystems and the validity of the partitioning. The decoupled BDF2 formula is analyzed within the same framework. Finally, the analysis is used in the design of a decoupled numerical integration algorithm with variable step size to control the local error and adaptive selection of partitionings. Two versions of the algorithm with decoupled implicit Euler and BDF2, respectively, are used in examples where a realistic problem is solved. The examples compare the results from the decoupled implicit Euler and BDF2 formulas, and compare them with results from the corresponding classical formulas. [ABSTRACT FROM AUTHOR]

This paper presents a new algorithm for accurate floating-point computation of the singular value decomposition (SVD) of the product A = BτSC, where B ∈ Rp×m, C ∈ Rq×n, S ∈ Rp×q, and p ≤ m, q ≤ n. The new algorithm uses diagonal scalings, the LU factorization with complete pivoting, the QR factorization with column pivoting, and matrix multiplication to replace A by A′ = B′τS′C′, where A and A′ have the same singular values and the matrix A′ is computed explicitly. The singular values of A′ are computed using the Jacobi SVD algorithm. It is shown that the accuracy of the new algorithm is determined by (i) the accuracy of the QR factorizations of Bτ and Cτ; (ii) the accuracy of the LU factorization with complete pivoting of S; and (iii) the accuracy of the computation of the SVD of a matrix A′ with moderate minD=diagκ2(A′D). Theoretical analysis and numerical evidence show that, in the case of rank(B) = rank(C) = p and full rank S, the accuracy of the new algorithm is unaffected by replacing B, S, C with, respectively, D1B, D2SD3, D4C, where Di, i = 1,…,4, are arbitrary diagonal matrices. As an application, the paper proposes new accurate algorithms for computing the (H,K)­SVD and (H-1,K)­SVD of S. [ABSTRACT FROM AUTHOR]

Max-plus methods have previously been used to solve deterministic control problems. The methods are based on max-plus (or min-plus) expansions and can yield curse-of-dimensionalityfree numerical methods. In this paper, we explore min-plus methods for continuous-time stochastic control on a finite-time horizon. We first approximate the original value function via timediscretization. By generalizing the min-plus distributive property to continuum spaces, we obtain an algorithm for recursive computation of the time-discretized values, which we refer to as the idempotent distributed dynamic programming principle (IDDPP). Under the IDDPP, the value function at each step can be represented as an infimum of functions in a certain class. This is a min-plus expansion for the value function. For the specific class of problems considered here, we see that the class can be taken as that consisting of the quadratic functions. A means for reducing the numbers of constituent quadratic functions is discussed. [ABSTRACT FROM AUTHOR]

In this paper, we introduce the uncertain optimal control problem of determining a control that minimizes the expectation of an objective functional for a system with parameter uncertainty in both dynamics and objective. We present a computational framework for the numerical solution of this problem, wherein an independently drawn random sample is taken from the space of uncertain parameters, and the expectation in the objective functional is approximated by a sample average. The result is a sequence of approximating standard optimal control problems that can be solved using existing techniques. To analyze the performance of this computational framework, we develop necessary conditions for both the original and approximate problems and show that the approximation based on sample averages is consistent in the sense of Polak [Optimization: Algorithms and Consistent Approximations, Springer, New York, 1997]. This property guarantees that accumulation points of a sequence of global minimizers (stationary points) of the approximate problem are global minimizers (stationary points) of the original problem. We show that the uncertain optimal control problem can further be approximated in a consistent manner by a sequence of nonlinear programs under mild regularity assumptions. In numerical examples, we demonstrate that the framework enables the solution of optimal search and optimal ensemble control problems. [ABSTRACT FROM AUTHOR]

This paper considers an implicitly restarted Krylov subspace method that approximates a stable, linear transfer function f(s) of order n by one of order m, where n >> m. It is well known that oblique projections onto a Krylov subspace may generate unstable partial realizations. To remedy this situation, the oblique projectors obtained via classical Krylov subspace methods are supplemented with further projectors which enable the formation of stable partial realizations directly from f(s). A key feature of this process is that it may be incorporated into an implicit restart scheme. A second difficulty arises from the fact that Krylov subspace methods often generate partial realizations that contain nonessential modes. To this end, balanced truncation may be employed to discard the unwanted part of the reduced-order model. This paper proposes oblique projection methods for large-scale model reduction that simultaneously compute stable reduced-order models while discarding all nonessential modes. It is shown that both of these tasks may be effected by a single oblique projection process. Furthermore, the process is shown to naturally fit into an implicit restart framework. The theoretical properties of these methods are thoroughly investigated, and exact low-dimensional expressions for the L∞-norm of the residual errors are derived. Finally, the behavior of the algorithm is illustrated on two large-scale examples. [ABSTRACT FROM AUTHOR]

The generalized null space decomposition (GNSD) is a unitary reduction of a general matrix A of order n to a block upper triangular form that reveals the structure of the Jordan blocks of A corresponding to a zero eigenvalue. The reduction was introduced by Kublanovskaya. It was extended first by Ruhe and then by Golub and Wilkinson, who based the reduction on the singular value decomposition. Unfortunately, if A has large Jordan blocks, the complexity of these algorithms can approach the order of n4. This paper presents an alternative algorithm, based on repeated updates of a QR decomposition of A, that is guaranteed to be of order n³. Numerical experiments confirm the stability of this algorithm, which turns out to produce essentially the same form as that of Golub and Wilkinson. The effect of errors in A on the ability to recover the original structure is investigated empirically. Several applications are discussed, including the computation of the Drazin inverse. [ABSTRACT FROM AUTHOR]

In this paper we discuss the deflation criterion used in the extended QR algorithm based on the chasing of rotations. We provide absolute and relative perturbation bounds for this deflation criterion. Further, we present a generalization of aggressive early deflation to the extended QR algorithms. Aggressive early deflation is the key technique for the identification and deflation of already converged, but hidden, eigenvalues. Often these possibilities for deflation are not detected by the standard technique. We present numerical results underpinning the power of aggressive early deflation also in the context of extended QR algorithms. We further generalize these ideas by the transcription of middle deflations. [ABSTRACT FROM AUTHOR]

This paper reports the multiscale analysis and numerical algorithms for a class of distributed elliptic optimal control problems with constraints. Multiscale asymptotic expansions are presented, and an explicit rate of convergence is derived. A multiscale numerical method is developed. Numerical experiments are carried out to validate the predicted convergence results. [ABSTRACT FROM AUTHOR]

In this paper we study the problem of finding a weighted geometric mean sufficiently close to the weighted Karcher mean of positive definite matrices. This problem arises in various numerical algorithms for computing the weighted Karcher mean and in a physical problem of averaging on the Riemannian manifold of positive definite matrices. We prove that this is possible for any "convex" geometric means satisfying the Jensen-type inequality for geodesically convex functions varying over an open set of weights in the simplex of positive probability vectors. This follows from a useful estimation of the Riemannian distance between the Karcher mean and any convex geometric means. A better upper bound between the Karcher mean and the Riemannian convex combination admitting an explicit formula is presented. [ABSTRACT FROM AUTHOR]

Stochastic models of chemical systems are often analyzed by solving the corresponding Fokker--Planck equation, which is a drift-diffusion partial differential equation for the probability distribution function. Efficient numerical solution of the Fokker--Planck equation requires adaptive mesh refinements. In this paper, we present a mesh refinement approach which makes use of a stochastic simulation of the underlying chemical system. By observing the stochastic trajectory for a relatively short amount of time, the areas of the state space with nonnegligible probability density are identified. By refining the finite element mesh in these areas, and coarsening elsewhere, a suitable mesh is constructed and used for the computation of the stationary probability density. Numerical examples demonstrate that the presented method is competitive with existing a posteriori methods. [ABSTRACT FROM AUTHOR]

In this paper a new O(N log³ N) solver for N x 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 2x2 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. [ABSTRACT FROM AUTHOR]

Mesh generation by Delaunay refinement is a widely used technique for constructing guaranteed quality triangular and tetrahedral meshes. The quality guarantees are usually provided in terms of the bounds on circumradius-to-shortest-edge ratio and on the grading of the resulting mesh. Traditionally circumcenters of skinny elements and middle points of boundary faces and edges are used for the positions of inserted points. However, recently variations of the traditional algorithms are being proposed that are designed to achieve certain optimization objectives by inserting new points in neighborhoods of the center points. In this paper we propose a general approach to the selection of point positions by defining one-, two-, and three-dimensional selection regions such that any point insertion strategy based on these regions is automatically endowed with the theoretical guarantees proven here. In particular, for the input models defined by planar linear complexes under the assumption that no input angle is less than 90°, we prove the termination of the proposed generalized algorithm, as well as the fidelity and good grading of the resulting meshes. [ABSTRACT FROM AUTHOR]

Recently a generalization of Francis's implicitly shifted QR algorithm was proposed, notably widening the class of matrices admitting low-cost implicit QR steps. This unifying framework covered the methods and theory for Hessenberg and inverse Hessenberg matrices and furnished also new, single-shifted, QR-type methods for, e.g., CMV matrices. Convergence of this approach was only suggested by numerical experiments. No theoretical claims supporting the results were presented. In this paper we present multishift variants of these new algorithms. We also provide a convergence theory that shows that the new algorithm performs nested subspace iterations on rational Krylov subspaces. Numerical experiments confirm the validity of the theory. [ABSTRACT FROM AUTHOR]

Modeling atomic and molecular systems requires computation-intensive quantum mechanical methods such as, but not limited to, density functional theory [R. A. Friesner, Proc. Natl. Acad. Sci. USA, 102 (2005), pp. 6648-6653]. These methods have been successful in predicting various properties of chemical systems at atomistic scales. Due to the inherent nonlocality of quantum mechanics, the scalability of these methods ranges from O(N³)to O (N7) depending on the method used and approximations involved. This significantly limits the size of simulated systems to a few thousand atoms, even on large scale parallel platforms. On the other hand, classical approximations of quantum systems, although computationally (relatively) easy to implement, yield simpler models that lack essential chemical properties such as reactivity and charge transfer. The recent work of van Duin et al. [J. Phys. Chem. A, 105 (2001), pp. 9396-9409] overcomes the limitations of nonreactive classical molecular dynamics (MD) approximations by carefully incorporating limited nonlocality (to mimic quantum behavior) through an empirical bond order potential. This reactive classical MD method, called ReaxFF, achieves essential quantum properties, while retaining the computational simplicity of classical MD, to a large extent. Implementation of reactive force fields presents significant algorithmic challenges. Since these methods model bond breaking and formation, efficient implementations must rely on complex dynamic data structures. Charge transfer in these methods is accomplished by minimizing electrostatic energy through charge equilibration. This requires the solution of large linear systems (108 degrees of freedom and beyond) with shielded electrostatic kernels at each time-step. Individual time-steps are themselves typically in the range of tenths of femtoseconds, requiring optimizations within and across time-steps to scale simulations to nanoseconds and beyond, where interesting phenomena may be observed. In this paper, we present implementation details of sPuReMD (serial Purdue reactive molecular dynamics program), a unique reactive classical MD code. We describe various data structures, and the charge equilibration solver at the core of the simulation engine. This Krylov subspace solver relies on a preconditioner based on incomplete LU factorization with thresholds (ILUT), specially targeted to our application. We comprehensively validate the performance and accuracy of sPuReMD on a variety of hydrocarbon systems. In particular, we show excellent per-time-step time, linear time scaling in system size, and a low memory footprint. sPuReMD is a freely distributed software with GPL and is currently being used to model diverse systems ranging from oxidative stress in biomembranes to strain relaxation in Si-Ge nanorods. [ABSTRACT FROM AUTHOR]

The dqds algorithm computes all the singular values of an n x n bidiagonal matrix to high relative accuracy in O(n²) cost. Its efficient implementation is now available as a LAPACK subroutine and is the preferred algorithm for this purpose. In this paper we incorporate into dqds a technique called aggressive early deflation, which has been applied successfully to the Hessenberg QR algorithm. Extensive numerical experiments show that aggressive early deflation often reduces the dqds runtime significantly. In addition, our theoretical analysis suggests that with aggressive early deflation, the performance of dqds is largely independent of the shift strategy. We confirm through experiments that the zero-shift version is often as fast as the shifted version. We give a detailed error analysis to prove that with our proposed deflation strategy, dqds computes all the singular values to high relative accuracy. [ABSTRACT FROM AUTHOR]

This paper is on the construction of energy-minimizing coarse spaces that obey certain functional constraints and can thus be used, for example, to build robust coarse spaces for elliptic problems with large variations in the coefficients. In practice they are built by patching together solutions to appropriate local saddle point or eigenvalue problems. We develop an abstract framework for such constructions, akin to an abstract Bramble--Hilbert-type lemma, and then apply it in the design of coarse spaces for discretizations of PDEs with highly varying coefficients. The stability and approximation bounds of the constructed interpolant are in the weighted L2 norm and are independent of the variations in the coefficients. Such spaces can be used, for example, in two-level overlapping Schwarz algorithms for elliptic PDEs with large coefficient jumps generally not resolved by a standard coarse grid or for numerical upscaling purposes. Some numerical illustration is provided. [ABSTRACT FROM AUTHOR]

We present two iterative algorithms for solving real nonsymmetric and complex nonHermitian linear systems of equations and that were developed from variants of the nonsymmetric Lanczos method. In this paper, we give the theoretical background of the two iterative methods and discuss their main computational aspects. Using a large number of numerical experiments, we analyze their convergence properties, and we also compare them with other popular nonsymmetric iterative solvers in use today. [ABSTRACT FROM AUTHOR]

Anisotropic diffusion methods are well known for giving good qualitative results for image denoising. This paper gives a review of the anisotropic diffusion methodology and its application to image denoising. We propose a fixed-point iteration using a multigrid solver to solve a regularized anisotropic diffusion equation, which is not only well-posed, but also has a nontrivial steady-state solution. A new regularization parameter-choice method (Brent-NCP), combining Brent's method and the normalized cumulative periodogram information of the misfit, is also introduced. We test our algorithm on several common test images with different noise levels. The experimental results demonstrate the effectiveness of the anisotropic diffusion with a multigrid approach and the broad applicability of the Brent-NCP parameter-choice algorithm. [ABSTRACT FROM AUTHOR]

This paper investigates backscatter data for the inverse obstacle problem in impedance tomography, when the obstacles are small. It is shown that under this circumstance the backscatter data are close approximations of a rational function that has second or fourth order poles at the locations of the obstacles. Furthermore, a numerical method is presented to locate the obstacles via the poles of certain Laurent-Padé approximations. Numerical experiments explore the potential of this algorithm also for extended obstacles, taking into consideration that the problem is severely ill-posed. [ABSTRACT FROM AUTHOR]

Linear advection of a scalar quantity by a specified velocity field arises in a number of different applications. Of particular interest here is the transport of species and energy in low Mach number models for combustion, atmospheric flows, and astrophysics, as well as contaminant transport in Darcy models of saturated subsurface flow. An important characteristic of these problems is that the velocity field is not known analytically. Instead, an auxiliary equation is solved to compute averages of the velocities over faces in a finite volume discretization. In this paper, we present a customized three-dimensional finite volume advection scheme for this class of problems that provides accurate resolution for smooth problems while avoiding undershoot and overshoot for nonsmooth profiles. The method is an extension of an algorithm by Bell, Dawson, and Shubin (BDS), which was developed for a class of scalar conservation laws arising in porous media flows in two dimensions. The original BDS algorithm is a variant of unsplit. higher-order Godunov methods based on construction of a limited bilinear profile within each computational cell. Here we present a three-dimensional extension of the original BDS algorithm that is based on a limited trilinear profile within each cell. We compare this new method to several other unsplit approaches, including piecewise linear methods, piecewise parabolic methods, and wave propagation schemes. [ABSTRACT FROM AUTHOR]

This paper presents an efficient multiscale stochastic framework for uncertainty quantification in modeling of flow through porous media with multiple rock types. The governing equations are based on Darcy's law with nonstationary stochastic permeability represented as a sum of local Karhunen-Loève expansions. The approximation uses stochastic collocation on either a tensor product or a sparse grid, coupled with a domain decomposition algorithm known as the multiscale mortar mixed finite element method. The latter method requires solving a coarse scale mortar interface problem via an iterative procedure. The traditional implementation requires the solution of local fine scale linear systems on each iteration. We employ a recently developed modification of this method that precomputes a multiscale flux basis to avoid the need for subdomain solves on each iteration. In the stochastic setting, the basis is further reused over multiple realizations, leading to collocation algorithms that are more efficient than the traditional implementation by orders of magnitude. Error analysis and numerical experiments are presented. [ABSTRACT FROM AUTHOR]