Preconditioned Krylov subspace method for the solution of least-squares problems

We consider preconditioned Krylov subspace iteration methods, e.g., CG, LSQR and GMRES, for the solution of large sparse least-squares problems min ∥Ax – b ∥₂, with A ∈ R ^m×n, based on the Krylov subspaces K_k (BA, Br) and K_k (AB, r), respectively, where B ∈ R ^n×m is the preconditioning matrix. More concretely, we propose and implement a class of incomplete QR factorization preconditioners based on the Givens rotations and analyze in detail the efficiency and robustness of the correspondingly preconditioned Krylov subspace iteration methods. A number of numerical experiments are used to further examine their numerical behaviour. It is shown that for both overdetermined and underdetermined least-squares problems, the preconditioned GMRES methods are the best for large, sparse and ill-conditioned matrices in terms of both CPU time and iteration step. Also, comparisons with the diagonal scaling and the RIF preconditioners are given to show the superiority of the newly-proposed GMRES-type methods. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Time decay of solutions of semilinear strongly damped generalized wave equations

We study the asymptotic behaviour in time of the solutions of dissipative perturbations of wave-type equations in ?^N, u_tt + Bu_t + Au + G(u) = 0, with commuting positive operators A, B and a power like non-linearity G(u). First we give some (pseudo) conformal invariants of the linear operator in the equation. This allows us to derive optimal decay rates for the solutions of the linearized problems. We then prove some decay estimates for the non-linear problems using the tools of scattering theory and the aforementioned conformal invariants.     

Computation of Morse decompositions for semilinear elliptic PDEs

A numerical method for computing all solutions of an elliptic boundary value problem Au + g[u, λ] = 0 and their Morse indices as steady‐states of the parabolic problem u_t + Au + g[u, λ] = 0 is presented. Morse decompositions are also determined. The method uses a finite element approach that is based on the method of alternative problems. Error estimates for the finite element approximations are verified and examples are given. © John Wiley & Sons, Inc. Numer Methods Partial Differential Eq 17: 290–312, 2001     

A linear system solver based on a modified Krylov subspace method for breakdown recovery

Despite its usefulness in solving eigenvalue problems and linear systems of equations, the nonsymmetric Lanczos method is known to suffer from a potential breakdown problem. Previous and recent approaches for handling the Lanczos exact and near-breakdowns include, for example, the look-ahead schemes by Parlett-Taylor-Liu [23], Freund-Gutknecht-Nachtigal [9], and Brezinski-Redivo Zaglia-Sadok [4]; the combined look-ahead and restart scheme by Joubert [18]; and the low-rank modified Lanczos scheme by Huckle [17]. In this paper, we present yet another scheme based on a modified Krylov subspace approach for the solution of nonsymmetric linear systems. When a breakdown occurs, our approach seeks a modified dual Krylov subspace, which is the sum of the original subspace and a new Krylov subspaceK _m (w _j ,A ^T ) wherew _j is a newstart vector (this approach has been studied by Ye [26] for eigenvalue computations). Based on this strategy, we have developed a practical algorithm for linear systems called the MLAN/QM algorithm, which also incorporates the residual quasi-minimization as proposed in [12]. We present a few convergence bounds for the method as well as numerical results to show its effectiveness.Research supported by Natural Sciences and Engineering Research Council of Canada.     

Application of a Mickens finite‐difference scheme to the cylindrical Bratu‐Gelfand problem

The boundary value problem Δu + λe^u = 0 where u = 0 on the boundary is often referred to as “the Bratu problem.” The Bratu problem with cylindrical radial operators, also known as the cylindrical Bratu‐Gelfand problem, is considered here. It is a nonlinear eigenvalue problem with two known bifurcated solutions for λ < λ_c, no solutions for λ > λ_c and a unique solution when λ = λ_c. Numerical solutions to the Bratu‐Gelfand problem at the critical value of λ_c = 2 are computed using nonstandard finite‐difference schemes known as Mickens finite differences. Comparison of numerical results obtained by solving the Bratu‐Gelfand problem using a Mickens discretization with results obtained using standard finite differences for λ < 2 are given, which illustrate the superiority of the nonstandard scheme. © 2004 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 20: 327–337, 2004     

Single‐cell discretization of O(kh2 + h4) for ∂u/∂n for three‐space dimensional mildly quasi‐linear parabolic equation

In this article, using a single computational cell, we report some stable two‐level explicit finite difference approximations of O(kh² + h⁴) for ?u/?n for three‐space dimensional quasi‐linear parabolic equation, where h > 0 and k > 0 are mesh sizes in space and time directions, respectively. When grid lines are parallel to x‐, y‐, and z‐coordinate axes, then ?u/?n at an internal grid point becomes ?u/?x, ?u/?y, and ?u/?z, respectively. The proposed methods are also applicable to the polar coordinates problems. The proposed methods have the simplicity in nature and use the same marching type of technique of solution. Stability analysis of a linear difference equation and computational efficiency of the methods are discussed. The results of numerical experiments are compared with exact solutions. © 2003 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 19: 327–342, 2003.     

Spectral Criteria of Abstract Sequences; Application to Difference Equations

We suppose that M is a closed subspace of l ^∞(J, X), the space of all bounded sequences {x(n)} _n?J ? X, where J ? {Z⁺,Z} and X is a complex Banach space. We define the M-spectrum σ_M (u) of a sequence u ? l ^∞(J,X). Certain conditions will be supposed on both M and σ_M (u) to insure the existence of u ? M. We prove that if u is ergodic, such that σ_M (u,) is at most countable and, for every λ ? σ_M (u), the sequence e^?iλnu(n) is ergodic, then u ? M. We apply this result to the operator difference equationu(n + 1) = Au(n) + ψ(n), n ? J,and to the infinite order difference equation Σ ^r _k=1 a_k (u(n + k) ? u(n)) + Σ _{s ? Z}?(n ? s)u(s) = h(n), n?J, where ψ?l ^∞(Z,X) such that ψ| _J ? M, A is the generator of a C ₀-semigroup of linear bounded operators {T(t)} _t>0 on X, h ? M, ? ? l ¹(Z) and a_k ?C. Certain conditions will be imposed to guarantee the existence of solutions in the class M.     

A second-order,chaos-free,explicit method for the numerical solution of a cubic reaction problem in neurophysiology

A second order explicit method is developed for the numerical solution of the initialvalue problem w′(t) ≡ dw(t)/dt = ?(w), t > 0, w(0) = W⁰, in which the function ?(w) = αw(1 ? w) (w ? a), with α and a real parameters, is the reaction term in a mathematical model of the conduction of electrical impulses along a nerve axon. The method is based on four first-order methods that appeared in an earlier paper by Twizell, Wang, and Price [Proc. R. Soc. (London) A 430 , 541–576 (1990)]. In addition to being chaos free and of higher order, the method is seen to converge to one of the correct steady-state solutions at w = 0 or w = 1 for any positive value of α. Convergence is monotonic or oscillatory depending on W⁰, α, a, and l, the parameter in the discretization of the independent variable t. The approach adopted is extended to obtain a numerical method that is second order in both space and time for solving the initial-value boundary-value problem ?u/?t = κ?²u/?x² + αu(1 ? u)(u ? a) in which u = u(x,t). The numerical method so developed obtained the solution by solving a single linear algebraic system at each time step. © 1993 John Wiley & Sons, Inc.     

Method for Solving the Multidimensional n‐Wave Resonant Equations and Geometry of Generalized Darboux‐Manakov‐Zakharov Systems

The intrinsic geometric properties of generalized Darboux‐Manakov‐Zakharov systems of semilinear partial differential equations (1) for a real‐valued function u(x₁, …, x_n) are studied with particular reference to the linear systems in this equation class. System (1) is overdetermined and will not generally be involutive in the sense of Cartan: its coefficients will be constrained by complicated nonlinear integrability conditions. We derive tools for explicitly constructing involutive systems of the form (1) , essentially solving the integrability conditions. Specializing to the linear case provides us with a novel way of viewing and solving the multidimensional n‐wave resonant interaction system and its modified version. For each integer n≥ 3 and nonnegative integer k, our procedure constructs solutions of the n‐wave resonant interaction system depending on at least k arbitrary functions each of one variable. The construction of these solutions relies only on differentiation, linear algebra, and the solution of ordinary differential equations.     

Efficient and reliable iterative methods for linear systems

The approximate solutions in standard iteration methods for linear systems Ax=b, with A an n by n nonsingular matrix, form a subspace. In this subspace, one may try to construct better approximations for the solution x. This is the idea behind Krylov subspace methods. It has led to very powerful and efficient methods such as conjugate gradients, GMRES, and Bi-CGSTAB. We will give an overview of these methods and we will discuss some relevant properties from the user's perspective view.The convergence of Krylov subspace methods depends strongly on the eigenvalue distribution of A, and on the angles between eigenvectors of A. Preconditioning is a popular technique to obtain a better behaved linear system. We will briefly discuss some modern developments in preconditioning, in particular parallel preconditioners will be highlighted: reordering techniques for incomplete decompositions, domain decomposition approaches, and sparsified Schur complements.     

Inverse subspace problems with applications

Given a square matrix A, the inverse subspace problem is concerned with determining a closest matrix to A with a prescribed invariant subspace. When A is Hermitian, the closest matrix may be required to be Hermitian. We measure distance in the Frobenius norm and discuss applications to Krylov subspace methods for the solution of large‐scale linear systems of equations and eigenvalue problems as well as to the construction of blurring matrices. Extensions that allow the matrix A to be rectangular and applications to Lanczos bidiagonalization, as well as to the recently proposed subspace‐restricted SVD method for the solution of linear discrete ill‐posed problems, also are considered.Copyright © 2013 John Wiley & Sons, Ltd.     

Extinction of solutions to a class of fast diffusion systems with nonlinear sources

Y. Wang In this paper, the finite time extinction of solutions to the fast diffusion system u_t=div(|?u|^{p ? 2}?u) + v^m, v_t=div(|?v|^{q ? 2}?v) + uⁿ is investigated, where 1 < p,q < 2, m,n > 0 and is a bounded smooth domain. After establishing the local existence of weak solutions, the authors show that if mn > (p ? 1)(q ? 1), then any solution vanishes in finite time provided that the initial data are ‘comparable’; if mn = (p ? 1)(q ? 1) and Ω is suitably small, then the existence of extinction solutions for small initial data is proved by using the De Giorgi iteration process and comparison method. On the other hand, for 1 < p = q < 2 and mn < (p ? 1)², the existence of at least one non‐extinction solution for any positive smooth initial data is proved. Copyright © 2015 John Wiley & Sons, Ltd.     

Krylov eigenvalue strategy using the FEAST algorithm with inexact system solves

The FEAST eigenvalue algorithm is a subspace iteration algorithm that uses contour integration to obtain the eigenvectors of a matrix for the eigenvalues that are located in any user‐defined region in the complex plane. By computing small numbers of eigenvalues in specific regions of the complex plane, FEAST is able to naturally parallelize the solution of eigenvalue problems by solving for multiple eigenpairs simultaneously. The traditional FEAST algorithm is implemented by directly solving collections of shifted linear systems of equations; in this paper, we describe a variation of the FEAST algorithm that uses iterative Krylov subspace algorithms for solving the shifted linear systems inexactly. We show that this iterative FEAST algorithm (which we call IFEAST) is mathematically equivalent to a block Krylov subspace method for solving eigenvalue problems. By using Krylov subspaces indirectly through solving shifted linear systems, rather than directly using them in projecting the eigenvalue problem, it becomes possible to use IFEAST to solve eigenvalue problems using very large dimension Krylov subspaces without ever having to store a basis for those subspaces. IFEAST thus combines the flexibility and power of Krylov methods, requiring only matrix–vector multiplication for solving eigenvalue problems, with the natural parallelism of the traditional FEAST algorithm. We discuss the relationship between IFEAST and more traditional Krylov methods and provide numerical examples illustrating its behavior.     

OpenMG: A new multigrid implementation in Python

In many large‐scale computations, systems of equations arise in the form Au = b, where A is a linear operation to be performed on the unknown data u, producing the known right‐hand side, b, which represents some constraint of known or assumed behavior of the system being modeled. Because such systems can be very large, solving them directly can be too slow. In contrast, a multigrid method removes different components of the error at different resolutions using smoothers that reduce high‐frequency components of the error more readily than low. Here, we present an open‐source multigrid solver written only in Python. OpenMG is a pure Python experimentation environment for testing multigrid concepts, not a production solver. The particular restriction method implemented is for ‘standard’ multigrid. By making the code simple and modular, we make the algorithmic details clear. The resulting solver is tested on an implicit pressure reservoir simulation problem with satisfactory results.Copyright © 2013 John Wiley & Sons, Ltd.     

Decay rate of the range component of solutions to some semilinear evolution equations

We examine the rate of decay to 0, as t → +∞., of the projection on the range of A of the solutions of an equation of the form u′ + Au + |u| ^p−1 u = 0 or u′′ + u′ + Au + |u| ^p−1 u = 0 in a bounded domain of ^N , where A = −Δ with Neumann boundary conditions or A = −Δ − λ₁ I with Dirichlet boundary conditions. In general this decay is much faster than the decay of the projection on the kernel; it is often exponential, but apparently not always.     

Partly zero eigenvectors and matrices that satisfy Au = b

For given vectors u,b ∈ F ⁿ (where F is a field with at least 3 elements), we establish criteria for deciding whether a digraph allows in its pattern class a matrix A which satisfies Au = b. As corollaries to this we give necessary and sufficient conditions for a pattern class to allow a matrix which has an eigenvector with a particular zero/nonzero pattern. Moreover we establish whether or not that eigenvector can correspond to a zero or nonzero eigenvalue. We use these results to establish the analogous results for loopfree digraphs, and thus we obtain results additional to work already done by Maybee, Olesky, and Van Den Driessche in this area. We then use bipartite graphs to generalize our criteria for solutions to Au = b to the rectangular case.     

Convergence conditions for a restarted GMRES method augmented with eigenspaces

We consider the GMRES(m,k) method for the solution of linear systems Ax=b, i.e. the restarted GMRES with restart m where to the standard Krylov subspace of dimension m the other subspace of dimension k is added, resulting in an augmented Krylov subspace. This additional subspace approximates usually an A‐invariant subspace. The eigenspaces associated with the eigenvalues closest to zero are commonly used, as those are thought to hinder convergence the most. The behaviour of residual bounds is described for various situations which can arise during the GMRES(m,k) process. The obtained estimates for the norm of the residual vector suggest sufficient conditions for convergence of GMRES(m,k) and illustrate that these augmentation techniques can remove stagnation of GMRES(m) in many cases. All estimates are independent of the choice of an initial approximation. Conclusions and remarks assessing numerically the quality of proposed bounds conclude the paper. Copyright © 2004 John Wiley & Sons, Ltd.     

Analyticity for Kuramoto–Sivashinsky‐type equations in two spatial dimensions

I. Stratis In this work, we investigate the analyticity properties of solutions of Kuramoto–Sivashinsky‐type equations in two spatial dimensions, with periodic initial data. In order to do this, we explore the applicability in three‐dimensional models of a spectral method, which was developed by the authors for the one‐dimensional Kuramoto–Sivashinsky equation. We introduce a criterion, which provides a sufficient condition for analyticity of a periodic function u∈C^∞, involving the rate of growth of ?ⁿu, in suitable norms, as n tends to infinity. This criterion allows us to establish spatial analyticity for the solutions of a variety of systems, including Topper–Kawahara, Frenkel–Indireshkumar, and Coward–Hall equations and their dispersively modified versions, once we assume that these systems possess global attractors. Copyright © 2015 John Wiley & Sons, Ltd.     

New convergence results and preconditioning strategies for the conjugate gradient method

Instead of the standard estimate in terms of the spectral condition number we develop a new CG iteration number estimate depending on the quantity B = 1/ntr M/(det M)^1/n, where M is an n × n preconditioned matrix. A new family of iterative methods for solving symmetric positive definite systems based on B-reducing strategies is described. Numerical results are presented for the new algorithms and compared with several well-known preconditioned CG methods.     

Spectral condition‐number estimation of large sparse matrices

We describe a randomized Krylov‐subspace method for estimating the spectral condition number of a real matrix A or indicating that it is numerically rank deficient. The main difficulty in estimating the condition number is the estimation of the smallest singular value ${\sigma }_{\min }$ of A. Our method estimates this value by solving a consistent linear least squares problem with a known solution using a specific Krylov‐subspace method called LSQR. In this method, the forward error tends to concentrate in the direction of a right singular vector corresponding to ${\sigma }_{\min }$. Extensive experiments show that the method is able to estimate well the condition number of a wide array of matrices. It can sometimes estimate the condition number when running dense singular value decomposition would be impractical due to the computational cost or the memory requirements. The method uses very little memory (it inherits this property from LSQR), and it works equally well on square and rectangular matrices.