Newton Matrix Research Articles

Jacobian-free implicit MDRK methods for stiff systems of ODEs

In this work, an approximate family of implicit multiderivative Runge-Kutta (MDRK) time integrators for stiff initial value problems is presented. The approximation procedure is based on the recent Approximate Implicit Taylor method (Baeza et al., 2020 [7]). As a Taylor method can be written in MDRK format, the novel family constitutes a multistage generalization. Two different alternatives are investigated for the computation of the higher order derivatives: either directly as part of the stage equation, or either as a separate formula for each derivative added on top of the stage equation itself. From linearizing through Newton's method, it turns out that the conditioning of the Newton matrix behaves significantly different for both cases. We show that direct computation results in a matrix with a conditioning that is highly dependent on the stiffness, increasing exponentially in the stiffness parameter with the amount of derivatives. Adding separate formulas has a more favorable behavior, the matrix conditioning being linearly dependent on the stiffness, regardless of the amount of derivatives. Despite increasing the Newton system significantly in size, through several numerical results it is demonstrated that doing so can be considerably beneficial.

Efficient and practical Newton solvers for non-linear Stokes systems in geodynamic problems

Many problems in geodynamic modelling result in a non-linear Stokes problem in which the viscosity depends on the strain rate and pressure (in addition to other variables). After discretization, the resulting non-linear system is most commonly solved using a Picard fixedpoint iteration. However, it is well understood that Newton's method - when augmented by globalization strategies to ensure convergence even from points far from the solution - can be substantially more efficient and accurate than a Picard solver. In this contribution, we evaluate how a straightforward Newton method must be modified to allow for the kinds of rheologies common in geodynamics. Specifically, we show that the Newton step is not actually well posed for strain rate-weakening models without modifications to the Newton matrix. We derive modifications that guarantee well posedness and that also allow for efficient solution strategies by ensuring that the top left block of the Newton matrix is symmetric and positive definite.We demonstrate the applicability and relevance of thesemodifications with a sequence of benchmarks and a test case of realistic complexity.

Globally Convergent Levenberg-Marquardt Method for Phase Retrieval

In this paper, we consider a nonlinear least squares model for the phase retrieval problem. Since the Hessian matrix may not be positive definite and the Gauss–Newton (GN) matrix is singular at any optimal solution, we propose a modified Levenberg–Marquardt (LM) method, where the Hessian is substituted by a summation of the GN matrix and a regularization term. Similar to the well-known Wirtinger flow (WF) algorithm under certain assumptions, we start from an initial point provably close to the set of the global optimal solutions. Global linear convergence and local quadratic convergence to the global solution set are established by estimating the smallest nonzero eigenvalues of the GN matrix, establishing local error bound properties and constructing a modified regularization condition. The computational cost becomes tractable if a preconditioned conjugate gradient method is applied to solve the LM equation inexactly. Specifically, the pre-conditioner is constructed from the expectation of the LM coefficient matrix by assuming the independence between the measurements and iteration point. The initialization can be dropped by slightly modifying the algorithm. Preliminary numerical experiments show that our algorithm is robust and it is often faster than the WF method on both random examples and natural image recovery.

Line search globalization of a semismooth Newton method for operator equations in Hilbert spaces with applications in optimal control

We consider the numerical solution of nonlinear and nonsmooth operator equations in Hilbert spaces. A semismooth Newton method is used for search direction generation. The operator equation is solved by a globalized semismooth Newton method that is equipped with an Armijo linesearch using a semismooth merit function. We prove that an accumulation point of the globalized algorithm is a solution and transition to fast local convergence under a directional Hadamard-like continuity assumption on the Newton matrix. In particular, no auxiliary descent directions or smoothing steps are required. Finally, we apply this method to a control-constrained and also to a regularized state-constrained optimal control problem subject to partial differential equations.

Solving the ghost–gluon system of Yang–Mills theory on GPUs

We solve the ghost–gluon system of Yang–Mills theory using graphics processing units (GPUs). Working in the Landau gauge, we use the Dyson–Schwinger formalism for the mathematical description, as this approach is well suited to directly benefit from the computing power of the GPUs. With the help of a Chebyshev expansion for the dressing functions and a subsequent appliance of a Newton–Raphson method, the non-linear system of coupled integral equations is linearized. The resulting Newton matrix is generated in parallel using OpenMPI and CUDA™. Our results show that it is possible to cut down the run time by two orders of magnitude as compared to a sequential version of the code. This makes the proposed techniques well suited for Dyson–Schwinger calculations on more complicated systems where the Yang–Mills sector of QCD serves as a starting point. In addition, the computation of Schwinger functions using GPU devices is studied.

Matrices and spectra satisfying the Newton inequalities

An n -by- n real matrix is called a Newton matrix (and its eigenvalues a Newton spectrum) if the normalized coefficients of its characteristic polynomial satisfy the Newton inequalities. A number of basic observations are made about Newton matrices, including closure under inversion, and then it is shown that a Newton matrix with nonnegative coefficients remains Newton under right translations. Those matrices that become (and stay) Newton under translation are characterized. In particular, Newton spectra in low dimensions are characterized.

On the solution of large-scale SDP problems by the modified barrier method using iterative solvers

The limiting factors of second-order methods for large-scale semidefinite optimization are the storage and factorization of the Newton matrix. For a particular algorithm based on the modified barri...

Efficient solution techniques for implicit finite element schemes with flux limiters

AbstractThe algebraic flux correction (AFC) paradigm is equipped with efficient solution strategies for implicit time‐stepping schemes. It is shown that Newton‐like techniques can be applied to the nonlinear systems of equations resulting from the application of high‐resolution flux limiting schemes. To this end, the Jacobian matrix is approximated by means of first‐ or second‐order finite differences. The edge‐based formulation of AFC schemes can be exploited to devise an efficient assembly procedure for the Jacobian. Each matrix entry is constructed from a differential and an average contribution edge by edge. The perturbation of solution values affects the nodal correction factors at neighbouring vertices so that the stencil for each individual node needs to be extended. Two alternative strategies for constructing the corresponding sparsity pattern of the resulting Jacobian are proposed. For nonlinear governing equations, the contribution to the Newton matrix which is associated with the discrete transport operator is approximated by means of divided differences and assembled edge by edge. Numerical examples for both linear and nonlinear benchmark problems are presented to illustrate the superiority of Newton methods as compared to the standard defect correction approach. Copyright © 2007 John Wiley & Sons, Ltd.

On the solution of large-scale SDP problems by the modified barrier method using iterative solvers

The limiting factors of second-order methods for large-scale semidefinite optimization are the storage and factorization of the Newton matrix. For a particular algorithm based on the modified barrier method, we propose to use iterative solvers instead of the routinely used direct factorization techniques. The preconditioned conjugate gradient method proves to be a viable alternative for problems with a large number of variables and modest size of the constrained matrix. We further propose to avoid explicit calculation of the Newton matrix either by an implicit scheme in the matrix---vector product or using a finite-difference formula. This leads to huge savings in memory requirements and, for certain problems, to further speed-up of the algorithm.

Armijo Newton method for convex best interpolation

Newton's method with Armijo line search (Armijo Newton method) has been practically known extremely efficient for the problem of convex best interpolation and numerical experiment strongly indicates its global convergence. However, similar to the classical Newton method, the Newton matrix far from the solution may be singular or near singular, posing a great deal of difficulties in proving the global convergence of Armijo Newton method. By employing the objective function of Lagrange dual problem, it is observed that whenever the Newton matrix is near singular at some point, one can easily [at cost of O(N), N is the dimension of the problem] find a nearby point which has a well-conditioned Newton matrix and lower function value. In this case, Armijo Newton method starts from this nearby point. We prove that this slightly modified Armijo Newton method is globally as well as locally quadratically convergent, and in an important case, it also has finite termination property. Numerical results demonstrate the efficiency of the proposed method.

Symbolic and Numeric Methods for Exploiting Structure in Constructing Resultant Matrices

Resultants characterize the existence of roots of systems of multivariate nonlinear polynomial equations, while their matrices reduce the computation of all common zeros to a problem in linear algebra. Sparse elimination theory has introduced the sparse (or toric) resultant, which takes into account the sparse structure of the polynomials. The construction of sparse resultant, or Newton, matrices is the critical step in the computation of the multivariate resultant and the solution of a nonlinear system. We reveal and exploit the quasi-Toeplitz structure of the Newton matrix, thus decreasing the time complexity of constructing such matrices by roughly one order of magnitude to achieve quasi-quadratic complexity in the matrix dimension. The space complexity is also decreased analogously. These results imply similar improvements in the complexity of computing the resultant polynomial itself and of solving zero-dimensional systems. Our approach relies on fast vector-by-matrix multiplication and uses the following two methods as building blocks. First, a fast and numerically stable method for determining the rank of rectangular matrices, which works exclusively over floating point arithmetic. Second, exact polynomial arithmetic algorithms that improve upon the complexity of polynomial multiplication under our model of sparseness, offering bounds linear in the number of variables and the number of non-zero terms.

A subdivision-based algorithm for the sparse resultant

Multivariate resultants generalize the Sylvester resultant of two polynomials and characterize the solvability of a polynomial system. They also reduce the computation of all common roots to a problem in linear algebra. We propose a determinantal formula for the sparse resultant of an arbitrary system of n + 1 polynomials in n variables. This resultant generalizes the classical one and has significantly lower degree for polynomials that are sparse in the sense that their mixed volume is lower than their Bézout number. Our algorithm uses a mixed polyhedral subdivision of the Minkowski sum of the Newton polytopes in order to construct a Newton matrix. Its determinant is a nonzero multiple of the sparse resultant and the latter equals the GCD of at most n + 1 such determinants. This construction implies a restricted version of an effective sparse Nullstellensatz. For an arbitrary specialization of the coefficients, there are two methods that use one extra variable and yield the sparse resultant. This is the first algorithm to handle the general case with complexity polynomial in the resultant degree and simply exponential in n . We conjecture its extension to producing an exact rational expression for the sparse resultant.

A Nonlinear Finite Element Eigenanalysis of Singular Plane Stress Fields in Bimaterial Wedges Including Complex Eigenvalues

A finite element formulation is developed for the analysis of variable-separable singular stress fields in power law hardening materials under conditions of plane stress. The displacement field within a sectorial element is assumed to be quadratic in the angular coordinate and of the power type in the radial direction as measured from the singular point. An iteration scheme that combines the Newton method and matrix singular value decomposition is used to solve the nonlinear homogeneous eigenvalue problem, where the eigenvalues and eigenfunctions are obtained simultaneously. The formulation and iteration scheme apply when the eigenvalue is complex. The examples considered include the single material crack and wedge to demonstrate convergence, and the bimaterial interface crack and the bimaterial wedge to demonstrate geometric versatility and the ability to handle complex eigenvalues. It is found that the real part of the complex eigenvalue for the interface crack agrees with the HRR value. In this case the associated complex eigenfunction is converted into an approximate real-valued eigenfunction that is valid for any mode-mix. In addition, the behavior of separable solutions near certain 'wedge paradox' geometries where non-separable solutions occur is investigated.

Using Implicit ODE Methods with Iterative Linear Equation Solvers in Spectral Methods

The use of Chebyshev spectral methods in space on an evolution partial differential equation results, in many cases, in a stiff system of ordinary differential equations (ODEs). ODE solvers based on explicit methods will therefore be inefficient for such equations. The use of conventional implicit ODE solvers is difficult since the Jacobian matrix of the ODE system is full and large.It is shown how to apply ODE solvers with iterative linear equation solvers to ODEs coming from spectral discretizations. These ODE solvers do not use the Jacobian explicitly, but good preconditioners for the Newton matrix are needed for their efficient operation. Preconditioners of the tensor product type are developed for certain classes of hyperbolic partial differential equations (PDEs), and numerical experiments show that these preconditioners give a good performance of the ODE solver and are a substantial improvement over the performance of explicit solvers.

An approach to solving multiparticle diffusion exhibiting nonlinear stiff coupling

A methodology for handling a class of stiff multiparticle parabolic PDE's in one and two dimensions is presented. The particular example considered in this work is the interaction and diffusion of two point defects in silicon, interstitials and vacancies. Newton's method, latency techniques, and second-order time-stepping approaches all contribute in significantly reduced computation times. A general class of diffusion-reaction problems is defined and conditions under which the corresponding Newton matrix is invertible and Newton's method converges to a globally unique solution are derived. The convergence properties of the purely reactive system are also derived and compared to those given by a Picard iteration. Application of basic iterative matrix techniques for the general diffusion-reaction system is discussed and specific numerical examples of point defect kinetics are given.

An Approach to Solving Multiparticle Diffusion Exhibiting Nonlinear Stiff Coupling

A methodology for handling a class of stiff multiparticle parabolic PDE's in one and two dimensions is presented. The particular example considered in this work is the interaction and diffusion of two point defects in silicon, interstitials and vacancies. Newton's method, latency techniques, and second-order time-stepping approaches all contribute in significantly reduced computation times. A general class of diffusion-reaction problems is defined and conditions under which the corresponding Newton matrix is invertible and Newton's method converges to a globally unique solution are derived. The convergence properties of the purely reactive system are also derived and compared to those given by a Picard iteration. Application of basic iterative matrix techniques for the general diffusion-reaction system is discussed and specific numerical examples of point defect kinetics are given.