Numerical Gradient Research Articles

Development of the FMO/RI-MP2 Fully Analytic Gradient Using a Hybrid-Distributed/Shared Memory Programming Model.

The fully analytic gradient of the second-order Møller-Plesset perturbation theory (MP2) with the resolution-of-the-identity (RI) approximation in the fragment molecular orbital (FMO) framework is derived and implemented using a hybrid multilevel parallel programming model, a combination of the general distributed data interface (GDDI) and the OpenMP API. The FMO/MP2 analytic gradient contains three parts, i.e., the internal fragment component, the electrostatic potential (ESP) component, and the response terms. The RI approximation is applied to the internal fragment MP2 gradient term, whose MP2 densities and monomer MP2 Lagrangians are shared with the ESP and the response terms. The FMO/RI-MP2 analytic gradient implementation is validated against the numerical gradient (with errors ∼10-6-10-5 Hartree/Bohr) and the energy conservation in molecular dynamics (MD) simulations using NVE ensembles. The RI approximation introduces an error of ∼10-5 Hartree/Bohr with a speedup of 4.0-8.0× compared with the currently available GDDI FMO/MP2 gradient. The node linear scaling of the fragmentation framework due to multilevel parallelism is well-preserved and is demonstrated in single-point gradient calculations of large water clusters (e.g., 1120 and 2165 molecules) using 300-800 KNL compute nodes with a parallel efficiency of more than 90%.

Introduction to Dynamic Stochastic Computing

Stochastic computing (SC) is an old but reviving computing paradigm for its simple data path that can perform various arithmetic operations. It allows for low power implementation, which would otherwise be complex using the conventional positional binary coding. In SC, a number is encoded by a random bit stream of '0's and '1's with an equal weight for every bit. However, a long bit stream is usually required to achieve a high accuracy. This requirement inevitably incurs a long latency and high energy consumption in an SC system. In this article, we present a new type of stochastic computing that uses dynamically variable bit streams, which is, therefore, referred to as dynamic stochastic computing (DSC). In DSC, a random bit is used to encode a single value from a digital signal. A sequence of such random bits is referred to as a dynamic stochastic sequence. Using a stochastic integrator, DSC is well suited for implementing accumulation-based iterative algorithms such as numerical integration and gradient descent. The underlying mathematical models are formulated for functional analysis and error estimation. A DSC system features a higher energy efficiency than conventional computing using a fixed-point representation with a power consumption as low as conventional SC. It is potentially useful in a broad spectrum of applications including signal processing, numerical integration and machine learning.

Expedited antenna optimization with numerical derivatives and gradient change tracking

PurposeThe purpose of this study is to propose a framework for expedited antenna optimization with numerical derivatives involving gradient variation monitoring throughout the optimization run and demonstrate it using a benchmark set of real-world wideband antennas. A comprehensive analysis of the algorithm performance involving multiple starting points is provided. The optimization results are compared with a conventional trust-region (TR) procedure, as well as the state-of-the-art accelerated TR algorithms.Design/methodology/approachThe proposed algorithm is a modification of the TR gradient-based algorithm with numerical derivatives in which a monitoring of changes of the system response gradients is performed throughout the algorithm run. The gradient variations between consecutive iterations are quantified by an appropriately developed metric. Upon detecting stable patterns for particular parameter sensitivities, the costly finite differentiation (FD)-based gradient updates are suppressed; hence, the overall number of full-wave electromagnetic (EM) simulations is significantly reduced. This leads to considerable computational savings without compromising the design quality.FindingsMonitoring of the antenna response sensitivity variations during the optimization process enables to detect the parameters for which updating the gradient information is not necessary at every iteration. When incorporated into the TR gradient-search procedures, the approach permits reduction of the computational cost of the optimization process. The proposed technique is dedicated to expedite direct optimization of antenna structures, but it can also be applied to speed up surrogate-assisted tasks, especially solving sub-problems that involve performing numerous evaluations of coarse-discretization models.Research limitations/implicationsThe introduced methodology opens up new possibilities for future developments of accelerated antenna optimization procedures. In particular, the presented routine can be combined with the previously reported techniques that involve replacing FD with the Broyden formula for directions that are satisfactorily well aligned with the most recent design relocation and/or performing FD in a sparse manner based on relative design relocation (with respect to the current search region) in consecutive algorithm iterations.Originality/valueBenchmarking against a conventional TR procedure, as well as previously reported methods, confirms improved efficiency and reliability of the proposed approach. The applications of the framework include direct EM-driven design closure, along with surrogate-based optimization within variable-fidelity surrogate-assisted procedures. To the best of the authors’ knowledge, no comparable approach to antenna optimization has been reported elsewhere. Particularly, it surmounts established methodology by carrying out constant supervision of the antenna response gradient throughout successive algorithm iterations and using gathered observations to properly guide the optimization routine.

All-in-one silicon photonic polarization processor

Abstract With the great developments in optical communication technology and large-scale optical integration technology, it is imperative to realize the traditional functions of polarization processing on an integration platform. Most of the existing polarization devices, such as polarization multiplexers/demultiplexers, polarization controllers, polarization analyzers, etc., perform only a single function. Definitely, integrating all these polarization functions on a chip will increase function flexibility and integration density and also cut the cost. In this article, we demonstrate an all-in-one chip-scale polarization processor based on a linear optical network. The polarization functions can be configured by tuning the array of phase shifters on the chip. We demonstrate multiple polarization processing functions, including those of a multiple-input-multiple-output polarization descrambler, polarization controller, and polarization analyzer, which are the basic building blocks of polarization processing. More functions can be realized by using an additional two-dimensional output grating. A numerical gradient descent algorithm is employed to self-configure and self-optimize these functions. Our demonstration suggests great potential for chip-scale, reconfigurable, and fully programmable photonic polarization processors with the artificial intelligence algorithm.

Superconvergence of the gradient approximation for weak Galerkin finite element methods on nonuniform rectangular partitions

This article presents a superconvergence for the gradient approximation of the second order elliptic equation discretized by weak Galerkin finite element methods on nonuniform rectangular partitions. The result shows a convergence of O(hr), 1.5≤r≤2, for the numerical gradient obtained from the lowest order weak Galerkin element consisting of piecewise linear and constant functions. For this numerical scheme, the optimal order of error estimate is O(h) for the gradient approximation. The superconvergence reveals a superior performance of the weak Galerkin finite element methods. Some computational results are included to numerically validate the superconvergence theory.

Two-order approximate and large stepsize numerical direction based on the quadratic hypothesis and fitting method

Many effective optimization algorithms require partial derivatives of objective functions, while some optimization problems ʼ objective functions have no derivatives. According to former research studies, some search directions are obtained using the quadratic hypothesis of objective functions. Based on derivatives, quadratic function assumptions, and directional derivatives, the computational formulas of numerical first-order partial derivatives, second-order partial derivatives, and numerical second-order mixed partial derivatives were constructed. Based on the coordinate transformation relation, a set of orthogonal vectors in the fixed coordinate system was established according to the optimization direction. A numerical algorithm was proposed, taking the second order approximation direction as an example. A large stepsize numerical algorithm based on coordinate transformation was proposed. Several algorithms were validated by an unconstrained optimization of the two-dimensional Rosenbrock objective function. The numerical second order approximation direction with the numerical mixed partial derivatives showed good results. Its calculated amount is 0.2843 % of that of without second-order mixed partial derivative. In the process of rotating the local coordinate system 360 ° , because the objective function is more complex than the quadratic function, if the numerical direction derivative is used instead of the analytic partial derivative, the optimization direction varies with a range of 103.05 ° . Because theoretical error is in the numerical negative gradient direction, the calculation with the coordinate transformation is 94.71 % less than the calculation without coordinate transformation. If there is no theoretical error in the numerical negative gradient direction or in the large-stepsize numerical optimization algorithm based on the coordinate transformation, the sawtooth phenomenon occurs. When each numerical mixed partial derivative takes more than one point, the optimization results cannot be improved. The numerical direction based on the quadratic hypothesis only requires the objective function to be obtained, but does not require derivability and does not take into account truncation error and rounding error. Thus, the application scopes of many optimization methods are extended.

Wavefront propagation based on the ray transfer matrix and numerical orthogonal Zernike gradient polynomials.

The aberrated wavefront propagates along its normal. Both the magnitude and boundary change after the propagation. Wavefronts characterized by Zernike coefficients and a normalized pupil radius can also be represented by a bundle of feature rays normal to the local surface. A ray transfer matrix parameterized by the pupil radius and propagation distance is proposed to transfer these feature rays to obtain the slope and position data of the propagated feature rays. Numerical orthogonal Zernike gradient polynomials are derived to reconstruct the wavefront from the discrete data by using a numerical method. Two aberrated wavefronts are performed as examples to validate the accuracy and flexibility of the proposed numerical method.

Shape optimisation of blended-wing-body underwater gliders based on free-form deformation

ABSTRACTAs a type of Autonomous Underwater Vehicles (AUV), blended-wing-body underwater gliders (BWBUGs) have drawn much attention due to their unique performance advantages in long-distance navigation missions. Because the configuration is the direct component to influence the navigation efficiency, it is significant to perform shape optimisation of BWBUGs. In this paper, an optimisation framework is presented, where Free-Form Deformation (FFD) is adopted for geometric parameterisation. Besides, an Euler-based CFD solver with a discrete adjoint method is applied for numerical simulation and gradients calculation, and Sequential Quadratic Programming (SQP) is employed as the optimiser. With the help of the proposed framework, an optimised BWBUG is regarded as the initial shape and four shape optimisation cases are carried out for different design purposes. All the objective functions aim to increase the lift-drag ratio under the thickness and volume constraints. The simulation results show that all the cases achieve a higher lift-drag ratio.

A full additive QM/MM scheme for the computation of molecular crystals with extension to many-body expansions.

An additive quantum mechanics/molecular mechanics (QM/MM) model for the theoretical investigation of molecular crystals (AC-QM/MM) is presented. At the one-body level, a single molecule is chosen as the QM region. The MM region around it consists of a finite cluster of explicit MM atoms, represented by point charges and Lennard-Jones potentials, with additional background charges to mimic periodic electrostatics. Cluster charges are QM-derived and calculated self-consistently to ensure a polarizable embedding. We have also considered the extension to many-body QM corrections, calculating the interactions of a central molecule to neighboring units in the crystal. Full gradient expressions have been derived, also including symmetry information. The scheme allows for the calculation of molecular properties as well as unconstrained optimizations of the molecular geometry and cell parameters with respect to the lattice energy. Benchmarking the approach with the X23 reference set confirms the convergence pattern of the many-body extension although a comparison to plane-wave density functional theory reveals a systematic overestimation of cohesive energies by 6-16 kJ mol-1. While the scheme primarily aims to provide an inexpensive and flexible way to model a molecule in a crystal environment, it can also be used to reach highly accurate cohesive energies by the straightforward application of wave function correlated approaches. Calculations with local coupled cluster with singles, doubles, and perturbative triples, albeit limited to numerical gradients, show an impressive agreement with experimental estimates for small molecular crystals.

Superconvergence of numerical gradient for weak Galerkin finite element methods on nonuniform Cartesian partitions in three dimensions

A superconvergence error estimate for the gradient approximation of the second order elliptic problem in three dimensions is analyzed by using weak Galerkin finite element scheme on the uniform and non-uniform cubic partitions. Due to the loss of the symmetric property from two dimensions to three dimensions, this superconvergence result in three dimensions is not a trivial extension of the recent superconvergence result in two dimensions Li et al. (0000) from rectangular partitions to cubic partitions. The error estimate for the numerical gradient in the L2-norm arrives at a superconvergence order of O(hr)(1.5≤r≤2) when the lowest order weak Galerkin finite elements consisting of piecewise linear polynomials in the interior of the elements and piecewise constants on the faces of the elements are employed. A series of numerical experiments are illustrated to confirm the established superconvergence theory in three dimensions.

Supercloseness of Linear DG-FEM and Its Superconvergence Based on the Polynomial Preserving Recovery for Helmholtz Equation

In this paper we study the supercloseness property of the linear discontinuous Galerkin (DG) finite element method and its superconvergence behavior after post-processing by the polynomial preserving recovery (PPR). The error estimate with explicit dependence on the wave number k, the penalty parameter $$\mu $$ and the mesh condition parameter $$\alpha $$ is derived. We prove the supercloseness between the DG finite element solution and the linear interpolation and the superconvergence for the recovered gradient by the PPR under the assumption $$k(kh)^2\le C_0$$ (h is the mesh size) and certain mesh conditions. Furthermore, we estimate the error between the DG numerical gradient and recovered gradient, which motivates us to define the a posteriori error estimator and design a Richardson extrapolation to post-process the recovered gradient by PPR. Finally, some numerical examples are provided to confirm the theoretical results of superconvergence analysis.

Numerical Gradient Schemes for Heat Equations Based on the Collocation Polynomial and Hermite Interpolation

As is well-known, the advantage of the high-order compact difference scheme (H-OCD) is that it is unconditionally stable and convergent on the order O ( τ 2 + h 4 ) (where τ is the time step size and h is the mesh size), under the maximum norm for a class of nonlinear delay partial differential equations with initial and Dirichlet boundary conditions. In this article, a new numerical gradient scheme based on the collocation polynomial and Hermite interpolation is presented. The convergence order of this kind of method is also O ( τ 2 + h 4 ) under the discrete maximum norm when the spatial step size is twice the one of H-OCD, which accelerates the computational process. In addition, some corresponding analyses are made and the Richardson extrapolation technique is also considered in the time direction. The results of numerical experiments are consistent with the theoretical analysis.

Convergence of the Numerical Scheme for Regularised Riemannian Mean Curvature Flow Equation

Abstract The aim of the paper is to study problem of image segmentation and missing boundaries completion introduced in [Mikula, K.—Sarti, A.––Sgallarri, A.: Co-volume method for Riemannian mean curvature flow in subjective surfaces multiscale segmentation, Comput. Vis. Sci. 9 (2006), 23–31], [Mikula, K.—Sarti, A.—Sgallari, F.: Co-volume level set method in subjective surface based medical image segmentation, in: Handbook of Medical Image Analysis: Segmentation and Registration Models (J. Suri et al., eds.), Springer, New York, 583–626, 2005], [Mikula, K.—Ramarosy, N.: Semi-implicit finite volume scheme for solving nonlinear diffusion equations in image processing, Numer. Math. 89 (2001), 561–590] and [Tibenský, M.: VyužitieMetód Založených na Level Set Rovnici v Spracovaní Obrazu, Faculty of mathematics, physics and informatics, Comenius University, Bratislava, 2016]. We generalize approach presented in [Eymard, R.—Handlovičová, A.—Mikula, K.: Study of a finite volume scheme for regularised mean curvature flow level set equation, IMA J. Numer. Anal. 31 (2011), 813–846] and apply it in the field of image segmentation. The so called regularised Riemannian mean curvature flow equation is presented and the construction of the numerical scheme based on the finite volume method approach is explained. The principle of the level set, for the first time given in [Osher, S.—Sethian, J. A.: Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988), 12–49] is used. Based on the ideas from [Eymard, R.—Handlovičová, A.– –Mikula, K.: Study of a finite volume scheme for regularised mean curvature flow level set equation, IMA J. Numer. Anal. 31 (2011), 813–846] we prove the stability estimates on the numerical solution and the uniqueness of the numerical solution. In the last section, there is a proof of the convergence of the numerical scheme to the weak solution of the regularised Riemannian mean curvature flow equation and the proof of the convergence of the approximation of the numerical gradient is mentioned as well.

The dynamics of open populations: integration of top–down, bottom–up and supply–side influences on intertidal oysters

Most communities are structured not by a single process but by some combination of top–down, bottom–up and supply–side (i.e. juvenile recruitment) factors. However, establishing how multiple processes interact remains a fundamental challenge. For example, the recruitment, growth, and mortality of estuarine species can vary along the steep and numerous environmental and biological gradients typical of these habitats, but the relative importance of those gradients is generally unknown. We took a novel approach to this question by coupling long‐term field observations of the Olympia oyster Ostrea lurida in a central California estuary with a state–space integral projection model. This approach revealed that the most parsimonious description of oyster population dynamics involved spatial variation in growth and adult mortality – but not juvenile mortality – as well as spatiotemporal variation in recruitment. These patterns match the available short‐term estimates of each of those processes from field studies, and reveal a synthetic view of oyster population dynamics. Larval recruitment has an interannual ‘boom and bust’ pattern, and during good recruitment years most larvae settle in the inner bay where water residence time is highest. Adult oyster mortality is also highest in the inner bay, where several invasive predators are abundant and lowest in the mid‐bay, where oyster growth is greatest (due to bottom–up factors), likely leading to a size refuge from native predators. Surprisingly, juvenile mortality was constant across the bay, possibly because of a lack of size refuge from native and invasive predators. Our research approach represents an important advance in disentangling the contributions of spatio–temporal variation in top–down, bottom–up and supply–side forces to the dynamics of populations with open recruitment.

Optimization of fixed-order controllers using exact gradients

Finding good controller settings that satisfy complex design criteria is not trivial. This is also the case for simple fixed-order controllers including the three parameter pid controller. To be rigorous, we formulate the design problem into an optimization problem. However, the algorithm may fail to converge to the optimal solution because of inaccuracies in the estimation of the gradients needed for this optimization. In this paper we derive exact gradients for the problem of optimizing performance (iae for input and output disturbances) with constraints on robustness (MS, MT). Exact gradients give increased accuracy and better convergence properties than numerical gradients, including forward finite differences. The approach may be easily extended to other control objectives and other fixed-order controllers, including Smith Predictor and pidf controllers.

Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz

The variational quantum eigensolver (VQE) algorithm combines the ability of quantum computers to efficiently compute expectation values with a classical optimization routine in order to approximate ground state energies of quantum systems. In this paper, we study the application of VQE to the simulation of molecular energies using the unitary coupled cluster (UCC) ansatz. We introduce new strategies to reduce the circuit depth for the implementation of UCC and improve the optimization of the wavefunction based on efficient classical approximations of the cluster amplitudes. Additionally, we propose an analytical method to compute the energy gradient that reduces the sampling cost for gradient estimation by several orders of magnitude compared to numerical gradients. We illustrate our methodology with numerical simulations for a system of four hydrogen atoms that exhibit strong correlation and show that the circuit depth of VQE using a UCC ansatz can be reduced without introducing significant loss of accuracy in the final wavefunctions and energies.

Analysis of solute-solvent interactions using the solvation model density combined with the fragment molecular orbital method

The energy and its analytic gradient are developed for the solvation model density (SMD) combined with the fragment molecular orbital (FMO) method. The accuracy of the energy is evaluated in comparison to full results without fragmentation for a set of neutral and charged polypeptides. The accuracy of the gradient is computed in comparison to numerical gradient. The components to the solvation energy in SMD are compared to the polarizable continuum model (PCM). FMO with SMD and PCM is applied to analyze protein–ligand binding in solution.

Symmetric critical knots for O'Hara's energies

We prove the existence of symmetric critical torus knots for O'Hara's knot energy family Eα, α∈(2,3) using Palais' classic principle of symmetric criticality. It turns out that in every torus knot class there are at least two smooth Eα-critical knots, which supports experimental observations using numerical gradient flows.

Validation of the Octavius 4D measuring system in verifying advanced external beams radiotherapy techniques

The development of advanced techniques as Intensity modulated radiotherapy (IMRT) and volumetric modulated arc therapy (VMAT) are characterized by their numerous steep dose gradient using the dynamic MLC modality and these requires specific quality assurance (QA) to meet this complexity. The conventional dosimetery techniques, such as ionization chamber, point dose measurements and film dosimetry, are gradually being replaced by detector arrays which can produce accurate results immediately. There are many commercial detector arrays one of them is 2D-ARRAY seven29 (PTW, Freiburg, Germany). We used the Octavius 4 D measuring system with its detector array, ionization chamber. The accuracy of the synchronizing motion of the phantom was validated. It was found that the detector array response is stable after warming up with 800 MU, the dose rate measurements and the dose response is linear with R2=0.929. The gamma index for the clinical IMRT plans were found to be 93.3%, 98.0 %, 95.5% and 96.2% for head and neck, prostate, breast and brain cases, and for VMAT plans are 96.3%, 99.4 %, 98.1% and 98.6% for head and neck, prostate, breast and brain, respectively, using gamma criteria 3%/3 mm. The aim of this work is to characterize and validate the Octavius 4D measuring system for 6MV beams. Also, to verify clinical IMRT and VMAT plans, the volumetric dose matrix of the clinical plans, dose line profiles and to study the factors affect the gamma index.

Quasi-Newton algorithm for optimal approximate linear regression design: Optimization in matrix space

Given a linear regression model and an experimental region for the independent variable, the problem of finding an optimal approximate design calls for minimizing a convex optimality criterion over a convex set of information matrices of feasible approximate designs. For numerical solution pure gradient methods are often used by design theorists, as vertex direction, vertex exchange, multiplicative algorithms, or combinations hereof. These methods have two major deficiencies: a slow convergence rate after a quick but rough approximation to the optimum, and often a large support of the obtained nearly optimal design. The latter feature is related to the fact that the methods optimize in the space of design measures which is usually of high or even infinite dimension, whereas the dimension of the information matrices is often small or moderate. For such situations a quasi-Newton method is revisited which was originally established by Gaffke & Heiligers (1996). In the present paper new possibilities of its application are demonstrated. The algorithm optimizes in matrix space. It shows a good global and an excellent local convergence behavior resulting in an accurate approximation of the optimum. A crucial subroutine solves convex quadratic minimization over the set of information matrices via repeated linear minimization over that set, providing thus the quasi-Newton step of the algorithm. This may also be of interest as a tool for computing an approximate design from a given information matrix and such that the support size of the design keeps Carathéodory’s bound. Illustrations are given for D- and I-optimality in particular multivariate random coefficient regression models and for T-optimal discriminating design in univariate polynomial models. Moreover, the behavior of the algorithm is tested for cases of larger dimensions: D- and I-optimal design for a third order polynomial model in several variables on a discretized cube.