Search

Let $\mathbf S \in \mathbb R^{n \times n}$ satisfy $\|\mathbf 1-\mathbf S\|_2\le\epsilon n$, where $\mathbf 1$ is the all ones matrix and $\|\cdot\|_2$ is the spectral norm. It is well-known that there exists such an $\mathbf S$ with just $O(n/\epsilon^2)$ non-zero entries: we can let $\mathbf S$ be the scaled adjacency matrix of a Ramanujan expander graph. We show that such an $\mathbf S$ yields a $universal$ $sparsifier$ for any positive semidefinite (PSD) matrix. In particular, for any PSD $\mathbf A \in \mathbb{R}^{n\times n}$ with entries bounded in magnitude by $1$, $\|\mathbf A - \mathbf A\circ\mathbf S\|_2 \le \epsilon n$, where $\circ$ denotes the entrywise (Hadamard) product. Our techniques also give universal sparsifiers for non-PSD matrices. In this case, letting $\mathbf S$ be the scaled adjacency matrix of a Ramanujan graph with $\tilde O(n/\epsilon^4)$ edges, we have $\|\mathbf A - \mathbf A \circ \mathbf S \|_2 \le \epsilon \cdot \max(n,\|\mathbf A\|_1)$, where $\|\mathbf A\|_1$ is the nuclear norm. We show that the above bounds for both PSD and non-PSD matrices are tight up to log factors. Since $\mathbf A \circ \mathbf S$ can be constructed deterministically, our result for PSD matrices derandomizes and improves upon known results for randomized matrix sparsification, which require randomly sampling ${O}(\frac{n \log n}{\epsilon^2})$ entries. We also leverage our results to give the first deterministic algorithms for several problems related to singular value approximation that run in faster than matrix multiplication time. Finally, if $\mathbf A \in \{-1,0,1\}^{n \times n}$ is PSD, we show that $\mathbf{\tilde A}$ with $\|\mathbf A - \mathbf{\tilde A}\|_2 \le \epsilon n$ can be obtained by deterministically reading $\tilde O(n/\epsilon)$ entries of $\mathbf A$. This improves the $1/\epsilon$ dependence on our result for general PSD matrices and is near-optimal., Comment: 41 pages

Artificial neural networks (ANNs) are powerful machine learning methods used in many modern applications such as facial recognition, machine translation, and cancer diagnostics. A common issue with ANNs is that they usually have millions or billions of trainable parameters, and therefore tend to overfit to the training data. This is especially problematic in applications where it is important to have reliable uncertainty estimates. Bayesian neural networks (BNN) can improve on this, since they incorporate parameter uncertainty. In addition, latent binary Bayesian neural networks (LBBNN) also take into account structural uncertainty by allowing the weights to be turned on or off, enabling inference in the joint space of weights and structures. In this paper, we will consider two extensions to the LBBNN method: Firstly, by using the local reparametrization trick (LRT) to sample the hidden units directly, we get a more computationally efficient algorithm. More importantly, by using normalizing flows on the variational posterior distribution of the LBBNN parameters, the network learns a more flexible variational posterior distribution than the mean field Gaussian. Experimental results show that this improves predictive power compared to the LBBNN method, while also obtaining more sparse networks. We perform two simulation studies. In the first study, we consider variable selection in a logistic regression setting, where the more flexible variational distribution leads to improved results. In the second study, we compare predictive uncertainty based on data generated from two-dimensional Gaussian distributions. Here, we argue that our Bayesian methods lead to more realistic estimates of predictive uncertainty., Comment: 24 pages, 10 figures

Bayesian neural networks (BNNs) have recently regained a significant amount of attention in the deep learning community due to the development of scalable approximate Bayesian inference techniques. There are several advantages of using a Bayesian approach: Parameter and prediction uncertainties become easily available, facilitating rigorous statistical analysis. Furthermore, prior knowledge can be incorporated. However, so far, there have been no scalable techniques capable of combining both structural and parameter uncertainty. In this paper, we apply the concept of model uncertainty as a framework for structural learning in BNNs and hence make inference in the joint space of structures/models and parameters. Moreover, we suggest an adaptation of a scalable variational inference approach with reparametrization of marginal inclusion probabilities to incorporate the model space constraints. Experimental results on a range of benchmark datasets show that we obtain comparable accuracy results with the competing models, but based on methods that are much more sparse than ordinary BNNs., Comment: arXiv admin note: text overlap with arXiv:1903.07594

For simulating incompressible flows by projection methods. it is generally accepted that the pressure-correction stage is the most time-consuming part of the flow solver. The objective of the present work is to develop a fast hybrid pressure-correction algorithm for numerical simulation of incompressible flows around obstacles in the context of projection methods. The key idea is to adopt different numerical methods/discretizations in the sub-steps of projection methods. Here, a classical second-order time-marching projection method which consists of two sub-steps is chosen for the purpose of demonstration. In the first sub-step, the momentum equations are discretized on unstructured grids and solved by conventional numerical methods, here, a meshless method. In the second sub-step (pressure-correction), the proposed algorithm adopts a double discretization system and combines the weighted least squares approximation with the essence of immersed boundary methods. Such a design allows us to develop a FFT-based solver to speed up the solution of the pressure Poisson equation for flow cases with obstacles, while keeping the implementation of boundary conditions for the momentum equations as easy as conventional numerical methods do with unstructured grids. Numerical experiments of five test cases have been performed to verify and validate the proposed hybrid algorithm and evaluate its computational performance. The results show that the new FFT-based hybrid algorithm is working and robust, and it is significantly faster than the multigrid-based reference method. The hybrid algorithm opens an avenue for the development of next-generation high-performance parallel computational fluid dynamics solvers for incompressible flows., Comment: 40 pages, 21 figures

We show that a very simple randomised algorithm for numerical integration can produce a near optimal rate of convergence for integrals of functions in the $d$-dimensional weighted Korobov space. This algorithm uses a lattice rule with a fixed generating vector and the only random element is the choice of the number of function evaluations. For a given computational budget $n$ of a maximum allowed number of function evaluations, we uniformly pick a prime $p$ in the range $n/2 < p \le n$. We show error bounds for the randomised error, which is defined as the worst case expected error, of the form $O(n^{-\alpha - 1/2 + \delta})$, with $\delta > 0$, for a Korobov space with smoothness $\alpha > 1/2$ and general weights. The implied constant in the bound is dimension-independent given the usual conditions on the weights. We present an algorithm that can construct suitable generating vectors \emph{offline} ahead of time at cost $O(d n^4 / \ln n)$ when the weight parameters defining the Korobov spaces are so-called product weights. For this case, numerical experiments confirm our theory that the new randomised algorithm achieves the near optimal rate of the randomised error., Comment: 29 pages, 2 figures

Given two nonempty and disjoint intersections of closed and convex subsets, we look for a best approximation pair relative to them, i.e., a pair of points, one in each intersection, attaining the minimum distance between the disjoint intersections. We propose an iterative process based on projections onto the subsets which generate the intersections. The process is inspired by the Halpern-Lions-Wittmann-Bauschke algorithm and the classical alternating process of Cheney and Goldstein, and its advantage is that there is no need to project onto the intersections themselves, a task which can be rather demanding. We prove that under certain conditions the two interlaced subsequences converge to a best approximation pair. These conditions hold, in particular, when the space is Euclidean and the subsets which generate the intersections are compact and strictly convex. Our result extends the one of Aharoni, Censor and Jiang ["Finding a best approximation pair of points for two polyhedra'', Computational Optimization and Applications 71 (2018), 509--523] who considered the case of finite-dimensional polyhedra., Comment: Added Mathematics Subject Classification, DOI and journal reference, a few additional very slight modifications

Graphs arising in statistical problems, signal processing, large networks, combinatorial optimization, and data analysis are often dense, which causes both computational and storage bottlenecks. One way of \textit{sparsifying} a \textit{weighted} graph, while sharing the same vertices as the original graph but reducing the number of edges, is through \textit{spectral sparsification}. We study this problem through the perspective of RandNLA. Specifically, we utilize randomized matrix multiplication to give a clean and simple analysis of how sampling according to edge weights gives a spectral approximation to graph Laplacians. Through the $CR$-MM algorithm, we attain a simple and computationally efficient sparsifier whose resulting Laplacian estimate is unbiased and of minimum variance. Furthermore, we define a new notion of \textit{additive spectral sparsifiers}, which has not been considered in the literature.

As we reach the limit of Moore's Law, researchers are exploring different paradigms to achieve unprecedented performance. Approximate Computing (AC), which relies on the ability of applications to tolerate some error in the results to trade-off accuracy for performance, has shown significant promise. Despite the success of AC in domains such as Machine Learning, its acceptance in High-Performance Computing (HPC) is limited due to stringent requirements for accuracy. We need tools and techniques to identify regions of code that are amenable to approximations and their impact on the application output quality to guide developers to employ selective approximation. To this end, we propose CHEF-FP, a flexible, scalable, and easy-to-use source-code transformation tool based on Automatic Differentiation (AD) for analyzing approximation errors in HPC applications. CHEF-FP uses Clad, an efficient AD tool built as a plugin to the Clang compiler and based on the LLVM compiler infrastructure, as a backend and utilizes its AD abilities to evaluate approximation errors in C++ code. CHEF-FP works at the source by injecting error estimation code into the generated adjoints. This enables the error-estimation code to undergo compiler optimizations resulting in improved analysis time and reduced memory usage. We also provide theoretical and architectural augmentations to source code transformation-based AD tools to perform FP error analysis. This paper primarily focuses on analyzing errors introduced by mixed-precision AC techniques. We also show the applicability of our tool in estimating other kinds of errors by evaluating our tool on codes that use approximate functions. Moreover, we demonstrate the speedups CHEF-FP achieved during analysis time compared to the existing state-of-the-art tool due to its ability to generate and insert approximation error estimate code directly into the derivative source., Comment: 11 pages, to appear in the 2023 IEEE International Parallel and Distributed Processing Symposium (IPDPS'23)

Recent developments in surrogate construction predominantly focused on two strategies to improve surrogate accuracy. Firstly, component-wise domain scaling informed by cross-validation. Secondly, regression to construct response surfaces using additional information in the form of additional function-values sampled from multi-fidelity models and gradients. Component-wise domain scaling reliably improves the surrogate quality at low dimensions but has been shown to suffer from high computational costs for higher dimensional problems. The second strategy, adding gradients to train surrogates, typically results in regression surrogates. Counter-intuitively, these gradient-enhanced regression-based surrogates do not exhibit improved accuracy compared to surrogates only interpolating function values. This study empirically establishes three main findings. Firstly, constructing the surrogate in poorly scaled domains is the predominant cause of deteriorating response surfaces when regressing with additional gradient information. Secondly, surrogate accuracy improves if the surrogates are constructed in a fully transformed domain, by scaling and rotating the original domain, not just simply scaling the domain. The domain transformation scheme should be based on the local curvature of the approximation surface and not its global curvature. Thirdly, the main benefit of gradient information is to efficiently determine the (near) optimal domain in which to construct the surrogate. This study proposes a foundational transformation algorithm that performs near-optimal transformations for lower dimensional problems. The algorithm consistently outperforms cross-validation-based component-wise domain scaling for higher dimensional problems. A carefully selected test problem set that varies between 2 and 16-dimensional problems is used to clearly demonstrate the three main findings of this study., Comment: 20 pages, 28 figures

We consider the problem of state estimation from a few linear measurements, where the state to recover is an element of the manifold $\mathcal{M}$ of solutions of a parameter-dependent equation. The state is estimated using prior knowledge on $\mathcal{M}$ coming from model order reduction. Variational approaches based on linear approximation of $\mathcal{M}$, such as PBDW, yields a recovery error limited by the Kolmogorov width of $\mathcal{M}$. To overcome this issue, piecewise-affine approximations of $\mathcal{M}$ have also been considered, that consist in using a library of linear spaces among which one is selected by minimizing some distance to $\mathcal{M}$. In this paper, we propose a state estimation method relying on dictionary-based model reduction, where a space is selected from a library generated by a dictionary of snapshots, using a distance to the manifold. The selection is performed among a set of candidate spaces obtained from a set of $\ell_1$-regularized least-squares problems. Then, in the framework of parameter-dependent operator equations (or PDEs) with affine parametrizations, we provide an efficient offline-online decomposition based on randomized linear algebra, that ensures efficient and stable computations while preserving theoretical guarantees., Comment: 29 pages, 7 figures

Parametric reduced-order modelling often serves as a surrogate method for hemodynamics simulations to improve the computational efficiency in many-query scenarios or to perform real-time simulations. However, the snapshots of the method require to be collected from the same discretisation, which is a straightforward process for physical parameters, but becomes challenging for geometrical problems, especially for those domains featuring unparameterised and unique shapes, e.g. patient-specific geometries. In this work, a data-driven surrogate model is proposed for the efficient prediction of blood flow simulations on similar but distinct domains. The proposed surrogate model leverages group surface registration to parameterise those shapes and formulates corresponding hemodynamics information into geometry-informed snapshots by the diffeomorphisms constructed between a reference domain and original domains. A non-intrusive reduced-order model for geometrical parameters is subsequently constructed using proper orthogonal decomposition, and a radial basis function interpolator is trained for predicting the reduced coefficients of the reduced-order model based on compressed geometrical parameters of the shape. Two examples of blood flowing through a stenosis and a bifurcation are presented and analysed. The proposed surrogate model demonstrates its accuracy and efficiency in hemodynamics prediction and shows its potential application toward real-time simulation or uncertainty quantification for complex patient-specific scenarios.

We consider a class of density-driven flow problems. We are particularly interested in the problem of the salinization of coastal aquifers. We consider the Henry saltwater intrusion problem with uncertain porosity, permeability, and recharge parameters as a test case. The reason for the presence of uncertainties is the lack of knowledge, inaccurate measurements, and inability to measure parameters at each spatial or time location. This problem is nonlinear and time-dependent. The solution is the salt mass fraction, which is uncertain and changes in time. Uncertainties in porosity, permeability, recharge, and mass fraction are modeled using random fields. This work investigates the applicability of the well-known multilevel Monte Carlo (MLMC) method for such problems. The MLMC method can reduce the total computational and storage costs. Moreover, the MLMC method runs multiple scenarios on different spatial and time meshes and then estimates the mean value of the mass fraction. The parallelization is performed in both the physical space and stochastic space. To solve every deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion. We use the solution obtained from the quasi-Monte Carlo method as a reference solution., Comment: 24 pages, 3 tables, 11 figures

Stability and optimal convergence analysis of a non-uniform implicit-explicit L1 finite element method (IMEX-L1-FEM) is studied for a class of time-fractional linear partial differential/integro-differential equations with non-self-adjoint elliptic part having (space-time) variable coefficients. The proposed scheme is based on a combination of an IMEX-L1 method on graded mesh in the temporal direction and a finite element method in the spatial direction. With the help of a discrete fractional Gr\"{o}nwall inequality, global almost optimal error estimates in $L^2$- and $H^1$-norms are derived for the problem with initial data $u_0 \in H_0^1(\Omega)\cap H^2(\Omega)$. The novelty of our approach is based on managing the interaction of the L1 approximation of the fractional derivative and the time discrete elliptic operator to derive the optimal estimate in $H^1$-norm directly. Furthermore, a super convergence result is established when the elliptic operator is self-adjoint with time and space varying coefficients, and as a consequence, an $L^\infty$ error estimate is obtained for 2D problems that too with the initial condition is in $ H_0^1(\Omega)\cap H^2(\Omega)$. All results proved in this paper are valid uniformly as $\alpha\longrightarrow 1^{-}$, where $\alpha$ is the order of the Caputo fractional derivative. Numerical experiments are presented to validate our theoretical findings., Comment: In the revised version, we have consistently maintained the assumption that $u_0 \in H_0^1(\Omega)\cap H^2(\Omega)$, given the impracticality of assuming $u_0 \in \dot{H}^3(\Omega)$. It is important to highlight that all the derived estimates and bounds in this revised version remain valid as $\alpha\to 1^{-}$, where $\alpha$ denotes the order of the Caputo fractional derivative

Global seismicity on all three solar system's bodies with in situ measurements (Earth, Moon, and Mars) is due mainly to mechanical Rieger resonance (RR) of the solar wind's macroscopic flapping, driven by the well-known PRg=~154-day Rieger period and detected commonly in most heliophysical data types and the interplanetary magnetic field (IMF). Thus, InSight mission marsquakes rates are periodic with PRg as characterized by a very high (>>12) fidelity {\Phi}=2.8 10^6 and by being the only >99%-significant spectral peak in the 385.8-64.3-nHz (1-180-day) band of highest planetary energies; the longest-span (v.9) release of raw data revealed the entire RR, excluding a tectonically active Mars. For check, I analyze rates of Oct 2015-Feb 2019, Mw5.6+ earthquakes, and all (1969-1977) Apollo mission moonquakes. To decouple magnetosphere and IMF effects, I study Earth and Moon seismicity during traversals of Earth magnetotail vs. IMF. The analysis showed with >99-67% confidence and {\Phi}>>12 fidelity that (an unspecified majority of) moonquakes and Mw5.6+ earthquakes also recur at Rieger periods. About half of spectral peaks split but also into clusters that average to usual Rieger periodicities, where magnetotail reconnecting clears the signal. Moonquakes are mostly forced at times of solar-wind resonance and not just during tides, as previously and simplistically believed. Earlier claims that solar plasma dynamics could be seismogenic are confirmed. This result calls for reinterpreting the seismicity phenomenon and for reliance on global magnitude scales. Predictability of solar-wind macroscopic dynamics is now within reach, which paves the way for long-term physics-based seismic and space weather prediction and the safety of space missions. Gauss-Vanicek Spectral Analysis revolutionizes geophysics by computing nonlinear global dynamics directly (renders approximating of dynamics obsolete)., Comment: 46 printed pages, 8 figures, 5 tables; includes time-series of all Apollo natural moonquakes (1969-1977), all Mw5.6+ natural earthquakes between Oct 2015-Feb 2019, and all v.9 InSight natural marsquakes (Jan 2019-Sep 2021); pre-publishing proof version. Press Release: https://www.openpr.com/news/2982920.html

Rather than as a star classically assumed to feature elusive dynamo or a proverbial engine and impulsively alternating polarity, the Sun reveals itself in the 385.8-2.439-nHz (1-mo-13-yr) band of polar ({\phi}Sun>|70{\deg}|) wind's decadal dynamics, dominated by fast (>700 km s^-1) winds, as a globally completely vibrating revolving-field magnetic alternator at work at all times. Thus North-South separation of 1994-2008 Ulysses <10 nT wind polar samplings reveals spectral signatures of an entirely >99%-significant Sun-borne global sharp Alfven resonance (AR) governed by Ps=~11-yr Schwabe global damping (equilibrium) mode northside, ~10-yr degeneration equatorially, and ~9-yr southside. AR is accompanied by a symmetrical antiresonance P(-) whose both N/S tailing harmonics P(-17) are the well-known PR=~154-day (Ps/3/3/3 to +-0.1%) Rieger period dominating solar-system (planetary) dynamics and space weather. Upward drifting low-frequency trends reveal a rigid core; undertones a core offset away from apex. Core wobble at 2.2+/-0.1-yr triggers AR. Multiple total spectral symmetries of solar activity (historical solar-cycle lengths, sunspot and calcium numbers) expose the solar alternator moderating sunspots, nanoflares and CMEs that resemble machinery sparking. Unlike a resonating motor restrained from separating its casing, the cageless Sun lacks a stator and vibrates freely, resulting in all-spin and mass release (fast winds) in an axial shake-off beyond L1. The result was verified against remote data and the experiment, replacing dynamo with magnetoalternator and advancing basic knowledge on the >100 billion trillions solar-type stars. Shannon's theory-based Gauss-Vanicek spectral analysis revolutionizes astrophysics by rigorously simulating fleet formations from a single spacecraft and physics by computing nonlinear global dynamics directly (rendering spherical approximation obsolete)., Comment: 30 printed pages, 14 figures, 6 tables; pre-publishing proof version. Press Release: https://www.openpr.com/news/3337049.html

The problem faced in this paper concerns the recognition of simple and complex geometric primitives in point clouds resulting from scans of mechanical CAD objects. A large number of points, the presence of noise, outliers, missing or redundant parts and uneven distribution are the main problems to be addressed to meet this need. In this article we propose a solution, based on the Hough transform, that can recognize simple and complex geometric primitives and is robust to noise, outliers, and missing parts. Additionally, we can extract a series of geometric descriptors that uniquely characterize a primitive and, based on them, aggregate the output into maximal or compound primitives, thus reducing oversegmentation. The results presented in the paper demonstrate the robustness of the method and its competitiveness with respect to other solutions proposed in the literature.

Coded computing is a method for mitigating straggling workers in a centralized computing network, by using erasure-coding techniques. Federated learning is a decentralized model for training data distributed across client devices. In this work we propose approximating the inverse of an aggregated data matrix, where the data is generated by clients; similar to the federated learning paradigm, while also being resilient to stragglers. To do so, we propose a coded computing method based on gradient coding. We modify this method so that the coordinator does not access the local data at any point; while the clients access the aggregated matrix in order to complete their tasks. The network we consider is not centrally administrated, and the communications which take place are secure against potential eavesdroppers., Comment: arXiv admin note: substantial text overlap with arXiv:2207.06271

This paper deals with the derivation of Non-Intrusive Reduced Basis (NIRB) techniques for sensitivity analysis, more specifically the direct and adjoint state methods. For highly complex parametric problems, these two approaches may become too costly. To reduce computational times, Proper Orthogonal Decomposition (POD) and Reduced Basis Methods (RBMs) have already been investigated. The majority of these algorithms are however intrusive in the sense that the High-Fidelity (HF) code must be modified. To address this issue, non-intrusive strategies are employed. The NIRB two-grid method uses the HF code solely as a ``black-box'', requiring no code modification. Like other RBMs, it is based on an offline-online decomposition. The offline stage is time-consuming, but it is only executed once, whereas the online stage is significantly less expensive than an HF evaluation. In this paper, we propose new NIRB two-grid algorithms for both the direct and adjoint state methods. On the direct method, we prove on a classical model problem, the heat equation, that HF evaluations of sensitivities reach an optimal convergence rate in $L^{\infty}(0,T;H^1(\Omega))$, and then establish that these rates are recovered by the proposed NIRB approximation. These results are supported by numerical simulations. We then numerically demonstrate that a Gaussian process regression can be used to approximate the projection coefficients of the NIRB two-grid method. This further reduces the computational costs of the online step while only computing a coarse solution of the initial problem. All numerical results are run with the model problem as well as a more complex problem, namely the Brusselator system., Comment: 37 pages, 3 figures. arXiv admin note: text overlap with arXiv:2211.08897

Claims of paleodata periodicity are many and so controversial that superimposing Phanerozoic (0-541 My) mass-extinction periods renders life on Earth impossible. This period hunt coincided with geochronology modernization tying geological timescales to orbital frequencies, enabling separation of astronomical signals from harmonics. I thus show on diverse data (geomagnetic polarity, cratering, extinction episodes) as a proxy of planetary paleodynamics that many-body subharmonic entrainment induces Earth's resonant response to astronomical forcing so that 2{\pi}-phase-shifted axial precession p=26 ky and its Pi=2{\pi}p/i; i=1...n harmonics get resonantly responsible for paleodata periodicity. This quasiperiodic nature of strata is co-triggered by a p'/4-lockstep to p'=41-ky obliquity. For verification, residuals analysis after suppressing 2{\pi}p (and so Pi) in GPTS-95 polarity reversals timescale's calibration at the South Atlantic Anomaly, extending to end-Campanian (0-83 My), successfully detected weak signals of Earth-Mars planetary resonances reported previously from older epochs. The only residual signal is 26.5-My Rampino period - carrier wave of crushing deflections and transformative reversals. While the Earth's (2{\pi}p, Pi) resonant response is the long-sought energy transfer mechanism of Milankovitch theory, fundamental system properties - 2{\pi}-phase-shift, 1/4 lockstep to a forcer, and discrete time translation symmetry (multiplied/halved periods) - typical of a quantum time crystal, appear macroscopic. The surprising cross-scale outcome confirms planetary precession is a cataclysmic geodynamic phenomenon claimed previously as the mechanism for Earth's reversals and expansion. As the Earth-Moon-Sun resonant system sets the tone for all subscale resonances like the quantum, quantum physics in any interpretation is fundamentally different and variable outside Earth's vicinity., Comment: 37 pages, 6 figures, 4 tables; conclusions expanded, abstract and highlights updated, references expanded and updated, typos corrected

We develop two new sets of stable, rank-adaptive Dynamically Orthogonal Runge-Kutta (DORK) schemes that capture the high-order curvature of the nonlinear low-rank manifold. The DORK schemes asymptotically approximate the truncated singular value decomposition at a greatly reduced cost while preserving mode continuity using newly derived retractions. We show that arbitrarily high-order optimal perturbative retractions can be obtained, and we prove that these new retractions are stable. In addition, we demonstrate that repeatedly applying retractions yields a gradient-descent algorithm on the low-rank manifold that converges superlinearly when approximating a low-rank matrix. When approximating a higher-rank matrix, iterations converge linearly to the best low-rank approximation. We then develop a rank-adaptive retraction that is robust to overapproximation. Building off of these retractions, we derive two rank-adaptive integration schemes that dynamically update the subspace upon which the system dynamics are projected within each time step: the stable, optimal Dynamically Orthogonal Runge-Kutta (so-DORK) and gradient-descent Dynamically Orthogonal Runge-Kutta (gd-DORK) schemes. These integration schemes are numerically evaluated and compared on an ill-conditioned matrix differential equation, an advection-diffusion partial differential equation, and a nonlinear, stochastic reaction-diffusion partial differential equation. Results show a reduced error accumulation rate with the new stable, optimal and gradient-descent integrators. In addition, we find that rank adaptation allows for highly accurate solutions while preserving computational efficiency., Comment: 29 pages, 8 figures

An autonomous robot with a limited vision range finds a path to the goal in an unknown environment in 2D avoiding polygonal obstacles. In the process of discovering the environmental map, the robot has to return to some positions marked previously, the regions where the robot traverses to return are defined as sequences of bundles of line segments. This paper presents a novel algorithm for finding approximately shortest paths along the sequences of bundles of line segments based on the method of multiple shooting. Three factors of the approach including bundle partition, collinear condition, and update of shooting points are presented. We then show that if the collinear condition holds, the exactly shortest paths of the problems are determined, otherwise, the sequence of paths obtained by the update of the method converges to the shortest path. The algorithm is implemented in Python and some numerical examples show that the running time of path-planning for autonomous robots using our method is faster than that using the rubber band technique of Li and Klette in Euclidean Shortest Paths, Springer, 53-89 (2011)., Comment: 26 pages, 39 figures

We present two new algorithms for Householder QR factorization of Block Low-Rank (BLR) matrices: one that performs block-column-wise QR, and another that is based on tiled QR. We show how the block-column-wise algorithm exploits BLR structure to achieve arithmetic complexity of $\mathcal{O}(mn)$, while the tiled BLR-QR exhibits $\mathcal{O}(mn^{1.5})$ complexity. However, the tiled BLR-QR has finer task granularity that allows parallel task-based execution on shared memory systems. We compare the block-column-wise BLR-QR using fork-join parallelism with tiled BLR-QR using task-based parallelism. We also compare these two implementations of Householder BLR-QR with a block-column-wise Modified Gram-Schmidt (MGS) BLR-QR using fork-join parallelism, and a state-of-the-art vendor-optimized dense Householder QR in Intel MKL. For a matrix of size 131k $\times$ 65k, all BLR methods are more than an order of magnitude faster than the dense QR in MKL. Our methods are also robust to ill-conditioning and produce better orthogonal factors than the existing MGS-based method. On a CPU with 64 cores, our parallel tiled Householder and block-column-wise Householder algorithms show a speedup of 50 and 37 times, respectively., Comment: 23 pages, 12 figures

The Carleman linearization is one of the mainstream approaches to lift a finite-dimensional nonlinear dynamical system into an infinite-dimensional linear system with the promise of providing accurate approximations of the original nonlinear system over larger regions around the equilibrium for longer time horizons with respect to the conventional first-order linearization approach. Finite-section approximations of the lifted system has been widely used to study dynamical and control properties of the original nonlinear system. In this context, some of the outstanding problems are to determine under what conditions, as the finite-section order (i.e., truncation length) increases, the trajectory of the resulting approximate linear system from the finite-section scheme converges to that of the original nonlinear system and whether the time interval over which the convergence happens can be quantified explicitly. In this paper, we provide explicit error bounds for the finite-section approximation and prove that the convergence is indeed exponential with respect to the finite-section order. For a class of nonlinear systems, it is shown that one can achieve exponential convergence over the entire time horizon up to infinity. Our results are practically plausible as our proposed error bound estimates can be used to compute proper truncation lengths for a given application, e.g., determining proper sampling period for model predictive control and reachability analysis for safety verifications. We validate our theoretical findings through several illustrative simulations., Comment: 25 Pages, 10 figures

A cumbersome operation in many scientific fields, is inverting large full-rank matrices. In this paper, we propose a coded computing approach for recovering matrix inverse approximations. We first present an approximate matrix inversion algorithm which does not require a matrix factorization, but uses a black-box least squares optimization solver as a subroutine, to give an estimate of the inverse of a real full-rank matrix. We then present a distributed framework for which our algorithm can be implemented, and show how we can leverage sparsest-balanced MDS generator matrices to devise matrix inversion coded computing schemes. We focus on balanced Reed-Solomon codes, which are optimal in terms of computational load; and communication from the workers to the master server. We also discuss how our algorithms can be used to compute the pseudoinverse of a full-rank matrix, and how the communication is secured from eavesdroppers.

The paper proposes an algorithm for a discretization (sampled-time implementation) of a homogeneous control preserving the finite-time and nearly fixed-time stability property of the original (sampling-free) system. The sampling period is assumed to be constant. Both single-input and multiple-input cases are considered. The robustness (Input-to-State Stability) of the obtained sampled-time control system is studied as well. Theoretical results are supported by numerical simulations.

This paper explores how to identify a reduced order model (ROM) from a physical system. A ROM captures an invariant subset of the observed dynamics. We find that there are four ways a physical system can be related to a mathematical model: invariant foliations, invariant manifolds, autoencoders and equation-free models. Identification of invariant manifolds and equation-free models require closed-loop manipulation of the system. Invariant foliations and autoencoders can also use off-line data. Only invariant foliations and invariant manifolds can identify ROMs, the rest identify complete models. Therefore, the common case of identifying a ROM from existing data can only be achieved using invariant foliations. Finding an invariant foliation requires approximating high-dimensional functions. For function approximation, we use polynomials with compressed tensor coefficients, whose complexity increases linearly with increasing dimensions. An invariant manifold can also be found as the fixed leaf of a foliation. This only requires us to resolve the foliation in a small neighbourhood of the invariant manifold, which greatly simplifies the process. Combining an invariant foliation with the corresponding invariant manifold provides an accurate ROM. We analyse the ROM in case of a focus type equilibrium, typical in mechanical systems. The nonlinear coordinate system defined by the invariant foliation or the invariant manifold distorts instantaneous frequencies and damping ratios, which we correct. Through examples we illustrate the calculation of invariant foliations and manifolds, and at the same time show that Koopman eigenfunctions and autoencoders fail to capture accurate ROMs under the same conditions., Comment: 48 pages, 16 figures. Update: some bugs were fixed in the numerical calculations

A variety of optimization problems takes the form of a minimum norm optimization. In this paper, we study the change of optimal values between two incrementally constructed least norm optimization problems, with new measurements included in the second one. We prove an exact equation to calculate the change of optimal values in the linear least norm optimization problem. With the result in this paper, the change of the optimal values can be pre-calculated as a metric to guide online decision makings, without solving the second optimization problem as long the solution and covariance of the first optimization problem are available. The result can be extended to linear least distance optimization problems, and nonlinear least distance optimization with (nonlinear) equality constraints through linearizations. This derivation in this paper provides a theoretically sound explanation to the empirical observations shown in RA-L 2018 bai et al. As an additional contribution, we propose another optimization problem, i.e. aligning two trajectories at given poses, to further demonstrate how to use the metric. The accuracy of the metric is validated with numerical examples, which is quite satisfactory in general (see the experiments in RA-L 2018 bai et al.} as well), unless in some extremely adverse scenarios. Last but not least, calculating the optimal value by the proposed metric is at least one magnitude faster than solving the corresponding optimization problems directly., Comment: 6 pages on IEEE International Conference on Robotics and Automation (ICRA), 2020

This paper presents the methods that have participated in the SHREC 2022 track on the fitting and recognition of simple geometric primitives on point clouds. As simple primitives we mean the classical surface primitives derived from constructive solid geometry, i.e., planes, spheres, cylinders, cones and tori. The aim of the track is to evaluate the quality of automatic algorithms for fitting and recognising geometric primitives on point clouds. Specifically, the goal is to identify, for each point cloud, its primitive type and some geometric descriptors. For this purpose, we created a synthetic dataset, divided into a training set and a test set, containing segments perturbed with different kinds of point cloud artifacts. Among the six participants to this track, two are based on direct methods, while four are either fully based on deep learning or combine direct and neural approaches. The performance of the methods is evaluated using various classification and approximation measures.

The automatic creation of geometric models from point clouds has numerous applications in CAD (e.g., reverse engineering, manufacturing, assembling) and, more in general, in shape modelling and processing. Given a segmented point cloud representing a man-made object, we propose a method for recognizing simple geometric primitives and their interrelationships. Our approach is based on the Hough transform (HT) for its ability to deal with noise, missing parts and outliers. In our method we introduce a novel technique for processing segmented point clouds that, through a voting procedure, is able to provide an initial estimate of the geometric parameters characterizing each primitive type. By using these estimates, we localize the search of the optimal solution in a dimensionally-reduced parameter space thus making it efficient to extend the HT to more primitives than those that are generally found in the literature, i.e. planes and spheres. Then, we extract a number of geometric descriptors that uniquely characterize a segment, and, on the basis of these descriptors, we show how to aggregate parts of primitives (segments). Experiments on both synthetic and industrial scans reveal the robustness of the primitive fitting method and its effectiveness for inferring relations among segments.

The exponential function maps the imaginary axis to the unit circle and, for many applications, this unitarity property is also desirable from its approximations. We show that this property is conserved not only by the (k,k)-rational barycentric interpolant of the exponential on the imaginary axis, but also by (k,k)-rational barycentric approximants that minimize a linearized approximation error. These results are a consequence of certain properties of singular vectors of Loewner-type matrices associated to linearized approximation errors. Prominent representatives of this class are rational approximants computed by the adaptive Antoulas--Anderson (AAA) method and the AAA--Lawson method. Our results also lead to a modified procedure with improved numerical stability of the unitarity property and reduced computational cost., Comment: 23 pages, 2 figures

The block-iterative projections (BIP) method of Aharoni and Censor [Block-iterative projection methods for parallel computation of solutions to convex feasibility problems, Linear Algebra and its Applications 120, (1989), 165--175] is an iterative process for finding asymptotically a point in the nonempty intersection of a family of closed convex subsets. It employs orthogonal projections onto the individual subsets in an algorithmic regime that uses "blocks" of operators and has great flexibility in constructing specific algorithms from it. We extend this algorithmic scheme to handle a family of continuous cutter operators and to find a common fixed point of them. Since the family of continuous cutters includes several important specific operators, our generalized scheme, which ensures global convergence and retains the flexibility of BIP, can handle, in particular, metric (orthogonal) projectors and continuous subgradient projections, which are very important in applications. We also allow a certain kind of adaptive perturbations to be included, and along the way we derive a perturbed Fej\'er monotonicity lemma which is of independent interest., Comment: Minor typographical corrections; added DOI; to appear in the Journal of Global Optimization

Data are often sampled irregularly in time. Dealing with this using Recurrent Neural Networks (RNNs) traditionally involved ignoring the fact, feeding the time differences as additional inputs, or resampling the data. All these methods have their shortcomings. We propose an elegant straightforward alternative approach where instead the RNN is in effect resampled in time to match the time of the data or the task at hand. We use Echo State Network (ESN) and Gated Recurrent Unit (GRU) as the basis for our solution. Such RNNs can be seen as discretizations of continuous-time dynamical systems, which gives a solid theoretical ground to our approach. Our Task-Synchronized ESN (TSESN) and GRU (TSGRU) models allow for a direct model time setting and require no additional training, parameter tuning, or computation (solving differential equations or interpolating data) compared to their regular counterparts, thus retaining their original efficiency. We confirm empirically that our models can effectively compensate for the time-non-uniformity of the data and demonstrate that they compare favorably to data resampling, classical RNN methods, and alternative RNN models proposed to deal with time irregularities on several real-world nonuniform-time datasets. We open-source the code at https://github.com/oshapio/task-synchronized-RNNs ., Comment: The 1st version was written in May 2019 and double-blind reviewed for a prominent conference. A major update. We changed the name of the article and methods to an arguably more precise one, and because a very similar title has been published in the meantime. We've rewritten much of the text, connected to the current literature, redone some experiments, figures, discussion, published source code

This work deals with the design optimization of electrical machines under the consideration of manufacturing uncertainties. In order to efficiently quantify the uncertainty, blackbox machine learning methods are employed. A multi-objective optimization problem is formulated, maximizing simultaneously the reliability, i.e., the yield, and further performance objectives, e.g., the costs. A permanent magnet synchronous machine is modeled and simulated in commercial finite element simulation software. Four approaches for solving the multi-objective optimization problem are described and numerically compared, namely: epsilon-constraint scalarization, weighted sum scalarization, a multi-start weighted sum approach and a genetic algorithm.

Significant effort has been made to solve computationally expensive optimization problems in the past two decades, and various optimization methods incorporating surrogates into optimization have been proposed. Most research focuses on either exploiting the surrogate by defining a utility optimization problem or customizing an existing optimization method to use one or multiple approximation models. However, only a little attention has been paid to generic concepts applicable to different types of algorithms and optimization problems simultaneously. Thus this paper proposes a generalized probabilistic surrogate-assisted framework (GPSAF), applicable to a broad category of unconstrained and constrained, single- and multi-objective optimization algorithms. The idea is based on a surrogate assisting an existing optimization method. The assistance is based on two distinct phases, one facilitating exploration and another exploiting the surrogates. The exploration and exploitation of surrogates are automatically balanced by performing a probabilistic knockout tournament among different clusters of solutions. A study of multiple well-known population-based optimization algorithms is conducted with and without the proposed surrogate assistance on single- and multi-objective optimization problems with a maximum solution evaluation budget of 300 or less. The results indicate the effectiveness of applying GPSAF to an optimization algorithm and the competitiveness with other surrogate-assisted algorithms.

B-spline models are a powerful way to represent scientific data sets with a functional approximation. However, these models can suffer from spurious oscillations when the data to be approximated are not uniformly distributed. Model regularization (i.e., smoothing) has traditionally been used to minimize these oscillations; unfortunately, it is sometimes impossible to sufficiently remove unwanted artifacts without smoothing away key features of the data set. In this article, we present a method of model regularization that preserves significant features of a data set while minimizing artificial oscillations. Our method varies the strength of a smoothing parameter throughout the domain automatically, removing artifacts in poorly-constrained regions while leaving other regions unchanged. The behavior of our method is validated on a collection of two- and three-dimensional data sets produced by scientific simulations., Comment: To appear in Proceedings of the International Conference on Computational Science, 2022. 15 pages, 5 figures

Rectified-linear-unit (ReLU) neural networks, which play a prominent role in deep learning, generate continuous and piecewise-linear (CPWL) functions. While they provide a powerful parametric representation, the mapping between the parameter and function spaces lacks stability. In this paper, we investigate an alternative representation of CPWL functions that relies on local hat basis functions. It is predicated on the fact that any CPWL function can be specified by a triangulation and its values at the grid points. We give the necessary and sufficient condition on the triangulation (in any number of dimensions) for the hat functions to form a Riesz basis, which ensures that the link between the parameters and the corresponding CPWL function is stable and unique. In addition, we provide an estimate of the $\ell_2\rightarrow L_2$ condition number of this local representation. Finally, as a special case of our framework, we focus on a systematic parametrization of $\mathbb{R}^d$ with control points placed on a uniform grid. In particular, we choose hat basis functions that are shifted replicas of a single linear box spline. In this setting, we prove that our general estimate of the condition number is optimal. We also relate our local representation to a nonlocal one based on shifts of a causal ReLU-like function.

It is common practice to use Laplace approximations to compute marginal likelihoods in Bayesian versions of generalised linear models (GLM). Marginal likelihoods combined with model priors are then used in different search algorithms to compute the posterior marginal probabilities of models and individual covariates. This allows performing Bayesian model selection and model averaging. For large sample sizes, even the Laplace approximation becomes computationally challenging because the optimisation routine involved needs to evaluate the likelihood on the full set of data in multiple iterations. As a consequence, the algorithm is not scalable for large datasets. To address this problem, we suggest using a version of a popular batch stochastic gradient descent (BSGD) algorithm for estimating the marginal likelihood of a GLM by subsampling from the data. We further combine the algorithm with Markov chain Monte Carlo (MCMC) based methods for Bayesian model selection and provide some theoretical results on the convergence of the estimates. Finally, we report results from experiments illustrating the performance of the proposed algorithm., Comment: 33 pages, 17 figures, tables

We consider the Cauchy problem for a second-order nonlinear evolution equation in a Hilbert space. This equation represents the abstract generalization of the Ball integro-differential equation. The general nonlinear case with respect to terms of the equation which include a square of a norm of a gradient is considered. A three-layer semi-discrete scheme is proposed in order to find an approximate solution. In this scheme, the approximation of nonlinear terms that are dependent on the gradient is carried out by using an integral mean. We show that the solution of the nonlinear discrete problem and its corresponding difference analogue of a first-order derivative is uniformly bounded. For the solution of the corresponding linear discrete problem, it is obtained high-order a priori estimates by using two-variable Chebyshev polynomials. Based on these estimates we prove the stability of the nonlinear discrete problem. For smooth solutions, we provide error estimates for the approximate solution. An iteration method is applied in order to find an approximate solution for each temporal step. The convergence of the iteration process is proved., Comment: 24 pages

We propose new approximate alternating projection methods, based on randomized sketching, for the low-rank nonnegative matrix approximation problem: find a low-rank approximation of a nonnegative matrix that is nonnegative, but whose factors can be arbitrary. We calculate the computational complexities of the proposed methods and evaluate their performance in numerical experiments. The comparison with the known deterministic alternating projection methods shows that the randomized approaches are faster and exhibit similar convergence properties., Comment: Accepted version

This paper aims to interpret the mechanism of feedforward ReLU networks by exploring their solutions for piecewise linear functions, through the deduction from basic rules. The constructed solution should be universal enough to explain some network architectures of engineering; in order for that, several ways are provided to enhance the solution universality. Some of the consequences of our theories include: Under affine-geometry background, the solutions of both three-layer networks and deep-layer networks are given, particularly for those architectures applied in practice, such as multilayer feedforward neural networks and decoders; We give clear and intuitive interpretations of each component of network architectures; The parameter-sharing mechanism for multi-outputs is investigated; We provide an explanation of overparameterization solutions in terms of affine transforms; Under our framework, an advantage of deep layers compared to shallower ones is natural to be obtained. Some intermediate results are the basic knowledge for the modeling or understanding of neural networks, such as the classification of data embedded in a higher-dimensional space, the generalization of affine transforms, the probabilistic model of matrix ranks, and the concept of distinguishable data sets., Comment: v2-v5,v7:hyperlink mode modified or typos corrected; v6:the writing revised; v8:one definition added; v9:corollary 7 refined

In this work, we propose a method for speeding up linear regression distributively, while ensuring security. We leverage randomized sketching techniques, and improve straggler resilience in asynchronous systems. Specifically, we apply a random orthonormal matrix and then subsample in \textit{blocks}, to simultaneously secure the information and reduce the dimension of the regression problem. In our setup, the transformation corresponds to an encoded encryption in an \textit{approximate} gradient coding scheme, and the subsampling corresponds to the responses of the non-straggling workers; in a centralized coded computing network. We focus on the special case of the \textit{Subsampled Randomized Hadamard Transform}, which we generalize to block sampling; and discuss how it can be used to secure the data. We illustrate the performance through numerical experiments., Comment: 3 figures, 5 pages excluding appendices

In this contribution, we introduce a multilevel approximation method with T-splines for fitting scattered point clouds iteratively, with an application to land remote sensing. This new procedure provides a local surface approximation by an explicit computation of the control points and is called a multilevel T-splines approximation (MTA). It is computationally efficient compared with the traditional global least-squares (LS) approach, which may fail when there is an unfavourable point density from a given refinement level. We validate our approach within a simulated framework and apply it to two real datasets: (i) a surface with holes scanned with a terrestrial laser scanner, and (ii) a patch on a sand-dune in the Netherlands. Both examples highlight the potential of the MTA for rapidly fitting large and noisy point clouds with variable point density and with similar results compared to the global LS approximation., Comment: 28 pages, 10 figures

In this paper, we prove a geometrical inequality which states that for any four points on a hemisphere with the unit radius, the largest sum of distances between the points is 4+4*sqrt(2). In our method, we have constructed a rectangular neighborhood of the local maximum point in the feasible set, which size is explicitly determined, and proved that (1): the objective function is bounded by a quadratic polynomial which takes the local maximum point as the unique critical point in the neighborhood, and (2): the rest part of the feasible set can be partitioned into a finite union of a large number of very small cubes so that on each small cube the conjecture can be verified by estimating the objective function with exact numerical computation., Comment: In Proceedings ADG 2021, arXiv:2112.14770

We explore analytically and numerically agglomeration driven by advection and localized source. The system is inhomogeneous in one dimension, viz. along the direction of advection. We analyze a simplified model with mass-independent advection velocity, diffusion coefficient, and reaction rates. We also examine a model with mass-dependent coefficients describing aggregation with sedimentation. For the simplified model, we obtain an exact solution for the stationary spatially dependent agglomerate densities. In the model describing aggregation with sedimentation, we report a new conservation law and develop a scaling theory for the densities. For numerical efficiency we exploit the low-rank approximation technique; this dramatically increases the computational speed and allows simulations of large systems. The numerical results are in excellent agreement with the predictions of our theory., Comment: 10 pages, 6 figures

Controlling nonlinear dynamics is a long-standing problem in engineering. Harnessing known physical information to accelerate or constrain stochastic learning pursues a new paradigm of scientific machine learning. By linearizing nonlinear systems, traditional control methods cannot learn nonlinear features from chaotic data for use in control. Here, we introduce Physics-Informed Deep Operator Control (PIDOC), and by encoding the control signal and initial position into the losses of a physics-informed neural network (PINN), the nonlinear system is forced to exhibit the desired trajectory given the control signal. PIDOC receives signals as physics commands and learns from the chaotic data output from the nonlinear van der Pol system, where the output of the PINN is the control. Applied to a benchmark problem, PIDOC successfully implements control with higher stochasticity for higher-order terms. PIDOC has also been proven to be capable of converging to different desired trajectories based on case studies. Initial positions slightly affect the control accuracy at the beginning stage yet do not change the overall control quality. For highly nonlinear systems, PIDOC is not able to execute control with high accuracy compared with the benchmark problem. The depth and width of the neural network structure do not greatly change the convergence of PIDOC based on case studies of van der Pol systems with low and high nonlinearities. Surprisingly, enlarging the control signal does not help to improve the control quality. The proposed framework can potentially be applied to many nonlinear systems for nonlinear controls., Comment: This framework is not proposed for actual control implementation but for simply simulating the control process using deep learning

Light deviation around a rotating black hole is calculated using the Kerr metric for both small and large deviation angles. For small angles the Lindstedt-Poincar\'e method is employed to get well-behaved solutions, as well as Pad\'e approximants. For large deviation angles numerical integration has been used., Comment: 29 pages, 17 figures

In-Stent Restenosis is a recurrence of coronary artery narrowing due to vascular injury caused by balloon dilation and stent placement. It may lead to the relapse of angina symptoms or to an acute coronary syndrome. An uncertainty quantification of a model for In-Stent Restenosis with four uncertain parameters (endothelium regeneration time, the threshold strain for smooth muscle cells bond breaking, blood flow velocity and the percentage of fenestration in the internal elastic lamina) is presented. Two quantities of interest were studied, namely the average cross-sectional area and the maximum relative area loss in a vessel. Due to the computational intensity of the model and the number of evaluations required in the uncertainty quantification, a surrogate model, based on Gaussian process regression with proper orthogonal decomposition, was developed which subsequently replaced the original In-Stent Restenosis model in the uncertainty quantification. A detailed analysis of the uncertainty propagation and sensitivity analysis is presented. Around 11% and 16% of uncertainty are observed on the average cross-sectional area and maximum relative area loss respectively, and the uncertainty estimates show that a higher fenestration mainly determines uncertainty in the neointimal growth at the initial stage of the process. On the other hand, the uncertainty in blood flow velocity and endothelium regeneration time mainly determine the uncertainty in the quantities of interest at the later, clinically relevant stages of the restenosis process. The uncertainty in the threshold strain is relatively small compared to the other uncertain parameters.

We consider the problem of minimizing a strongly convex function that depends on an uncertain parameter $\theta$. The uncertainty in the objective function means that the optimum, $x^*(\theta)$, is also a function of $\theta$. We propose an efficient method to compute $x^*(\theta)$ and its statistics. We use a chaos expansion of $x^*(\theta)$ along a truncated basis and study first-order methods that compute the optimal coefficients. We establish the convergence rate of the method as the number of basis functions, and hence the dimensionality of the optimization problem is increased. We give the first non-asymptotic rates for the gradient descent and the accelerated gradient descent methods. Our analysis exploits convexity and does not rely on a diminishing step-size strategy. As a result, it is much faster than the state-of-the-art both in theory and in our preliminary numerical experiments. A surprising side-effect of our analysis is that the proposed method also acts as a variance reduction technique to the problem of estimating $x^*(\theta)$., Comment: Minor presentation edits only. Submitted to SIAM UQ Journal

This paper presents a spectral scheme for the numerical solution of nonlinear conservation laws in non-periodic domains under arbitrary boundary conditions. The approach relies on the use of the Fourier Continuation (FC) method for spectral representation of non-periodic functions in conjunction with smooth localized artificial viscosity assignments produced by means of a Shock-Detecting Neural Network (SDNN). Like previous shock capturing schemes and artificial viscosity techniques, the combined FC-SDNN strategy effectively controls spurious oscillations in the proximity of discontinuities. Thanks to its use of a localized but smooth artificial viscosity term, whose support is restricted to a vicinity of flow-discontinuity points, the algorithm enjoys spectral accuracy and low dissipation away from flow discontinuities, and, in such regions, it produces smooth numerical solutions -- as evidenced by an essential absence of spurious oscillations in level set lines. The FC-SDNN viscosity assignment, which does not require use of problem-dependent algorithmic parameters, induces a significantly lower overall dissipation than other methods, including the Fourier-spectral versions of the previous entropy viscosity method. The character of the proposed algorithm is illustrated with a variety of numerical results for the linear advection, Burgers and Euler equations in one and two-dimensional non-periodic spatial domains., Comment: 40 pages, 17 figures, 1 table, 1 algorithm

Non-cleared bilateral OTC derivatives between two financial firms or systemically important non-financial entities are subject to regulations that require the posting of initial and variation margin. The ISDA standard approach (SIMM) provides a way for computing the initial margin. It involves computing sensitivities of the contracts with respect to several market factors. In this paper, the authors extend the well known LSMC technique to efficiently estimate the sensitivities required in the ISDA SIMM methodology., Comment: 14 pages