Iterative Linear Algebra Research Articles

Windowed Green function method for wave scattering by periodic arrays of 2D obstacles

AbstractThis paper introduces a novel boundary integral equation (BIE) method for the numerical solution of problems of planewave scattering by periodic line arrays of two‐dimensional penetrable obstacles. Our approach is built upon a direct BIE formulation that leverages the simplicity of the free‐space Green function but in turn entails evaluation of integrals over the unit‐cell boundaries. Such integrals are here treated via the window Green function method. The windowing approximation together with a finite‐rank operator correction—used to properly impose the Rayleigh radiation condition—yield a robust second‐kind BIE that produces superalgebraically convergent solutions throughout the spectrum, including at the challenging Rayleigh–Wood anomalies. The corrected windowed BIE can be discretized by means of off‐the‐shelf Nyström and boundary element methods, and it leads to linear systems suitable for iterative linear algebra solvers as well as standard fast matrix–vector product algorithms. A variety of numerical examples demonstrate the accuracy and robustness of the proposed methodology.

High-Dimensional Gaussian Sampling: A Review and a Unifying Approach Based on a Stochastic Proximal Point Algorithm

Efficient sampling from a high-dimensional Gaussian distribution is an old but high-stake issue. Vanilla Cholesky samplers imply a computational cost and memory requirements which can rapidly become prohibitive in high dimension. To tackle these issues, multiple methods have been proposed from different communities ranging from iterative numerical linear algebra to Markov chain Monte Carlo (MCMC) approaches. Surprisingly, no complete review and comparison of these methods have been conducted. This paper aims at reviewing all these approaches by pointing out their differences, close relations, benefits and limitations. In addition to this state of the art, this paper proposes a unifying Gaussian simulation framework by deriving a stochastic counterpart of the celebrated proximal point algorithm in optimization. This framework offers a novel and unifying revisit of most of the existing MCMC approaches while extending them. Guidelines to choose the appropriate Gaussian simulation method for a given sampling problem in high dimension are proposed and illustrated with numerical examples.

Variational Density Fitting with a Krylov Subspace Method.

In this work, we present the implementation of a variational density fitting methodology that uses iterative linear algebra for solving the associated system of linear equations. It is well known that most difficulties with this system arise from the fact that the coefficient matrix is in general ill-conditioned and, due to finite precision round-off errors, it may not be positive definite. The dimensionality, given by the number of auxiliary functions, also poses a challenge in terms of memory and time demand since the coefficient matrix is dense. The methodology presented is based on a preconditioned Krylov subspace method able to deal with indefinite ill-conditioned equation systems. To assess its potential, it has been combined with double asymptotic electron repulsion integral expansions as implemented in the deMon2k package. A numerical study on a set of problems with up to 130,000 auxiliary functions shows its effectiveness to alleviate the abovementioned problematic. A comparison with the default methodology used in deMon2k based on a truncated eigenvalue decomposition of the coefficient matrix indicates that the proposed method exhibits excellent robustness and scalability when implemented in a parallel setting.

Implementation of an interior point method with basis preconditioning

The implementation of a linear programming interior point solver is described that is based on iterative linear algebra. The linear systems are preconditioned by a basis matrix, which is updated from one interior point iteration to the next to bound the entries in a certain tableau matrix. The update scheme is based on simplex-type pivot operations and is implemented using linear algebra techniques from the revised simplex method. An initial basis is constructed by a crash procedure after a few interior point iterations. The basis at the end of the interior point solve provides the starting basis for a crossover method which recovers a basic solution to the linear program. Results of a computational study on a diverse set of medium to large-scale problems are discussed.

Quantum Monte Carlo with the Langevin Equation: Coupled Bose-Fermi Systems

We describe the algorithmic details and a performance evaluation of a Langevin approach to a strongly interacting electron-phonon system, and show it has a near linear scaling with lattice size Ns. Many of the limitations of previous attempts to employ such methods to condensed matter lattice Hamiltonians are absent. In particular, the iterative linear algebra solution remains well behaved at strong coupling and low temperatures. The use of Fourier Acceleration is crucial for efficiency, and its use makes the method competitive with the widely-used local update methods, which scale as for on-site interactions and for long range electron-phonon coupling, even on rather small lattice sizes.

Newton projection with proportioning using iterative linear algebra for model predictive control with long prediction horizon

ABSTRACTThis paper presents an algorithm to solve a sparse Quadratic Programming (QP) problem. The QP problem is suitable for Model Predictive Control (MPC) applications in particular. MPC is a modern multivariable control method which requires the solution to a quadratic programming problem at each sampling instant. The proposed algorithm is an active-set based strategy which uses the proportioning test for the selection of the active-set reduction and expansion while utilizing the sparse nature of the problem by the preconditioned MINRES algorithm to solve the face problem. Numerical experiments illustrate the performance of the algorithm, and the results are compared with the state-of-the-art solvers.

A new inexact SQP algorithm for nonlinear systems of mixed equalities and inequalities

Traditional inexact SQP algorithm can only solve equality constrained optimization (Byrd et al. Math. Program. 122, 273–299 2010). In this paper, we propose a new inexact SQP algorithm with affine scaling technique for nonlinear systems of mixed equalities and inequalities, which arise in complementarity and variational inequalities. The nonlinear systems are transformed into a special nonlinear optimization with equality and bound constraints, and then we give a new inexact SQP algorithm for solving it. The new algorithm equipped with affine scaling technique does not require a quadratic programming subproblem with inequality constraints. The search direction is computed by solving one linear system approximately using iterative linear algebra techniques. Under mild assumptions, we discuss the global convergence. The preliminary numerical results show the effectiveness of the proposed algorithm.

Analytical Jacobian-vector products for the matrix-free time integration of partial differential equations

Many scientific and engineering applications require the solution of large systems of initial value problems arising from method of lines discretization of partial differential equations. For systems with widely varying time scales, or with complex physical dynamics, implicit time integration schemes are preferred due to their superior stability properties. These schemes solve at each step linear systems with matrices formed using the Jacobian of the right hand side function. For large applications iterative linear algebra methods, which make use of Jacobian-vector products, are employed. This paper studies the impact that the method of computing Jacobian-vector products has on the overall performance and accuracy of the time integration process. The analysis shows that the most beneficial approach is the direct computation of exact Jacobian-vector products in the context of matrix-free time integrators. This approach does not suffer from approximation errors, reuses the parallelism and data distribution already present in the right-hand side vector computations, and avoids storing or operating on the entire Jacobian matrix.

RT-CUDA: A Software Tool for CUDA Code Restructuring

Recent development in graphic processing units (GPUs) has opened a new challenge in harnessing their computing power as a new general purpose computing paradigm. However, porting applications to CUDA remains a challenge to average programmers, which have to package code in separate functions, explicitly manage data transfers between the host and device memories, and manually optimize GPU memory utilization. In this paper, we propose a restructuring tool (RT-CUDA) that takes a C-like program and some user directives as compiler hints to produce an optimized CUDA code. The tool strategy is based on efficient management of the memory system to minimize data motion by managing the transfer between host and device, maximizing bandwidth for device memory accesses, and enhancing data locality and re-use of cached data using shared-memory and registers. Enhanced resource utilization is implemented by re-writing code as parametric kernels and use of efficient auto-tuning. The tool enables calling numerical libraries (CuBLAS, CuSPARSE, etc.) to help implement applications in science simulation like iterative linear algebra solvers. For the above applications, the tool implement an inter-block global synchronization which allow the execution overall among a few iterations which is helpful to balance load and to avoid polling. Evaluation of RT-CUDA has been performed using a variety of basic linear algebra operators (Madd, MM, MV, VV, etc.) as well as the programming of iterative solvers for systems of linear equations like Jacobi and Conjugate Gradient algorithms. Significant speedup has been achieved over other compilers like PGI OpenACC and GPGPU compilers for the above applications. Evaluation shows that generated kernels efficiently call math libraries and enable implementing complete iterative solvers. The tool help scientists developing parallel simulators like reservoir simulators, molecular dynamics, etc. without exposing to complexity of GPU and CUDA programming. We have partnership with a group of researchers at the Saudi Aramco, a national company in Saudi Arabia. RT-CUDA is currently explored as a potential development tool for applications involving linear algebra solvers by the above group. In addition, RT-CUDA is being used by Senior and Graduate students at King Fahd University of Petroleum and Minerals in their projects as part of RT-CUDA continuous enhancement.

Verification of continuum drift kinetic equation solvers in NIMROD

Verification of continuum solutions to the electron and ion drift kinetic equations (DKEs) in NIMROD [C. R. Sovinec et al., J. Comp. Phys. 195, 355 (2004)] is demonstrated through comparison with several neoclassical transport codes, most notably NEO [E. A. Belli and J. Candy, Plasma Phys. Controlled Fusion 54, 015015 (2012)]. The DKE solutions use NIMROD's spatial representation, 2D finite-elements in the poloidal plane and a 1D Fourier expansion in toroidal angle. For 2D velocity space, a novel 1D expansion in finite elements is applied for the pitch angle dependence and a collocation grid is used for the normalized speed coordinate. The full, linearized Coulomb collision operator is kept and shown to be important for obtaining quantitative results. Bootstrap currents, parallel ion flows, and radial particle and heat fluxes show quantitative agreement between NIMROD and NEO for a variety of tokamak equilibria. In addition, velocity space distribution function contours for ions and electrons show nearly identical detailed structure and agree quantitatively. A Θ-centered, implicit time discretization and a block-preconditioned, iterative linear algebra solver provide efficient electron and ion DKE solutions that ultimately will be used to obtain closures for NIMROD's evolving fluid model.

Scalable Nonlinear Programming via Exact Differentiable Penalty Functions and Trust-Region Newton Methods

We present an approach for nonlinear programming based on the direct minimization of an exact differentiable penalty function using trust-region Newton techniques. The approach provides desirable features required for scalability: it can detect and exploit directions of negative curvature, it is superlinearly convergent, and it enables the scalable computation of the Newton step through iterative linear algebra. Moreover, it presents features that are desirable for parametric optimization problems that must be solved in a latency-limited environment, as is the case for model predictive control and mixed-integer nonlinear programming. These features are fast detection of activity, efficient warm starting, and progress on a primal-dual merit function at every iteration. We note that other algorithmic approaches fail to satisfy at least one of these features. We derive general convergence results for our approach and demonstrate its behavior through numerical studies.

A dwindling filter inexact projected Hessian algorithm for large scale nonlinear constrained optimization

In this paper, we propose a dwindling filter inexact projected Hessian algorithm for solving large scale nonlinear constrained optimization. For large-scale applications, inexact projected Hessian algorithm is needed to get search direction by solving one or more linear systems approximately using iterative linear algebra techniques. The envelope of the dwindling filter becomes thinner and thinner as the step size approaches zero so that the new filter has more flexibility for the acceptance of the trial step compared with traditional filter. Under mild conditions, global convergence and local superlinear convergence rate are obtained. The numerical experiments are reported to show the effectiveness of the proposed algorithm for large scale problems.

Education

Take a random polygon, find the the midpoints of its sides, and construct a new polygon from them. Repeat ad infinitum. What happens? This simple geometric iteration is at the heart of this issue's Education section offering, “From Random Polygon to Ellipse: An Eigenanalysis” by Adam N. Elmachtoub and Charles F. Van Loan. A cursory numerical experiment (although the traditionalist may wish to try this with a compass and straightedge) shows that the polygon shrinks to a point which the astute observer will realize is the centroid of the original vertices. If we now zoom in on the sequence of shrinking polygons, most easily accomplished by rescaling at each iteration, we see even more order. The rescaled vertices approach a binary oscillation between two polygons, both of which can be inscribed in the same ellipse. While students posed with these questions may think they are pursuing a geometry problem, the not-too-secret agenda is to introduce them to some fundamental ideas of numerical linear algebra, namely iteration, rescaling, and convergence. The authors guide us to write the problem as an iterative algorithm whose properties can be analyzed via eigenvalues and eigenvectors and whose convergence is akin to the power method for finding the eigenvalue of largest magnitude of a matrix. One might consider threading this example through a course in numerical analysis or numerical linear algebra, having the students explore the behavior at the beginning of the course (a great programming exercise), and eventually leading them through the analysis after introducing iterative methods for computing eigenvalues and eigenvectors. The problem opens the door to many generalizations, including exploring what happens for different geometric constructions and asking for geometric interpretations for other iterative linear algebra algorithms. The midpoint and rescaling construction is also simple enough to show to freshmen or even high school students; who knows, perhaps a few of them may be drawn to the study of mathematics by the allure of binary blinking ellipses?

A parallel computational model for sensitivity analysis in optimization for robustness

We present an efficient method for computing, in a parallel distributed environment, gradients of functionals depending on the solution sensitivities of dynamical systems described by ordinary differential equations. This work was motivated by the need to compute cost function gradients for dynamically constrained optimization for robustness, i.e. the problem of finding the model parameters of a given dynamical system that lead to minimum solution sensitivity. The proposed approach for this particular second-order sensi-tivity problem falls into the class of the so-called ‘adjoint-over-forward’ methods and is based on solving the continuous forward sensitivity equations (for evaluating the cost functional) and their continuous adjoint systems (for evaluation of the gradient of the cost functional). Efficiency is ensured by: (i) a suitable distribution over processors of the individual sensitivity systems, both forward and adjoint so that inter-process communication is minimized; and (ii) use of an iterative solver combined with a block-diagonal preconditioner for the solution of the linear systems arising in the implicit integration of the resulting ordinary differential equations (ODE) systems. The proposed algorithm was implemented as an extension of the cvodes solver in sundials [P. Brown, K. Grant, A. Hindmarsh, S. Lee, R.S.D. Shumaker, and C. Woodward. SUNDIALS: SUite of Nonlinear and DIfferential/ALgebraic equation Solvers, ACM Trans. Math. Software, 31(3) (2005), pp. 363–396], but it can be used in conjunction with any sensitivity-enabled ODE integrator that provides adjoint sensitivity capabilities and (if based on implicit integration methods) support for iterative linear algebra.

An Inexact SQP Method for Equality Constrained Optimization

We present an algorithm for large-scale equality constrained optimization. The method is based on a characterization of inexact sequential quadratic programming (SQP) steps that can ensure global convergence. Inexact SQP methods are needed for large-scale applications for which the iteration matrix cannot be explicitly formed or factored and the arising linear systems must be solved using iterative linear algebra techniques. We address how to determine when a given inexact step makes sufficient progress toward a solution of the nonlinear program, as measured by an exact penalty function. The method is globalized by a line search. An analysis of the global convergence properties of the algorithm and numerical results are presented.

Toward the reusability for iterative linear algebra software in distributed environment

LAKe (Linear Algebra Kernel) is a linear algebra class library developed using an object oriented approach in order to enable a good reuse code for sequential or parallel linear algebra applications. An application implemented using LAKe can be run either in sequential or in parallel mode using the same code. This paper proposes an extension to the LAKe library, which enables an implemented method to run in sequential and parallel mode simultaneously. This extension of LAKe allows a concurrent reuse of sequential and parallel components inside the same linear algebra application. To take advantage of such concurrent reusability, the sequential and parallel components have to be able to collaborate with each other. The hybrid methods need such properties to run. These methods are defined by a combination of several numerical techniques, or several copies of the same method parameterized differently in order to accelerate the convergence and/or to improve the accuracy of the solution of a given large linear algebra problem. In order to validate our approach, we present some numerical experiments, making use of a hybrid method with the proposed extension of LAKe.

Systolic array implementation of some iterative algorithms for solving systems of linear algebraic equations

The paper elucidates the internal properties of some iterative linear algebra algorithms through transformation of their lattice graphs.