Linear System Solver Research Articles

Optimal power flow solution via noise-resilient quantum interior-point methods

This paper presents quantum interior-point methods (QIPMs) tailored to tackle the DC optimal power flow (OPF) problem using noisy intermediate-scale quantum devices. The optimization model is redefined as a linearly constrained quadratic optimization. By incorporating the Harrow–Hassidim–Lloyd (HHL) quantum algorithm into the IPM framework, Newton’s direction is determined through the resolution of linear equation systems. To mitigate the impact of HHL error and quantum noise on Newton’s direction, we utilized a noise-tolerant quantum IPM. This approach provides high-quality OPF solutions even in scenarios where inexact solutions to the linear equation systems result in approximated Newton’s direction. To enhance performance in cases of slow convergence and uphold the feasibility of OPF upon convergence, we propose a classically augmented noise-tolerant QIPM. This technique is designed to expedite convergence relative to classical IPM while maintaining the accuracy of the solution. The proposed QIPM variants are studied through comprehensive simulations and error analyses on 3-bus, 5-bus, 118-bus, and 300-bus systems. By modeling the errors and incorporating quantum computer noise, we simulate the algorithms on Qiskit and classical computers to better understand their effectiveness and feasibility under realistic conditions.

Unconditionally convergent [formula omitted] splitting iterative methods for variable coefficient Riesz space fractional diffusion equations

In this paper, we consider fast solvers for discrete linear systems generated by Riesz space fractional diffusion equations. We extract a scalar matrix, a compensation matrix, and a τ matrix from the coefficient matrix, and use their sum to construct a class of τ splitting iterative methods. Additionally, we design a preconditioner for the conjugate gradient method. Theoretical analyses show that the proposed τ splitting iterative methods are unconditionally convergent with convergence rates independent of step-sizes. Numerical results are provided to demonstrate the effectiveness of the proposed iterative methods.

A Note on the Convergence of Multigrid Methods for the Riesz–Space Equation and an Application to Image Deblurring

In recent decades, a remarkable amount of research has been carried out regarding fast solvers for large linear systems resulting from various discretizations of fractional differential equations (FDEs). In the current work, we focus on multigrid methods for a Riesz–Space FDE whose theoretical convergence analysis of such multigrid methods is currently limited in the relevant literature to the two-grid method. Here we provide a detailed theoretical convergence study in the multilevel setting. Moreover, we discuss its use combined with a band approximation and we compare the result with both τ and circulant preconditionings. The numerical tests include 2D problems as well as the extension to the case of a Riesz–FDE with variable coefficients. Finally, we investigate the use of a Riesz–Space FDE in a variational model for image deblurring, comparing the performance of specific preconditioning strategies.

A nonsmooth primal-dual method with interwoven PDE constraint solver

We introduce an efficient first-order primal-dual method for the solution of nonsmooth PDE-constrained optimization problems. We achieve this efficiency through not solving the PDE or its linearisation on each iteration of the optimization method. Instead, we run the method interwoven with a simple conventional linear system solver (Jacobi, Gauss–Seidel, conjugate gradients), always taking only one step of the linear system solver for each step of the optimization method. The control parameter is updated on each iteration as determined by the optimization method. We prove linear convergence under a second-order growth condition, and numerically demonstrate the performance on a variety of PDEs related to inverse problems involving boundary measurements.

Recycling Krylov subspaces for efficient partitioned solution of aerostructural adjoint systems

Robust and efficient solvers for coupled-adjoint linear systems are crucial to successful aerostructural optimization. Monolithic and partitioned strategies can be applied. The monolithic approach is expected to offer better robustness and efficiency for strong fluid-structure interactions. However, it requires a higher implementation cost and convergence may depend on appropriate scaling and initialization strategies. On the other hand, the modularity of the partitioned method enables a straightforward implementation while its convergence may require relaxation. In addition, a partitioned solver often leads to a higher number of iterations to get the same level of convergence as the monolithic one.The objective of this paper is to accelerate the fluid-structure coupled-adjoint partitioned solver by considering techniques borrowed from approximate invariant subspace recycling strategies adapted to sequences of linear systems with varying right-hand sides. Indeed, in a partitioned framework, the structural source term attached to the fluid block of equations affects the right-hand side with the nice property of quickly converging to a constant value. For the fluid block of equations, we also consider deflation of approximate eigenvectors in conjunction with advanced inner-outer Krylov solvers.To demonstrate the effectiveness of our proposed recycling strategy, we compute the coupled derivatives for two emblematic problems in transonic viscous flow. The first one is an aeroelastic configuration of the ONERA-M6 fixed wing. For this exercise, the fluid grid was coupled to a structural model specifically designed to exhibit a high flexibility. The second one is a wing-body aeroelastic configuration of the Common Research Model. All computations are performed using a fully linearized one-equation Spalart-Allmaras turbulence model. Numerical simulations show up to 39% reduction in matrix-vector products for GCRO-DR and up to 25% for the nested FGCRO-DR solver.

An Optimal Linear-combination-of-unitaries-based Quantum Linear System Solver

Solving systems of linear equations is one of the most important primitives in many different areas, including in optimization, simulation, and machine learning. Quantum algorithms for solving linear systems have the potential to provide a quantum advantage for these problems. In this work, we recall the Chebyshev iterative method and the corresponding optimal polynomial approximation of the inverse. We show that the Chebyshev iteration polynomial can be efficiently evaluated both using quantum singular value transformation (QSVT) as well as linear combination of unitaries (LCU). We achieve this by bounding the 1-norm of the coefficients of the polynomial expressed in the Chebyshev basis. This leads to a considerable constant-factor improvement in the runtime of quantum linear system solvers that are based on LCU or QSVT (or, conversely, a several orders of magnitude smaller error with the same runtime/circuit depth).

Equivariant neural operators for gradient-consistent topology optimization

Abstract Most traditional methods for solving partial differential equations (PDEs) require the costly solving of large linear systems. Neural operators (NOs) offer remarkable speed-ups over classical numerical PDE solvers. Here, we conduct the first exploration and comparison of NOs for three-dimensional topology optimization. Specifically, we propose replacing the PDE solver within the popular Solid Isotropic Material with Penalization (SIMP) algorithm, which is its main computational bottleneck. For this, the NO not only needs to solve the PDE with sufficient accuracy but also has the additional challenge of providing accurate gradients which are necessary for SIMP’s density updates. To realize this, we do three things: (i) We introduce a novel loss term to promote gradient-consistency. (ii) We guarantee equivariance in our NOs to increase the physical correctness of predictions. (iii) We introduce a novel NO architecture called U-Net Fourier neural operator (U-Net FNO), which combines the multi-resolution properties of U-Nets with the Fourier neural operator (FNO)’s focus on local features in frequency space. In our experiments we demonstrate that the inclusion of the novel gradient loss term is necessary to obtain good results. Furthermore, enforcing group equivariance greatly improves the quality of predictions, especially on small training datasets. Finally, we show that in our experiments the U-Net FNO outperforms both a standard U-Net, as well as other FNO methods.

High-Throughput Hardware Design for Linear Equation System Solving of VVC Affine Prediction

The affine prediction is the main novelty in the inter-frame prediction of Versatile Video Coding (VVC) standard. However, implementing the affine prediction requires a huge computational effort that makes mandatory the use of dedicated hardware accelerators to achieve real-time processing and meet the constraints of area and power dissipation for mobile devices. In light of this, this work presents different approaches for a hardware design dedicated to solving the linear equation system, which is essential for refiningthe affine Motion Vector. Synthesis results show that the High- Throughput architecture, which explores the parallelism of internal operations, can reach an accuracy of 99.99%, requiring 333.9kgates to be implemented and presenting a power dissipation of 120.4mW when running at 540MHz, the operational frequency required to process 60 frames per second of HD 1080p videos.

A study of concurrent multi-frontal solvers for modern massively parallel architectures

Leveraging Trace Theory, we investigate the efficient parallelization of direct solvers for large linear equation systems. Our focus lies on a multi-frontal algorithm, and we present a methodology for achieving near-optimal scheduling on modern massively parallel machines. By employing trace theory with Diekert Graphs and Foata Normal Form, we rigorously validate the correctness of our proposed solution. To establish a strong link between the mesh and elimination tree of the multi-frontal solver, we conduct extensive testing on matrices derived from the Finite Element Method (FEM). Furthermore, we assess the performance of computations on both GPU and CPU platforms, employing practical implementation strategies.

A comparison of numerical schemes for the GPU-accelerated simulation of variably-saturated groundwater flow

The optimal strategy to numerically solve the Richards equation is problem-dependent, especially on GPUs because an efficient numerical scheme on CPU might not scale well on GPU. In this work, four numerical schemes are coded to investigate their parallel performance on CPU and GPU. The results indicate that the scaling of Richards solvers on GPU is affected by the numerical scheme, the linear system solver, the soil constitutive model, the code structure, the problem size and the adaptive time stepping strategies. Compared with CPU, parallel simulations on GPU exhibit stronger variance in the scaling of different code sections. The poorly-scaled components could deteriorate the overall scaling. Under all circumstances, using a GPU significantly enhances computational speed, especially for large problems. Clearly, GPU computing have significant potential in accelerating large-scale hydrological simulations, but care must be taken on the design and implementation of the model structure.

Second-order Rosenbrock-exponential (ROSEXP) methods for partitioned differential equations

In this paper, we introduce a new framework for deriving partitioned implicit-exponential integrators for stiff systems of ordinary differential equations and construct several time integrators of this type. The new approach is suited for solving systems of equations where the forcing term is comprised of several additive nonlinear terms. We analyze the stability, convergence, and efficiency of the new integrators and compare their performance with existing schemes for such systems using several numerical examples. We also propose a novel approach to visualizing the linear stability of the partitioned schemes, which provides a more intuitive way to understand and compare the stability properties of various schemes. Our new integrators are A-stable, second-order methods that require only one call to the linear system solver and one exponential-like matrix function evaluation per time step.

AI‐enhanced iterative solvers for accelerating the solution of large‐scale parametrized systems

SummaryRecent advances in the field of machine learning open a new era in high performance computing for challenging computational science and engineering applications. In this framework, the use of advanced machine learning algorithms for the development of accurate and cost‐efficient surrogate models of complex physical processes has already attracted major attention from scientists. However, despite their powerful approximation capabilities, surrogate model predictions are still far from being near to the ‘exact’ solution of the problem. To address this issue, the present work proposes the use of up‐to‐date machine learning tools in order to equip a new generation of iterative solvers of linear equation systems, capable of very efficiently solving large‐scale parametrized problems at any desired level of accuracy. The proposed approach consists of the following two steps. At first, a reduced set of model evaluations is performed using a standard finite element methodology and the corresponding solutions are used to establish an approximate mapping from the problem's parametric space to its solution space using a combination of deep feedforward neural networks and convolutional autoencoders. This mapping serves as a means of obtaining very accurate initial predictions of the system's response to new query points at negligible computational cost. Subsequently, an iterative solver inspired by the Algebraic Multigrid method in combination with Proper Orthogonal Decomposition, termed POD‐2G, is developed that successively refines the initial predictions of the surrogate model towards the exact solution. The application of POD‐2G as a standalone solver or as preconditioner in the context of preconditioned conjugate gradient methods is demonstrated on several numerical examples of large scale systems, with the results indicating its strong superiority over conventional iterative solution schemes.

Sparse Approximate Multifrontal Factorization with Composite Compression Methods

This article presents a fast and approximate multifrontal solver for large sparse linear systems. In a recent work by Liu et al., we showed the efficiency of a multifrontal solver leveraging the butterfly algorithm and its hierarchical matrix extension, HODBF (hierarchical off-diagonal butterfly) compression to compress large frontal matrices. The resulting multifrontal solver can attain quasi-linear computation and memory complexity when applied to sparse linear systems arising from spatial discretization of high-frequency wave equations. To further reduce the overall number of operations and especially the factorization memory usage to scale to larger problem sizes, in this article we develop a composite multifrontal solver that employs the HODBF format for large-sized fronts, a reduced-memory version of the nonhierarchical block low-rank format for medium-sized fronts, and a lossy compression format for small-sized fronts. This allows us to solve sparse linear systems of dimension up to 2.7 × larger than before and leads to a memory consumption that is reduced by 70% while ensuring the same execution time. The code is made publicly available in GitHub.

Randomized Block Adaptive Linear System Solvers

Randomized Block Adaptive Linear System Solvers

Resolution of Linear Systems Using Interval Arithmetic and Cholesky Decomposition

Resolution of Linear Systems Using Interval Arithmetic and Cholesky Decomposition

On Krylov methods for large-scale CBCT reconstruction.

Krylov subspace methods are a powerful family of iterative solvers for linear systems of equations, which are commonly used for inverse problems due to their intrinsic regularization properties. Moreover, these methods are naturally suited to solve large-scale problems, as they only require matrix-vector products with the system matrix (and its adjoint) to compute approximate solutions, and they display a very fast convergence. Even if this class of methods has been widely researched and studied in the numerical linear algebra community, its use in applied medical physics and applied engineering is still very limited. e.g. in realistic large-scale computed tomography (CT) problems, and more specifically in cone beam CT (CBCT). This work attempts to breach this gap by providing a general framework for the most relevant Krylov subspace methods applied to 3D CT problems, including the most well-known Krylov solvers for non-square systems (CGLS, LSQR, LSMR), possibly in combination with Tikhonov regularization, and methods that incorporate total variation regularization. This is provided within an open source framework: the tomographic iterative GPU-based reconstruction toolbox, with the idea of promoting accessibility and reproducibility of the results for the algorithms presented. Finally, numerical results in synthetic and real-world 3D CT applications (medical CBCT andμ-CT datasets) are provided to showcase and compare the different Krylov subspace methods presented in the paper, as well as their suitability for different kinds of problems.

An Adaptive Two-Grid Solver for DPG Formulation of Compressible Navier–Stokes Equations in 3D

Abstract We present an overlapping domain decomposition iterative solver for linear systems resulting from the discretization of compressible viscous flows with the Discontinuous Petrov–Galerkin (DPG) method in three dimensions. It is a two-grid solver utilizing the solution on the auxiliary coarse grid and the standard block-Jacobi iteration on patches of elements defined by supports of the coarse mesh base shape functions. The simple iteration defined in this way is used as a preconditioner for the conjugate gradient procedure. Theoretical analysis indicates that the condition number of the preconditioned system should be independent of the actual finite element mesh and the auxiliary coarse mesh, provided that they are quasiuniform. Numerical tests confirm this result. Moreover, they show that presence of strongly flattened or elongated elements does not slow the convergence. The finite element mesh is subject to adaptivity, i.e. dividing the elements with large errors until a required accuracy is reached. The auxiliary coarse mesh is adjusting to the nonuniform actual mesh.

Finite element based direct and iterative approach to investigate a magneto-micropolar flow through a rectangular channel

This article aims to numerically compute the magnetic induced flow in the framework of higher order continuum through its implementation using FreeFem++. The flow through a rectangular channel is assumed laminar, incompressible and under constant applied pressure. External magnetic field is applied through the boundary of the channel. To compute the physical response theory of micropolar continuum is taken into consideration and governing dynamics is described in the form of coupled partial differential equations in velocity, induction and microrotation field variables along with associated boundary conditions. Sparse linear system is thus obtained with the application of finite element method. The resulting system is then solved by sparse solver available in unsymmetric multifrontal package (UMFPACK) through FreeFem++. The variational form of the model is presented and implemented in FreeFem++. Pertinent to physical significance different parameters i.e., magnetic Reynolds number, Hartmann number, micropolar coupling constant, micropolar parameter are studied and results are shown through simulations and graphs.It is observed that with an increase in micropolar coupling parameter the magnetic induction as well as the microrotational velocity of particles decreases in the medium. The maximum magnitude of microrotational velocities is observed near walls of the channel. The microrotational spin in the channel flow slows down due to an increase in induced magnetic field. Moreover, it is found that direct solver for the asymmetric sparse linear system arise in this case does not lead to stable solution. This instability in the solution arises when the Hartmann number in simulations increases than a value of 1.5. Therefore, iterative solver based on generalized minimal residual method (GMRES) should be adopted to weed out the instabilities due to fine scale oscillations of translational as well as microrotational motions of the fluid particles in case of large Hartmann number.

PLSS: A Projected Linear Systems Solver

PLSS: A Projected Linear Systems Solver

Acceleration of Nuclear Reactor Simulation and Uncertainty Quantification Using Low-Precision Arithmetic

In recent years, interest in approximate computing has been increasing significantly in many disciplines in the context of saving energy and computation cost by trading off on the quality of numerical simulation. The hardware acceleration based on the low-precision floating-point arithmetic is anticipated by the upcoming generation of microprocessors and code compilers and has already proven to be beneficial for weather and climate modelling and neural network training. The present work illustrates the application of low-precision arithmetic for the nuclear reactor core uncertainty analysis. We studied the performance of an elementary transient reactor core model for the arbitrary precision of the floating-point multiplication in a direct linear system solver. Using this model, we calculated the reactor core transients initiated by the control rod ejection taking into account the uncertainty of the model input parameters. Then, we evaluated the round-off errors of the model outputs for different precision levels. The comparison of the round-off errors and the model uncertainty showed the model could be run using a 15-bit floating-point number precision with an acceptable degradation of the result’s accuracy. This precision corresponds to a gain of about 6× in the bit complexity of the linear system solution algorithm, which can be actualized in terms of reduced energy costs on low-precision hardware.