The conditioning of least‐squares problems in variational data assimilation

In variational data assimilation a least‐squares objective function is minimised to obtain the most likely state of a dynamical system. This objective function combines observation and prior (or background) data weighted by their respective error statistics. In numerical weather prediction, data assimilation is used to estimate the current atmospheric state, which then serves as an initial condition for a forecast. New developments in the treatment of observation uncertainties have recently been shown to cause convergence problems for this least‐squares minimisation. This is important for operational numerical weather prediction centres due to the time constraints of producing regular forecasts. The condition number of the Hessian of the objective function can be used as a proxy to investigate the speed of convergence of the least‐squares minimisation. In this paper we develop novel theoretical bounds on the condition number of the Hessian. These new bounds depend on the minimum eigenvalue of the observation error covariance matrix and the ratio of background error variance to observation error variance. Numerical tests in a linear setting show that the location of observation measurements has an important effect on the condition number of the Hessian. We identify that the conditioning of the problem is related to the complex interactions between observation error covariance and background error covariance matrices. Increased understanding of the role of each constituent matrix in the conditioning of the Hessian will prove useful for informing the choice of correlated observation error covariance matrix and observation location, particularly for practical applications.     

Matrix-free interior point method

In this paper we present a redesign of a linear algebra kernel of an interior point method to avoid the explicit use of problem matrices. The only access to the original problem data needed are the matrix-vector multiplications with the Hessian and Jacobian matrices. Such a redesign requires the use of suitably preconditioned iterative methods and imposes restrictions on the way the preconditioner is computed. A two-step approach is used to design a preconditioner. First, the Newton equation system is regularized to guarantee better numerical properties and then it is preconditioned. The preconditioner is implicit, that is, its computation requires only matrix-vector multiplications with the original problem data. The method is therefore well-suited to problems in which matrices are not explicitly available and/or are too large to be stored in computer memory. Numerical properties of the approach are studied including the analysis of the conditioning of the regularized system and that of the preconditioned regularized system. The method has been implemented and preliminary computational results for small problems limited to 1 million of variables and 10 million of nonzero elements demonstrate the feasibility of the approach.     

A Study on the Conditioning of Finite Element Equations with Arbitrary Anisotropic Meshes via a Density Function Approach

The linear finite element approximation of a general linear diffusion problem with arbitrary anisotropic meshes is considered. The conditioning of the resultant stiffness matrix and the Jacobi preconditioned stiffness matrix is investigated using a density function approach proposed by Fried in 1973. It is shown that the approach can be made mathematically rigorous for general domains and used to develop bounds on the smallest eigenvalue and the condition number that are sharper than existing estimates in one and two dimensions and comparable in three and higher dimensions. The new results reveal that the mesh concentration near the boundary has less influence on the condition number than the mesh concentration in the interior of the domain. This is especially true for the Jacobi preconditioned system where the former has little or almost no influence on the condition number. Numerical examples are presented.     

Convergence acceleration of preconditioned conjugate gradient solver based on error vector sampling for a sequence of linear systems

In this article, we focus on solving a sequence of linear systems that have identical (or similar) coefficient matrices. For this type of problem, we investigate subspace correction (SC) and deflation methods, which use an auxiliary matrix (subspace) to accelerate the convergence of the iterative method. In practical simulations, these acceleration methods typically work well when the range of the auxiliary matrix contains eigenspaces corresponding to small eigenvalues of the coefficient matrix. We develop a new algebraic auxiliary matrix construction method based on error vector sampling in which eigenvectors with small eigenvalues are efficiently identified in the solution process. We use the generated auxiliary matrix for convergence acceleration in the following solution step. Numerical tests confirm that both SC and deflation methods with the auxiliary matrix can accelerate the solution process of the iterative solver. Furthermore, we examine the applicability of our technique to the estimation of the condition number of the coefficient matrix. We also present the algorithm of the preconditioned conjugate gradient method with condition number estimation.     

Minimizing the condition number of a positive definite matrix by completion

Summary. We consider the problem of minimizing the spectral condition number of a positive definite matrix by completion: \noindent where is an Hermitian positive definite matrix, a matrix and is a free Hermitian matrix. We reduce this problem to an optimization problem for a convex function in one variable. Using the minimal solution of this problem we characterize the complete set of matrices that give the minimum condition number. Received October 15, 1993     

An example of numerical nonconvergence of a variable-metric method

An example is given where truncation error, caused by finite computer arithmetic, leads to the BFGS variable-metric method becoming stuck, despite the approximated Hessian matrix, the gradient vector, and the search vector satisfying analytical conditions for convergence. A restart criterion to eliminate the problem is suggested.     

Structure of the Hessian matrix and an economical implementation of Newton’s method in the problem of canonical approximation of tensors

A tensor given by its canonical decomposition is approximated by another tensor (again, in the canonical decomposition) of fixed lower rank. For this problem, the structure of the Hessian matrix of the objective function is analyzed. It is shown that all the auxiliary matrices needed for constructing the quadratic model can be calculated so that the computational effort is a quadratic function of the tensor dimensionality (rather than a cubic function as in earlier publications). An economical version of the trust region Newton method is proposed in which the structure of the Hessian matrix is efficiently used for multiplying this matrix by vectors and for scaling the trust region. At each step, the subproblem of minimizing the quadratic model in the trust region is solved using the preconditioned conjugate gradient method, which is terminated if a negative curvature direction is detected for the Hessian matrix.     

Eigenvalue estimates for preconditioned saddle point matrices

New accurate eigenvalue bounds for symmetric matrices of saddle point form are derived and applied for both unpreconditioned and preconditioned versions of the matrices. The estimates enable a better understanding of how preconditioners should be chosen. The preconditioners provide efficient iterative solution of the corresponding linear systems with, for some important applications, an optimal order of computational complexity. The methods are applied for Stokes problem and for linear elasticity problems. Copyright © 2005 John Wiley & Sons, Ltd.     

Preconditioned iterative methods for a class of nonlinear eigenvalue problems

This paper proposes new iterative methods for the efficient computation of the smallest eigenvalue of symmetric nonlinear matrix eigenvalue problems of large order with a monotone dependence on the spectral parameter. Monotone nonlinear eigenvalue problems for differential equations have important applications in mechanics and physics. The discretization of these eigenvalue problems leads to nonlinear eigenvalue problems with very large sparse ill-conditioned matrices monotonically depending on the spectral parameter. To compute the smallest eigenvalue of large-scale matrix nonlinear eigenvalue problems, we suggest preconditioned iterative methods: preconditioned simple iteration method, preconditioned steepest descent method, and preconditioned conjugate gradient method. These methods use only matrix-vector multiplications, preconditioner-vector multiplications, linear operations with vectors, and inner products of vectors. We investigate the convergence and derive grid-independent error estimates for these methods. Numerical experiments demonstrate the practical effectiveness of the proposed methods for a model problem.     

A survey of some estimates of eigenvalues and condition numbers for certain preconditioned matrices

Eigenvalue and condition number estimates for preconditioned iteration matrices provide the information required to estimate the rate of convergence of iterative methods, such as preconditioned conjugate gradient methods. In recent years various estimates have been derived for (perturbed) modified (block) incomplete factorizations. We survey and extend some of these and derive new estimates. In particular we derive upper and lower estimates of individual eigenvalues and of condition number. This includes a discussion that the condition number of preconditioned second order elliptic difference matrices is O(h⁻¹). Some of the methods are applied to compute certain parameters involved in the computation of the preconditioner.     

广义块Toeplitz特征值问题的基于sine变换的预处理子

在求块Toeplitz矩阵束（Amn，Bmn）特征值的Lanczos过程中，通过对移位块Toepltz矩阵Amn-ρBmn进行基于sine变换的块预处理，从而改进了位移块Toeplitz矩阵的谱分布，加速了Lanczos过程的收敛速度．该块预处理方法能通过快速算法有效快速执行．本文证明了预处理后Lanczos过程收敛迅速，并通过实验证明该算法求解大规模矩阵问题尤其有效．     

Stability of the Optimal Solution to a Problem of Variational Assimilation with Error Covariance Matrices of Observational Data for a Sea Thermodynamics Model

A mathematical model of sea thermodynamics developed at the Institute of Numerical Mathematics of the Russian Academy of Sciences is considered. A problem of variational assimilation of daily-averaged sea surface temperature (SST) data with observational error covariance matrices is formulated and investigated. On the basis of variational assimilation of satellite observation data, an inverse problem of restoring the heat flux on the sea surface is solved. The stability of the optimal solution to the problem of variational data assimilation is studied, and the results of numerical experiments with the model for Baltic Sea dynamics are presented.     

An incomplete Hessian Newton minimization method and its application in a chemical database problem

To efficiently solve a large scale unconstrained minimization problem with a dense Hessian matrix, this paper proposes to use an incomplete Hessian matrix to define a new modified Newton method, called the incomplete Hessian Newton method (IHN). A theoretical analysis shows that IHN is convergent globally, and has a linear rate of convergence with a properly selected symmetric, positive definite incomplete Hessian matrix. It also shows that the Wolfe conditions hold in IHN with a line search step length of one. As an important application, an effective IHN and a modified IHN, called the truncated-IHN method (T-IHN), are constructed for solving a large scale chemical database optimal projection mapping problem. T-IHN is shown to work well even with indefinite incomplete Hessian matrices. Numerical results confirm the theoretical results of IHN, and demonstrate the promising potential of T-IHN as an efficient minimization algorithm.     

Convergence of Rump's method for inverting arbitrarily ill-conditioned matrices

In this paper, the problem of inverting regular matrices with arbitrarily large condition number is treated in double precision defined by IEEE 754 floating point standard. In about 1984, Rump derived a method for inverting arbitrarily ill-conditioned matrices. The method requires the possibility to calculate a dot product in higher precision. Rump's method is of theoretical interest. Rump made it clear that inverting an arbitrarily ill-conditioned matrix in single or double precision does not produce meaningless numbers, but contains a lot of information in it. Rump's method uses such inverses as preconditioners. Numerical experiments exhibit that Rump's method converges rapidly for various matrices with large condition numbers. Why Rump's method is so efficient for inverting arbitrarily ill-conditioned matrices is a little mysterious. Thus, to prove its convergence is an interesting problem in numerical error analysis. In this article, a convergence theorem is presented for a variant of Rump's method.     

基于子结构法构造用非协调元解椭圆型问题的预处理器（Ⅰ）

基于子结构法构造用非协调元解椭圆型问题的预处理器（Ⅰ）顾金生，胡显承（清华大学应用数学系）ＴＨＥＣＯＮＳＴＲＵＣＴＩＯＮＯＦＰＲＥＣＯＮＤＩＴＩＯＮＥＲＳＦＯＲＥＬＬＩＰＴＩＣＰＲＯＢＬＥＭＳＤＩＳＣＲＥＴＩＺＥＤＢＹＮＯＮＣＯＮＦＯＲＭＩＮＧＦＩＮ...     

Convergence of quasi-Newton matrices generated by the symmetric rank one update

Quasi-Newton algorithms for unconstrained nonlinear minimization generate a sequence of matrices that can be considered as approximations of the objective function second derivatives. This paper gives conditions under which these approximations can be proved to converge globally to the true Hessian matrix, in the case where the Symmetric Rank One update formula is used. The rate of convergence is also examined and proven to be improving with the rate of convergence of the underlying iterates. The theory is confirmed by some numerical experiments that also show the convergence of the Hessian approximations to be substantially slower for other known quasi-Newton formulae.The work of this author was supported by the National Sciences and Engineering Research Council of Canada, and by the Information Technology Research Centre, which is funded by the Province of Ontario.     

Spectral estimates for unreduced symmetric KKT systems arising from Interior Point methods

We consider symmetrized Karush–Kuhn–Tucker systems arising in the solution of convex quadratic programming problems in standard form by Interior Point methods. Their coefficient matrices usually have 3 × 3 block structure, and under suitable conditions on both the quadratic programming problem and the solution, they are nonsingular in the limit. We present new spectral estimates for these matrices: the new bounds are established for the unpreconditioned matrices and for the matrices preconditioned by symmetric positive definite augmented preconditioners. Some of the obtained results complete the analysis recently given by Greif, Moulding, and Orban in [SIAM J. Optim., 24 (2014), pp. 49‐83]. The sharpness of the new estimates is illustrated by numerical experiments. Copyright © 2016 John Wiley & Sons, Ltd.     

The Global Convergence of Self-Scaling BFGS Algorithm with Nonmonotone Line Search for Unconstrained Nonconvex Optimization Problems

The self-scaling quasi-Newton method solves an unconstrained optimization problem by scaling the Hessian approximation matrix before it is updated at each iteration to avoid the possible large eigenvalues in the Hessian approximation matrices of the objective function. It has been proved in the literature that this method has the global and superlinear convergence when the objective function is convex （or even uniformly convex）. We propose to solve unconstrained nonconvex optimization problems by a self-scaling BFGS algorithm with nonmonotone linear search. Nonmonotone line search has been recognized in numerical practices as a competitive approach for solving large-scale nonlinear problems. We consider two different nonmonotone line search forms and study the global convergence of these nonmonotone self-scale BFGS algorithms. We prove that, under some weaker condition than that in the literature, both forms of the self-scaling BFGS algorithm are globally convergent for unconstrained nonconvex optimization problems.     

EXPLICIT BOUNDS OF EIGENVALUES FOR STIFFNESS MATRICES BY QUADRATIC HIERARCHICAL BASIS METHOD)

The bounds for the eigenvalues of the stiffness matrices in the finite element discretiza-tion corresponding to Lu:=-u″with zero boundary conditions by quadratic hierarchical basic are shown explicitly.The condition numberof the resulting system behaves like O(1/h)where h is the mesh size.We also analyze a main diagonal preconditioner of the stiffness matrix which reduces the condition number of the preconditioned system to O(1).     

Numerical study in stochastic homogenization for elliptic partial differential equations: Convergence rate in the size of representative volume elements

We describe the numerical scheme for the discretization and solution of 2D elliptic equations with strongly varying piecewise constant coefficients arising in the stochastic homogenization of multiscale composite materials. An efficient stiffness matrix generation scheme based on assembling the local Kronecker product matrices is introduced. The resulting large linear systems of equations are solved by the preconditioned conjugate gradient iteration with a convergence rate that is independent of the grid size and the variation in jumping coefficients (contrast). Using this solver, we numerically investigate the convergence of the representative volume element (RVE) method in stochastic homogenization that extracts the effective behavior of the random coefficient field. Our numerical experiments confirm the asymptotic convergence rate of systematic error and standard deviation in the size of RVE rigorously established in Gloria et al. The asymptotic behavior of covariances of the homogenized matrix in the form of a quartic tensor is also studied numerically. Our approach allows laptop computation of sufficiently large number of stochastic realizations even for large sizes of the RVE.