Numerical Experiments Research Articles

Optimal box contraction for solving linear systems via simulated and quantum annealing

Solving linear systems of equations is an important problem in engineering. Many quantum algorithms, such as the Harrow–Hassidim–Lloyd algorithm and the box algorithm, have been proposed for solving such systems. The focus of this article is on improving the efficiency of the box algorithm. The basic principle behind this algorithm is to transform the linear system into a series of quadratic unconstrained binary optimization (QUBO) problems, which are then solved on annealing machines. The computational efficiency of the box algorithm is entirely determined by the number of iterations, which, in turn, depends on the box contraction ratio, typically set to 0.5. Here, it is shown through theoretical analysis that a contraction ratio of 0.5 is sub-optimal and that a computational speed-up can be achieved with a contraction ratio of 0.2. This is confirmed through numerical experiments where a computational speed-up between 20 % to 60 % is observed when the optimal contraction ratio is used.

A trust region algorithm for finite-strain plasticity with strongly coupled hardening

As extensions to our Lagrangian finite-strain plasticity framework based on the approximate exponential integrator, we introduce two new algorithms: (i) a fixed-radius trust region root finder for the nonlinear system involving the elastic right Cauchy Green tensor and plastic multiplier (ii) a partitioned approach for the hardening variables, which are determined in a staggered form using the strongly-coupled concept typically adopted for multiphysics problems. This allows the use of intricate hyperelastic laws combined with recent yield functions, which otherwise would involve a laborious treatment, and the use of corresponding work-hardening. Work-hardening would introduce significant nonlinearities in the constitutive system if used in a fully-coupled form. For the partitioned approach, a dynamic relaxation algorithm is adopted. This allows the efficient solution of the two nonlinear equations without significant drifting. Results show that robustness and efficiency are significantly improved. Herein, algorithms are described in detail. Significant testing is performed with imposed strains up to 100 for a carbon steel. This contributes to the robustness of equilibrium iterations. Drifting is also assessed as a function of number of steps in the dynamic relaxation algorithm. Numerical experimentation is performed in the 3D tension test for an anisotropic yield function.

Assessing Different Two-Stage Stochastic Models for Optimizing Food Bank Networks’ Operations During Natural Disasters

Humanitarian logistics face significant challenges during natural disasters due to operational uncertainties. Humanitarian logistics networks such as food banks must manage both regular operations and disaster-induced supply and demand. The study aims to develop and assess two-stage stochastic models that support decision-making under these dual operations. We evaluate various decisional strategies through extensive numerical experiments inspired in the operation of the food bank network Bancos de Alimentos de México (BAMX), highlighting the importance of suitable timeframes for reactive and anticipative decisions. The findings offer valuable insights for managers in balancing routine and emergency responses efficiently.

Streaming First‐Order Optimization Methods With Momentum for Dominant Eigenspace Estimation

ABSTRACTIn this article, we consider the dominant eigenspace recovery for PCA in the streaming scenario, and design several new algorithms for this problem by incorporating various momentum schemes into two classical gradient based streaming eigensolvers, which emerges as a frequently‐used trick in improving gradient type methods. Theoretically, we establish the convergence guarantees of the proposed first‐order momentum methods, and prove their decreasing variance under favorable conditions. Extensive numerical experiments are conducted to corroborate the result that momentum does yield variance reduction, and also highlight the advantages of lowering the stepsize sensitivity. In addition, our methods with appropriately chosen momentum yield a dramatically enhanced performance in terms of both the convergence speed and accuracy.

Water budget of the Caspian Sea by numerical experiments with ocean circulation model INMIO-CICE in the last glacial maximum and pre-industrial period

We used the hydrodynamic model of the Caspian Sea, INMIO-CICE, to calculate equilibrium river runoff and evaporation from the sea surface over a wide range of sea levels (from –85 to +50 m asl) for different climatic conditions: the Last Glacial Maximum (about 21kyr) and pre-industrial climate (~1850 CE). Data from the climate model INMCM4.8 were used as boundary conditions. It was found that to maintain sea level at 35–50 m asl, corresponding to the maximum values of the Khvalynian transgression, a river runoff of about 400 km3/year was required in the Last Glacial Maximum. In the Last Glacial Maximum evaporation from the sea surface decreased by 105–170 mm (12–22%), and precipitation, according to the INMCM4.8 model, by 50–70 mm (15–30%). This caused the equilibrium runoff to decrease by about 10–20% compared to pre-industrial conditions. Smaller absolute and relative changes correspond to lower sea levels. The maximum decrease in evaporation occurred at 5 m asl.

Enhancing Accuracy and Efficiency in Stiff ODE Integration Using Variable Step Diagonal BBDF Approaches

Recent advancements in mathematical modelling have uncovered a growing number of systems exhibiting stiffness, a phenomenon that challenges the effectiveness of traditional numerical methods. Motivated by the need for more robust numerical techniques to address this issue, this paper presents an enhanced version of the Diagonally Block Backward Differentiation Formula (BBDF) that incorporates intermediate points, known as off-step points, to improve the accuracy and efficiency of solutions for stiff ordinary differential equations (ODEs). The new scheme leverages an adaptive step-size strategy to refine accuracy and efficiency between regular and off-grid integration steps. Theoretical analysis confirms that the proposed scheme is an A-stable and convergent method, as it satisfies the fundamental criteria of consistency, zero-stability, and A-stability. Numerical experiments on single and multivariable systems across varying time scales demonstrate significant improvements in solving stiff ODEs compared to existing techniques. Therefore, the new proposed method is an effective solver for stiff ODEs.

Noise-resilient designs and analysis for optical neural networks

Abstract All analog signal processing is fundamentally subject to noise, and this is also the case in next generation implementations of Optical Neural Networks (ONNs). Therefore, we propose the first hardware-based approach to mitigate noise in ONNs. A tree-like and an accordion-like design are constructed from a given Neural Network (NN) that one wishes to implement. Both designs have the capability that the resulting ONNs give outputs close to the desired solution. To establish the latter, we analyze the designs mathematically. Specifically, we investigate a probabilistic framework for the tree-like design that establishes the correctness of the design, i.e., for any feed-forward NN with Lipschitz continuous activation functions, an ONN can be constructed that produces output arbitrarily close to the original. ONNs constructed with the tree-like design thus also inherit the universal approximation property of NNs. For the accordion-like design, we restrict the analysis to NNs with linear activation functions and characterize the ONNs’ output distribution using exact formulas. Finally, we report on numerical experiments with LeNet ONNs that give insight into the number of components required in these designs for certain accuracy gains. The results indicate that adding just a few components and/or adding them only in the first (few) layers in the manner of either design can already be expected to increase the accuracy of ONNs considerably. To illustrate the effect we point to a specific simulation of a LeNet implementation, in which adding one copy of the layers components in each layer reduces the Mean Squared Error (MSE) by 59.1% for the tree-like design and by 51.5% for the accordion-like design. In this scenario, the gap in accuracy of prediction between the noiseless NN and the ONNs reduces even more: 93.3% for the tree-like design and 80% for the accordion-like design.

Time-lapse difference seismic inversion based on scattering theory

Time-lapse seismic monitoring technology is an effective tool for quantitatively estimating subsurface medium variations and plays a crucial role in resource and engineering explorations. Among existing time-lapse inversion methods, the difference inversion technique offers higher accuracy but requires linear inversion operators. During linearization, the background is often approximated as a homogeneous medium, resulting in the loss of effective information. Additionally, merely using reflection wave information is not conducive to improving inversion accuracy. However, the true structure of the subsurface is often complex and inhomogeneous, which can easily generate complex wavefields, such as scattering waves, thereby hindering the resolution of differences between time-lapse models. Therefore, incorporating more comprehensive wavefield information is essential for enhancing the accuracy of time-lapse seismic inversion. To address this issue, a wave-equation-based time-lapse difference inversion method is proposed to retrieve the production-induced property perturbation. In the proposed method, the forward operator is formulated through the scattering wave equation, the operator is linearized through the Born approximation, and the WKBJ method is used to approximate Green's function with high accuracy for the inhomogeneous background. In this study, time-lapse seismic parallel inversion and difference inversion are conducted using the proposed forward operator. Both numerical experiments and application of field seismic data indicate that the proposed time-lapse difference inversion method can realize a more accurate quantitative characterization of time-lapse variations.

Wasserstein task embedding for measuring task similarities

Measuring similarities between different tasks is critical in a broad spectrum of machine learning problems, including transfer, multi-task, continual, and meta-learning. Most current approaches to measuring task similarities are architecture-dependent: (1) relying on pre-trained models, or (2) training networks on tasks and using forward transfer as a proxy for task similarity. In this paper, we leverage the optimal transport theory and define a novel task embedding for supervised classification that is model-agnostic, training-free, and capable of handling (partially) disjoint label sets. In short, given a dataset with ground-truth labels, we perform a label embedding through multi-dimensional scaling and concatenate dataset samples with their corresponding label embeddings. Then, we define the distance between two datasets as the 2-Wasserstein distance between their updated samples. Lastly, we leverage the 2-Wasserstein embedding framework to embed tasks into a vector space in which the Euclidean distance between the embedded points approximates the proposed 2-Wasserstein distance between tasks. We show that the proposed embedding leads to a significantly faster comparison of tasks compared to related approaches like the Optimal Transport Dataset Distance (OTDD). Furthermore, we demonstrate the effectiveness of our embedding through various numerical experiments and show statistically significant correlations between our proposed distance and the forward and backward transfer among tasks on a wide variety of image recognition datasets.

A spatial Durbin model for interval-valued data with t-distribution

Interval-valued data have received much attention due to their wide applications in finance, econometrics, meteorology and medicine. However, most regression models developed for interval-valued data assume that the observations are mutually independent, which can not be adapted to scenarios where individuals are spatially correlated. We propose a new linear model to account for areal-type spatial dependencies present in interval-valued data. Specifically, two spatial Durbin models are separately built for the centers and log-ranges of the interval-valued response. To make the new model more robust to heavy-tailed noise, we assume the error term follows a t-distribution. The parameters are obtained using an expectation-maximization (EM) algorithm. Numerical experiments are designed to examine the performance of the proposed method. We also use a weather dataset to demonstrate the usefulness of our model.

RANS Simulation of Minimum Ignition Energy of Stoichiometric and Leaner CH4/Air Mixtures at Higher Pressures in Quiescent Conditions

Minimum ignition energy (MIE) has been extensively studied via experiments and simulations. However, our literature review reveals little quantitative consistency, with results varying from 0.324 to 1.349 mJ for ϕ = 1.0 and from 0.22 to 0.944 mJ for ϕ = 0.9. Therefore, there is a need to resolve these discrepancies. This RANS study aims to partially address this knowledge gap. Additionally, it presents other flame evolution parameters essential for robust combustion design. Using the reactingFOAM solver, we predict the threshold energy required to ignite the fuel mixture. For this, the single step using the Arrhenius law is selected to model ignition in the flame kernel of stochiometric and lean CH4/air mixtures, allowing it to develop into a self-sustained flame. The ignition power density, an energy quantity normalised with volume, is incrementally varied, keeping the kernel critical radius rs constant at 0.5 mm in the quiescent mixture of two equivalence ratios ϕ 0.9 and 1.0, for varied operating pressures of 1, 5, and 10 bar at the constant initial temperature of 300 K. The minimum ignition energy is validated with twelve independent 1-bar datasets both numerically and experimentally. The effect of pressure on MIEs, which diminish as pressure rises, is significant. At ϕ = 1.0 (and 0.9), the flame temperature reached 481.24 K (457.803 K) at 1 bar, 443.176 K (427.356 K) at 5 bar, and 385.56 K (382.688 K) at 10 bar. The minimum ignition energy was validated using twelve independent 1-bar datasets from both numerical simulations and experiments. The results show strong agreement with many experimental findings. Finally, a mathematical formulation of MIE is devised; a function of pressure and equivalence ratio shows a slightly curved relationship.

Hidden-Markov models for ordinal time series

AbstractA common approach for modeling categorical time series is Hidden-Markov models (HMMs), where the actual observations are assumed to depend on hidden states in their behavior and transitions. Such categorical HMMs are even applicable to nominal data but suffer from a large number of model parameters. In the ordinal case, however, the natural order among the categorical outcomes offers the potential to reduce the number of parameters while improving their interpretability at the same time. The class of ordinal HMMs proposed in this article link a latent-variable approach with categorical HMMs. They are characterized by parametric parsimony and allow the easy calculation of relevant stochastic properties, such as marginal and bivariate probabilities. These points are illustrated by numerical examples and simulation experiments, where the performance of maximum likelihood estimation is analyzed in finite samples. The developed methodology is applied to real-world data from a health application.

Are information criteria good enough to choose the right the number of regimes in hidden Markov models?

Selecting the number of regimes in Hidden Markov models is an important problem. Many criteria are used to do so, such as Akaike information criterion (AIC), Bayesian information criterion (BIC), integrated completed likelihood (ICL), deviance information criterion (DIC), and Watanabe-Akaike information criterion (WAIC), to name a few. In this article, we introduced goodness-of-fit tests for general Hidden Markov models with covariates, where the distribution of the observations is arbitrary, i.e. continuous, discrete, or a mixture of both. A selection procedure is proposed based on this goodness-of-fit test. The main aim of this article is to compare the classical information criterion with the new criterion when the outcome is either continuous, discrete or zero-inflated. Numerical experiments assess the finite sample performance of the goodness-of-fit tests, and comparisons between the different criteria are made.

Impacts of Local and Remote SST Warming on Summer Circulation Changes in the Western North Pacific

Abstract This study explores how future SST warming in remote ocean basins may affect the western North Pacific (WNP) wet season climate by applying a high-resolution atmospheric general circulation model to conduct a series of numerical experiments. A marked precipitation and tropical cyclone (TC) activity reduction, as well as enhanced anticyclonic circulation, in the WNP is projected in AMIP experiments forced by SST change in a future warming scenario. The sensitivity experiments reveal that various SST warming phenomena (e.g., in the global SST warming pattern, the tropical ocean belt, the Indian Ocean, the tropical Atlantic, and the subtropical northeast Pacific) and the increase in greenhouse gas concentration could weaken the precipitation, TC activity, and circulation. By contrast, the SST warming in the WNP and eastern equatorial Pacific has opposite and mixed effects, respectively, and tends to weakly offset the dominant influences of remote ocean warming. These results indicate that the WNP, being the epicenter of the global teleconnection of divergent and rotational flow, is susceptible to the influence of the SST warming in remote ocean basins. The remote forcing as projected in future scenarios would overwhelm the enhancing effect of local SST warming and weaken the circulation, convection, and TC activity in the WNP. These findings further the understanding of how the decreased precipitation and enhanced subtropical high in the WNP may be easily triggered by remote SST warming as revealed in the AMIP-type simulations. How this effect would be affected by air–sea coupling needs further investigation.

Robust and accurate numerical framework for multi-dimensional fractional-order telegraph equations using Jacobi/Jacobi-Romanovski spectral technique

This paper presents a novel spectral algorithm for the numerical solution of multi-dimensional fractional-order telegraph equations, a critical model used to capture the combined effects of diffusion and wave propagation. The core innovation of this work is the application of Jacobi-Romanovski polynomials as the basis functions for spectral discretization. These polynomials offer unique advantages, including the ability to handle nonstandard domains and boundary conditions, making them particularly suitable for partial differential equation (PDE) applications. A comprehensive error analysis is conducted, providing deep insights into the convergence rates and factors affecting the accuracy of the numerical solutions. Extensive numerical experiments further demonstrate the superior performance of the proposed spectral algorithm in solving a wide range of multi-dimensional fractional-order telegraph equation models. The results show a significant improvement in accuracy and computational efficiency compared to traditional numerical methods, such as finite difference or finite element techniques. This research advances the field of computational science by offering a robust, efficient, and versatile numerical framework for the precise solution of complex multi-dimensional PDEs.

Heavy-Tail Phenomenon in Decentralized SGD

Recent theoretical studies have shown that heavy-tails can emerge in stochastic optimization due to ‘multiplicative noise’, even under surprisingly simple settings, such as linear regression with Gaussian data. While these studies have uncovered several interesting phenomena, they consider conventional stochastic optimization problems, which exclude decentralized settings that naturally arise in modern machine learning applications. In this paper, we study the emergence of heavy-tails in decentralized stochastic gradient descent (DE-SGD), and investigate the effect of decentralization on the tail behavior. We first show that, when the loss function at each computational node is twice continuously differentiable and strongly convex outside a compact region, the law of the DE-SGD iterates converges to a distribution with polynomially decaying (heavy) tails. To have a more explicit control on the tail exponent, we then consider the case where the loss at each node is a quadratic function, and show that the tail-index can be estimated as a function of the step-size, batch-size, and the topological properties of the network of the computational nodes. Then, we provide theoretical and empirical results showing that DE-SGD has heavier tails than centralized SGD. We also compare DE-SGD to disconnected SGD where nodes distribute the data but do not communicate. Our theory uncovers an interesting interplay between the tails and the network structure: we identify two regimes of parameters (stepsize and network size), where DE-SGD can have lighter or heavier tails than disconnected SGD depending on the regime. Finally, to support our theoretical results, we provide numerical experiments conducted on linear regression and neural networks that suggest the heavier tails is correlated with better generalization in the decentralized setting.

The Reduced-Dimension Method for Crank–Nicolson Mixed Finite Element Solution Coefficient Vectors of the Extended Fisher–Kolmogorov Equation

A reduced-dimension (RD) method based on the proper orthogonal decomposition (POD) technology and the linearized Crank–Nicolson mixed finite element (CNMFE) scheme for solving the 2D nonlinear extended Fisher–Kolmogorov (EFK) equation is proposed. The method reduces CPU runtime and error accumulation by reducing the dimension of the unknown CNMFE solution coefficient vectors. For this purpose, the CNMFE scheme of the above EFK equation is established, and the uniqueness, stability and convergence of the CNMFE solutions are discussed. Subsequently, the matrix-based RDCNMFE scheme is derived by applying the POD method. Furthermore, the uniqueness, stability and error estimates of the linearized RDCNMFE solution are proved. Finally, numerical experiments are carried out to validate the theoretical findings. In addition, we contrast the RDCNMFE method with the CNMFE method, highlighting the advantages of the dimensionality reduction method.

Conformal symplectic and constraint-preserving model order reduction of constrained conformal Hamiltonian systems

Dissipative and constrained dynamical systems model various physical phenomena comprising damping and constraining forces. These forces give rise to specific geometric properties of the dynamical system's governing equations. In this paper, we reduce and solve high-dimensional parameterized constrained and damped nonlinear differential equations with structure-preserving algorithms. These algorithms preserve the constrained equation's geometric properties, such as conformal symplecticness and phase-space, through model order reduction and numerical integration. We develop a reduced model based on an orthosymplectic reduced basis and establish the well-posedness and conformal symplecticness of the reduced model. Moreover, we derive a priori error bounds for state space and constraints errors. Subsequently, we utilize two RATTLE-based integrators to simulate the full and reduced models while preserving their geometric properties. Numerical experiments on constrained damped oscillators and a constrained lattice grid justify theoretical results and demonstrate the advantages of the structure-preserving computational methods for constrained dynamical systems.

Managing Resources for Shared Micromobility: Approximate Optimality in Large-Scale Systems

We consider the problem of managing resources in shared micromobility systems (bike sharing and scooter sharing). An important task in managing such systems is periodic repositioning/recharging/sourcing of units to avoid stockouts or excess inventory at nodes with unbalanced flows. We consider a discrete-time model; each period begins with an initial inventory at each node in the network, and then, customers (demand) materialize at the nodes. Each customer picks up a unit at the origin node and drops it off at a randomly sampled destination node with an origin-specific probability distribution. We model the above network inventory management problem as an infinite horizon discrete-time discounted Markov decision process (MDP) and prove the asymptotic optimality of a novel mean-field approximation to the original MDP as the number of stations becomes large. To compute an approximately optimal policy for the mean-field dynamics, we provide an algorithm with a running time that is logarithmic in the desired optimality gap. Lastly, we compare the performance of our mean field-based policy with state-of-the-art heuristics via numerical experiments, including experiments using Austin scooter-sharing data. This paper was accepted by Jeannette Song, operations management. Supplemental Material: The online appendix and data files are available at https://doi.org/10.1287/mnsc.2022.02023 .

Strength of viscous subduction interfaces: A global compilation

The long- and short-term dynamics of subduction zones are strongly affected by interface rheology. Current estimates of interface viscosity are based on endmember flow laws or simple mixing models, but the variation of possible materials on the subduction interface theoretically permits viscosity variations of up to five orders of magnitude. To better constrain the range of strength along deep/viscous subduction interfaces, we compiled a global database of rock types and block-and-matrix distributions in exhumed deep subduction interface mélange shear zones. We applied numerical shear experiments to each shear zone to quantify their bulk shear strengths/viscosities. Our results show that natural subduction shear zones have viscosities ranging from ∼1018 to 1020 Pa⋅s and shear strengths of ∼1−100 MPa. The variation among them is dominantly controlled by temperature, but even for shear zones deformed at the same temperature, variations of up to a factor of 50 are observed. These variations are attributed to differences in matrix composition, shear zone width, and block distributions. Our results are consistent with paleopiezometric and force-balance−derived estimates of interface strength, and they confirm the importance of variation in geological inputs and their spatial distribution along the interface for subduction dynamics.