Computing the conditioning of the components of a linear least‐squares solution

In this paper, we address the accuracy of the results for the overdetermined full rank linear least‐squares problem. We recall theoretical results obtained in (SIAM J. Matrix Anal. Appl. 2007; 29 (2):413–433) on conditioning of the least‐squares solution and the components of the solution when the matrix perturbations are measured in Frobenius or spectral norms. Then we define computable estimates for these condition numbers and we interpret them in terms of statistical quantities when the regression matrix and the right‐hand side are perturbed. In particular, we show that in the classical linear statistical model, the ratio of the variance of one component of the solution by the variance of the right‐hand side is exactly the condition number of this solution component when only perturbations on the right‐hand side are considered. We explain how to compute the variance–covariance matrix and the least‐squares conditioning using the libraries LAPACK (LAPACK Users' Guide (3rd edn). SIAM: Philadelphia, 1999) and ScaLAPACK (ScaLAPACK Users' Guide. SIAM: Philadelphia, 1997) and we give the corresponding computational cost. Finally we present a small historical numerical example that was used by Laplace (Théorie Analytique des Probabilités. Mme Ve Courcier, 1820; 497–530) for computing the mass of Jupiter and a physical application if the area of space geodesy. Copyright © 2008 John Wiley & Sons, Ltd.     

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.     

Improved seed methods for symmetric positive definite linear equations with multiple right‐hand sides

We consider symmetric positive definite systems of linear equations with multiple right‐hand sides. The seed conjugate gradient (CG) method solves one right‐hand side with the CG method and simultaneously projects over the Krylov subspace thus developed for the other right‐hand sides. Then the next system is solved and used to seed the remaining ones. Rounding error in the CG method limits how much the seeding can improve convergence. We propose three changes to the seed CG method: only the first right‐hand side is used for seeding, this system is solved past convergence, and the roundoff error is controlled with some reorthogonalization. We will show that results are actually better with only one seeding, even in the case of related right‐hand sides. Controlling rounding error gives the potential for rapid convergence for the second and subsequent right‐hand sides. Polynomial preconditioning can help reduce storage needed for reorthogonalization. The new seed methods are applied to examples including matrices from quantum chromodynamics. Copyright © 2013 John Wiley & Sons, Ltd.     

Local mesh refinement for the discretization of Neumann boundary control problems on polyhedra

This paper deals with optimal control problems constrained by linear elliptic partial differential equations. The case where the right‐hand side of the Neumann boundary is controlled, is studied. The variational discretization concept for these problems is applied, and discretization error estimates are derived. On polyhedral domains, one has to deal with edge and corner singularities, which reduce the convergence rate of the discrete solutions, that is, one cannot expect convergence order two for linear finite elements on quasi‐uniform meshes in general. As a remedy, a local mesh refinement strategy is presented, and a priori bounds for the refinement parameters are derived such that convergence with optimal rate is guaranteed. As a by‐product, finite element error estimates in the H¹(Ω)‐norm, L²(Ω)‐norm and L²(Γ)‐norm for the boundary value problem are obtained, where the latter one turned out to be the main challenge. Copyright © 2015 John Wiley & Sons, Ltd.     

Optimal damping of selected eigenfrequencies using dimension reduction

We consider linear vibrational systems described by a system of second‐order differential equations of the form , where M and K are positive definite matrices, representing mass and stiffness, respectively. The damping matrix D is assumed to be positive semidefinite. We are interested in finding an optimal damping matrix that will damp a certain (critical) part of the eigenfrequencies. For this, we use an optimization criterion based on the minimization of the average total energy of the system. This is equivalent to the minimization of the trace of the solution of the corresponding Lyapunov equation AX + XA^T = ?GG^T, where A is the matrix obtained from linearizing the second‐order differential equation, and G depends on the critical part of the eigenfrequencies to be damped. The main result is the efficient approximation and the corresponding error bound for the trace of the solution of the Lyapunov equation obtained through dimension reduction, which includes the influence of the right‐hand side GG^T and allows us to control the accuracy of the trace approximation. This trace approximation yields a very accelerated optimization algorithm for determining the optimal damping. Copyright © 2011 John Wiley & Sons, Ltd.     

Solution to fuzzy system of linear equations with crisp coefficients

This paper presents a new and simple method to solve fuzzy real system of linear equations by solving two n × n crisp systems of linear equations. In an original system, the coefficient matrix is considered as real crisp, whereas an unknown variable vector and right hand side vector are considered as fuzzy. The general system is initially solved by adding and subtracting the left and right bounds of the vectors respectively. Then obtained solutions are used to get a final solution of the original system. The proposed method is used to solve five example problems. The results obtained are also compared with the known solutions and found to be in good agreement with them.     

Delta and singular delta locus for one‐dimensional systems of conservation laws

This work gives a condition for existence of singular and delta shock wave solutions to Riemann problem for 2×2 systems of conservation laws. For a fixed left‐hand side value of Riemann data, the condition obtained in the paper describes a set of possible right‐hand side values. The procedure is similar to the standard one of finding the Hugoniot locus. Fluxes of the considered systems are globally Lipschitz with respect to one of the dependent variables. The association in a Colombeau‐type algebra is used as a solution concept. Copyright © 2004 John Wiley &Sons, Ltd.     

A uniformly optimal‐order estimate for finite volume method for advection‐diffusion equations

We prove an optimal‐order error estimate in a weighted energy norm for finite volume method for two‐dimensional time‐dependent advection–diffusion equations on a uniform space‐time partition of the domain. The generic constants in the estimates depend only on certain norms of the true solution but not on the scaling parameter. These estimates, combined with a priori stability estimates of the governing partial differential equations with full regularity, yield a uniform estimate of the finite volume method, in which the generic constants depend only on the Sobolev norms of the initial and right side data but not on the scaling parameter. We use the interpolation of spaces and stability estimates to derive a uniform estimate for problems with minimal or intermediate regularity, where the convergence rates are proportional to certain Besov norms of the initial and right‐hand side data. © 2013 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 30: 17‐43, 2014     

Extending the eigCG algorithm to nonsymmetric Lanczos for linear systems with multiple right‐hand sides

The technique that was used to build the eigCG algorithm for sparse symmetric linear systems is extended to the nonsymmetric case using the BiCG algorithm. We show that, similar to the symmetric case, we can build an algorithm that is capable of computing a few smallest magnitude eigenvalues and their corresponding left and right eigenvectors of a nonsymmetric matrix using only a small window of the BiCG residuals while simultaneously solving a linear system with that matrix. For a system with multiple right‐hand sides, we give an algorithm that computes incrementally more eigenvalues while solving the first few systems and then uses the computed eigenvectors to deflate BiCGStab for the remaining systems. Our experiments on various test problems, including Lattice QCD, show the remarkable ability of eigBiCG to compute spectral approximations with accuracy comparable with that of the unrestarted, nonsymmetric Lanczos. Furthermore, our incremental eigBiCG followed by appropriately restarted and deflated BiCGStab provides a competitive method for systems with multiple right‐hand sides. Copyright © 2013 John Wiley & Sons, Ltd.     

模糊线性系统的扰动分析

使用谱范数分析了模糊线性系统在三种情形下的扰动: (1)右端模糊向量有扰动, 系数矩阵不变; (2)系数矩阵有扰动,右端模糊向量不变; (3)系数矩阵和右端模糊向量都有扰动,并通过数值实例验证给出的扰动界的估计.     

A new error in variables model for solving positive definite linear system using orthogonal matrix decompositions

The need to estimate a positive definite solution to an overdetermined linear system of equations with multiple right hand side vectors arises in several process control contexts. The coefficient and the right hand side matrices are respectively named data and target matrices. A number of optimization methods were proposed for solving such problems, in which the data matrix is unrealistically assumed to be error free. Here, considering error in measured data and target matrices, we present an approach to solve a positive definite constrained linear system of equations based on the use of a newly defined error function. To minimize the defined error function, we derive necessary and sufficient optimality conditions and outline a direct algorithm to compute the solution. We provide a comparison of our proposed approach and two existing methods, the interior point method and a method based on quadratic programming. Two important characteristics of our proposed method as compared to the existing methods are computing the solution directly and considering error both in data and target matrices. Moreover, numerical test results show that the new approach leads to smaller standard deviations of error entries and smaller effective rank as desired by control problems. Furthermore, in a comparative study, using the Dolan-Moré performance profiles, we show the approach to be more efficient.     

Pointwise Estimates for Laplace Equation. Applications to the Free Boundary of the Obstacle Problem with Dini Coefficients

In this paper we are interested in pointwise regularity of solutions to elliptic equations. In a first result, we prove that if the modulus of mean oscillation of Δu at the origin is Dini (in L ^p average), then the origin is a Lebesgue point of continuity (still in L ^p average) for the second derivatives D ² u. We extend this pointwise regularity result to the obstacle problem for the Laplace equation with Dini right hand side at the origin. Under these assumptions, we prove that the solution to the obstacle problem has a Taylor expansion up to the order 2 (in the L ^p average). Moreover we get a quantitative estimate of the error in this Taylor expansion for regular points of the free boundary. In the case where the right hand side is moreover double Dini at the origin, we also get a quantitative estimate of the error for singular points of the free boundary. Our method of proof is based on some decay estimates obtained by contradiction, using blow-up arguments and Liouville Theorems. In the case of singular points, our method uses moreover a refined monotonicity formula.      

A note on the efficiency of the conjugate gradient method for a class of time‐dependent problems

We discuss the efficiency of the conjugate gradient (CG) method for solving a sequence of linear systems; Auⁿ⁺¹ = uⁿ, where A is assumed to be sparse, symmetric, and positive definite. We show that under certain conditions the Krylov subspace, which is generated when solving the first linear system Au¹ = u⁰, contains the solutions {uⁿ} for subsequent time steps. The solutions of these equations can therefore be computed by a straightforward projection of the right‐hand side onto the already computed Krylov subspace. Our theoretical considerations are illustrated by numerical experiments that compare this method with the order‐optimal scheme obtained by applying the multigrid method as a preconditioner for the CG‐method at each time step. Copyright © 2007 John Wiley & Sons, Ltd.     

The mapping properties of the radiosity operator along an edge

In this article we study the radiosity operator along an edge between two adjacent half‐planes. First we show that the radiosity operator is invertible in a whole scale of anisotropic Sobolev spaces. In the absence of any shadows we are able to derive regularity properties of the solution, which depend only on the angle between the half‐planes, the reflectivity coefficients and the right‐hand side. This work can be considered as a supplement to the article of Rathsfeld (Mathematical Methods in the Applied Sciences 1999; 22 : 217–241). Copyright © 2002 John Wiley & Sons, Ltd.     

Limited memory preconditioners for symmetric indefinite problems with application to structural mechanics

This paper presents a class of limited memory preconditioners (LMP) for solving linear systems of equations with symmetric indefinite matrices and multiple right‐hand sides. These preconditioners based on limited memory quasi‐Newton formulas require a small number k of linearly independent vectors and may be used to improve an existing first‐level preconditioner. The contributions of the paper are threefold. First, we derive a formula to characterize the spectrum of the preconditioned operator. A spectral analysis of the preconditioned matrix shows that the eigenvalues are all real and that the LMP class is able to cluster at least k eigenvalues at 1. Secondly, we show that the eigenvalues of the preconditioned matrix enjoy interlacing properties with respect to the eigenvalues of the original matrix provided that the k linearly independent vectors have been prior projected onto the invariant subspaces associated with the eigenvalues of the original matrix in the open right and left half‐plane, respectively. Third, we focus on theoretical properties of the Ritz‐LMP variant, where Ritz information is used to determine the k vectors. Finally, we illustrate the numerical behaviour of the Ritz limited memory preconditioners on realistic applications in structural mechanics that require the solution of sequences of large‐scale symmetric saddle‐point systems. Numerical experiments show the relevance of the proposed preconditioner leading to a significant decrease in terms of computational operations when solving such sequences of linear systems. A saving of up to 43% in terms of computational effort is obtained on one of these applications. Copyright © 2016 John Wiley & Sons, Ltd.     

Right preconditioned MINRES for singular systems†

We consider solving large sparse symmetric singular linear systems. We first introduce an algorithm for right preconditioned minimum residual (MINRES) and prove that its iterates converge to the preconditioner weighted least squares solution without breakdown for an arbitrary right‐hand‐side vector and an arbitrary initial vector even if the linear system is singular and inconsistent. For the special case when the system is consistent, we prove that the iterates converge to a min‐norm solution with respect to the preconditioner if the initial vector is in the range space of the right preconditioned coefficient matrix. Furthermore, we propose a right preconditioned MINRES using symmetric successive over‐relaxation (SSOR) with Eisenstat's trick. Some numerical experiments on semidefinite systems in electromagnetic analysis and so forth indicate that the method is efficient and robust. Finally, we show that the residual norm can be further reduced by restarting the iterations.     

On solution of total least squares problems with multiple right-hand sides

Consider a linear approximation problem AX≈B with multiple right–hand sides. When errors in the data are confirmed both to B and A, the total least squares (TLS) concept is used to solve this problem. Contrary to the standard least squares approximation problem, a solution of the TLS problem may not exist. For a single (vector) right–hand side, the classical theory has been developed by G.H. Golub, C.F. Van Loan [2], and S. Van Huffel, J. Vandewalle [4], and then complemented recently by the core problem approach of C.C. Paige, Z. Strakoš [5,6,7]. Analysis of the problem with multiple right–hand sides is still under development. In this short contribution we present conditions for the existence of a TLS solution. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Reliable,goal‐oriented postprocessing for FE‐discretizations

In this article we present strategies to improve the quality of adaptive FE‐approximations measured in terms of linear functionals. The ideas are based on the so called dual‐weighted‐residual (DWR) approach to a posteriori error control for FE‐schemes. In more details, we exploit those parts of an underlying error representation, which are completely computable, to improve the FE‐solution. Furthermore, the remaining parts of the error identity can be estimated by well‐established a posteriori energy estimates yielding reliable error bounds for the postprocessed values. © 2004 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2005     

Spectral collocation methods for nonlinear weakly singular Volterra integro‐differential equations

In this paper, a spectral collocation approximation is proposed for neutral and nonlinear weakly singular Volterra integro‐differential equations (VIDEs) with non‐smooth solutions. We use some suitable variable transformations to change the original equation into a new equation, so that the solution of the resulting equation possesses better regularity, and the the Jacobi orthogonal polynomial theory can be applied conveniently. Under reasonable assumptions on the nonlinearity, we carry out a rigorous error analysis in L^∞ norm and weighted L² norm. To perform the numerical simulations, some test examples (linear and nonlinear) are considered with nonsmooth solutions, and numerical results are presented. Further more, the comparative study of the proposed methods with some existing numerical methods is provided.     

Die mehrstellenformeln für den Laplaceoperator

Summary Difference methods for the numerical solution of linear partial differential equations may often be improved by using a weighted right hand side instead of the original right hand side of the differential equation. Difference formulas, for which that is possible, are called Mehrstellenformeln or Hermitian formulas. In this paper the Hermitian formulas for the approximation of Laplace's operator are characterized by a very simple condition. We prove, that in two-dimensional case for a Hermitian formula of ordern at leastn+3 discretization points are necessary. We give examples of such optimal formulas of arbitrary high-order.