We introduce a new class of preconditioners to enable flexible GMRES to find a least-squares solution, and potentially the pseudoinverse solution, of large-scale sparse, asymmetric, singular, and potentially inconsistent systems. We develop the preconditioners based on a new observation that generalized inverses (i.e., $\boldsymbol{A}^{g}\in\{\boldsymbol{G}\mid\boldsymbol{A}\boldsymbol{G}\boldsymbol{A}=\boldsymbol{A}\}$) enable the preconditioned Krylov subspaces (KSP) to converge in a single step. We then compute an approximate generalized inverse (AGI) efficiently using a hybrid incomplete factorization (HIF), which combines multilevel incomplete LU with rank-revealing QR on its final Schur complement. We define the criteria of $\epsilon$-accuracy and stability of AGI to guarantee the convergence of preconditioned GMRES for consistent systems. For inconsistent systems, we fortify HIF with iterative refinement to obtain HIFIR, which allows accurate computations of the null-space vectors. By combining the two techniques, we then obtain a new solver, called PIPIT, for obtaining the pseudoinverse solutions for systems with low-dimensional null spaces. We demonstrate the robustness of HIF and HIFIR and show that they improve both accuracy and efficiency of the prior state of the art by orders of magnitude for systems with up to a million unknowns.
The binary $k$-dimensional simplex code is known to be a $2^{k-1}$-batch code and is conjectured to be a $2^{k-1}$-functional batch code. Here, we offer a simple, constructive proof of a result that is "in between" these two properties. Our approach is to relate these properties to certain (old and new) additive problems in finite abelian groups. We also formulate a conjecture for finite abelian groups that generalizes the above-mentioned conjecture.
This paper describes and compares some structure preserving techniques for the solution of linear discrete ill-posed problems with the t-product. A new randomized tensor singular value decomposition (R-tSVD) with a t-product is presented for low tubal rank tensor approximations. Regularization of linear inverse problems by truncated tensor eigenvalue decomposition (T-tEVD), truncated tSVD (T-tSVD), randomized T-tSVD (RT-tSVD), t-product Golub-Kahan bidiagonalization (tGKB) process, and t-product Lanczos (t-Lanczos) process are considered. A solution method that is based on reusing tensor Krylov subspaces generated by the tGKB process is described. The regularization parameter is the number of iterations required by each method. The discrepancy principle is used to determine this parameter. Solution methods that are based on truncated iterations are compared with solution methods that combine Tikhonov regularization with the tGKB and t-Lanczos processes. Computed examples illustrate the performance of these methods when applied to image and gray-scale video restorations. Our new RT-tSVD method is seen to require less CPU time and yields restorations of higher quality than the T-tSVD method.
The Poisson pressure solve resulting from the spectral element discretization of the incompressible Navier-Stokes equation requires fast, robust, and scalable preconditioning. In the current work, a parallel scaling study of Chebyshev-accelerated Schwarz and Jacobi preconditioning schemes is presented, with special focus on GPU architectures, such as OLCF's Summit. Convergence properties of the Chebyshev-accelerated schemes are compared with alternative methods, such as low-order preconditioners combined with algebraic multigrid. Performance and scalability results are presented for a variety of preconditioner and solver settings. The authors demonstrate that Chebyshev-accelerated-Schwarz methods provide a robust and effective smoothing strategy when using $p$-multigrid as a preconditioner in a Krylov-subspace projector. At the same time, optimal preconditioning parameters can vary for different geometries, problem sizes, and processor counts. This variance motivates the development of an autotuner to optimize solver parameters on-line, during the course of production simulations.
The discretization of robust quadratic optimal control problems under uncertainty using the finite element method and the stochastic collocation method leads to large saddle-point systems, which are fully coupled across the random realizations. Despite its relevance for numerous engineering problems, the solution of such systems is notoriusly challenging. In this manuscript, we study efficient preconditioners for all-at-once approaches using both an algebraic and an operator preconditioning framework. We show in particular that for values of the regularization parameter not too small, the saddle-point system can be efficiently solved by preconditioning in parallel all the state and adjoint equations. For small values of the regularization parameter, robustness can be recovered by the additional solution of a small linear system, which however couples all realizations. A mean approximation and a Chebyshev semi-iterative method are investigated to solve this reduced system. Our analysis considers a random elliptic partial differential equation whose diffusion coefficient $\kappa(x,\omega)$ is modeled as an almost surely continuous and positive random field, though not necessarily uniformly bounded and coercive. We further provide estimates on the dependence of the preconditioned system on the variance of the random field. Such estimates involve either the first or second moment of the random variables $1/\min_{x\in \overline{D}} \kappa(x,\omega)$ and $\max_{x\in \overline{D}}\kappa(x,\omega)$, where $D$ is the spatial domain. The theoretical results are confirmed by numerical experiments, and implementation details are further addressed.
A method is presented for the evaluation of integrals on tetrahedra where the integrand has a singularity at one vertex. The approach uses a transformation to spherical polar coordinates which explicitly eliminates the singularity and facilitates the evaluation of integration limits. The method can also be implemented in an adaptive form which gives convergence to a required tolerance. Results from the method are compared to the output from an exact analytical method and show high accuracy. In particular, when the adaptive algorithm is used, highly accurate results are found for poorly conditioned tetrahedra which normally present difficulties for numerical quadrature techniques.
Potts models, which can be used to analyze dependent observations on a lattice, have seen widespread application in a variety of areas, including statistical mechanics, neuroscience, and quantum computing. To address the intractability of Potts likelihoods for large spatial fields, we propose fast ordered conditional approximations that enable rapid inference for observed and hidden Potts models. Our methods can be used to directly obtain samples from the approximate joint distribution of an entire Potts field. The computational complexity of our approximation methods is linear in the number of spatial locations; in addition, some of the necessary computations are naturally parallel. We illustrate the advantages of our approach using simulated data and a satellite image.
Much recent interest has focused on the design of optimization algorithms from the discretization of an associated optimization flow, i.e., a system of differential equations (ODEs) whose trajectories solve an associated optimization problem. Such a design approach poses an important problem: how to find a principled methodology to design and discretize appropriate ODEs. This paper aims to provide a solution to this problem through the use of contraction theory. We first introduce general mathematical results that explain how contraction theory guarantees the stability of the implicit and explicit Euler integration methods. Then, we propose a novel system of ODEs, namely the Accelerated-Contracting-Nesterov flow, and use contraction theory to establish it is an optimization flow with exponential convergence rate, from which the linear convergence rate of its associated optimization algorithm is immediately established. Remarkably, a simple explicit Euler discretization of this flow corresponds to the Nesterov acceleration method. Finally, we present how our approach leads to performance guarantees in the design of optimization algorithms for time-varying optimization problems.
Quantum many-body systems involving bosonic modes or gauge fields have infinite-dimensional local Hilbert spaces which must be truncated to perform simulations of real-time dynamics on classical or quantum computers. To analyze the truncation error, we develop methods for bounding the rate of growth of local quantum numbers such as the occupation number of a mode at a lattice site, or the electric field at a lattice link. Our approach applies to various models of bosons interacting with spins or fermions, and also to both abelian and non-abelian gauge theories. We show that if states in these models are truncated by imposing an upper limit $\Lambda$ on each local quantum number, and if the initial state has low local quantum numbers, then an error at most $\epsilon$ can be achieved by choosing $\Lambda$ to scale polylogarithmically with $\epsilon^{-1}$, an exponential improvement over previous bounds based on energy conservation. For the Hubbard-Holstein model, we numerically compute a bound on $\Lambda$ that achieves accuracy $\epsilon$, obtaining significantly improved estimates in various parameter regimes. We also establish a criterion for truncating the Hamiltonian with a provable guarantee on the accuracy of time evolution. Building on that result, we formulate quantum algorithms for dynamical simulation of lattice gauge theories and of models with bosonic modes; the gate complexity depends almost linearly on spacetime volume in the former case, and almost quadratically on time in the latter case. We establish a lower bound showing that there are systems involving bosons for which this quadratic scaling with time cannot be improved. By applying our result on the truncation error in time evolution, we also prove that spectrally isolated energy eigenstates can be approximated with accuracy $\epsilon$ by truncating local quantum numbers at $\Lambda=\textrm{polylog}(\epsilon^{-1})$.
This paper addresses the problem of formally verifying desirable properties of neural networks, i.e., obtaining provable guarantees that neural networks satisfy specifications relating their inputs and outputs (robustness to bounded norm adversarial perturbations, for example). Most previous work on this topic was limited in its applicability by the size of the network, network architecture and the complexity of properties to be verified. In contrast, our framework applies to a general class of activation functions and specifications on neural network inputs and outputs. We formulate verification as an optimization problem (seeking to find the largest violation of the specification) and solve a Lagrangian relaxation of the optimization problem to obtain an upper bound on the worst case violation of the specification being verified. Our approach is anytime i.e. it can be stopped at any time and a valid bound on the maximum violation can be obtained. We develop specialized verification algorithms with provable tightness guarantees under special assumptions and demonstrate the practical significance of our general verification approach on a variety of verification tasks.
Network embedding has attracted considerable research attention recently. However, the existing methods are incapable of handling billion-scale networks, because they are computationally expensive and, at the same time, difficult to be accelerated by distributed computing schemes. To address these problems, we propose RandNE, a novel and simple billion-scale network embedding method. Specifically, we propose a Gaussian random projection approach to map the network into a low-dimensional embedding space while preserving the high-order proximities between nodes. To reduce the time complexity, we design an iterative projection procedure to avoid the explicit calculation of the high-order proximities. Theoretical analysis shows that our method is extremely efficient, and friendly to distributed computing schemes without any communication cost in the calculation. We demonstrate the efficacy of RandNE over state-of-the-art methods in network reconstruction and link prediction tasks on multiple datasets with different scales, ranging from thousands to billions of nodes and edges.