Large Systems Of Linear Equations Research Articles

A global Lanczos method for image restoration

Image restoration often requires the solution of large linear systems of equations with a very ill-conditioned, possibly singular, matrix and an error-contaminated right-hand side. The latter represents the available blur and noise-contaminated image, while the matrix models the blurring. Computation of a meaningful restoration of the available image requires the use of a regularization method. We consider the situation when the blurring matrix has a Kronecker product structure and an estimate of the norm of the desired image is available, and illustrate that efficient restoration of the available image can be achieved by Tikhonov regularization based on the global Lanczos method, and by using the connection of the latter to Gauss-type quadrature rules.

Solving equations through particle dynamics

The present work evaluates a recently developed particle method (DFPM). The basic idea behind this method is to utilize a Newtonian system of interacting particles that through dissipation solves mathematical problems. We find that this second order dynamical system results in an algorithm that is among the best methods known. The present work studies large systems of linear equations. Of special interest is the wide eigenvalue spectrum. This case is common as the discretization of the continuous problem becomes dense. The convergence rate of DFPM is shown to be in parity with that of the conjugate gradient method, both analytically and through numerical examples. However, an advantage with DFPM is that it is cheaper per iteration. Another advantage is that it is not restricted to symmetric matrices only, as is the case for the conjugate gradient method. The convergence properties of DFPM are shown to be superior to the closely related approach utilizing only a first order dynamical system, and also to several other iterative methods in numerical linear algebra. The performance properties are understood and optimized by taking advantage of critically damped oscillators in classical mechanics. Just as in the case of the conjugate gradient method, a limitation is that all eigenvalues (spring constants) are required to be of the same sign. DFPM has no other limitation such as matrix structure or a spectral radius as is common among iterative methods. Examples are provided to test the particle algorithm’s merits and also various performance comparisons with existent numerical algorithms are provided.

Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method

We have investigated the use of the adaptive high-order finite-element method (FEM) for geoelectromagnetic modeling. Because high-order FEM is challenging from the numerical and computational points of view, most published finite-element studies in geoelectromagnetics use the lowest order formulation. Solution of the resulting large system of linear equations poses the main practical challenge. We have developed a fully parallel and distributed robust and scalable linear solver based on the optimal block-diagonal and auxiliary space preconditioners. The solver was found to be efficient for high finite element orders, unstructured and nonconforming locally refined meshes, a wide range of frequencies, large conductivity contrasts, and number of degrees of freedom (DoFs). Furthermore, the presented linear solver is in essence algebraic; i.e., it acts on the matrix-vector level and thus requires no information about the discretization, boundary conditions, or physical source used, making it readily efficient for a wide range of electromagnetic modeling problems. To get accurate solutions at reduced computational cost, we have also implemented goal-oriented adaptive mesh refinement. The numerical tests indicated that if highly accurate modeling results were required, the high-order FEM in combination with the goal-oriented local mesh refinement required less computational time and DoFs than the lowest order adaptive FEM.

Wavelet-In-Time Multigrid-In-Space Preconditioning of Parabolic Evolution Equations

Two space-time variational formulations of linear parabolic evolution equations are considered: one is symmetric and elliptic on the trial space, while the other is not. In each case, a space-time Petrov--Galerkin discretization using suitable tensor product trial and test functions leads to a large linear system of equations. The well-posedness of this system with respect to parabolic norms induces a canonical preconditioner for the algebraic equations that arise after a choice of basis. For the iterative resolution of this algebraic system with parallelization in the temporal direction, we propose a sparse algebraic wavelet-in-time transformation on possibly nonuniform temporal meshes. This transformation approximately block-diagonalizes the preconditioner, and the individual spatial blocks can then be inverted, for instance, by standard spatial multigrid methods in parallel. The performance of the preconditioner is documented in a series of numerical experiments.

Restarted Pulay mixing for efficient and robust acceleration of fixed-point iterations

We present a variant of the restarted Pulay's Direct Inversion in the Iterative Subspace (DIIS) method for efficiently and robustly accelerating the convergence of fixed-point iterations. Specifically, we propose a simple modification of DIIS without any additional parameters, which we refer to as the r-Pulay method. We demonstrate the efficacy of r-Pulay in the context of the Jacobi iteration for solving large linear systems of equations, as well as in the Self Consistent Field (SCF) approach for Density Functional Theory (DFT) calculations. Overall, we find r-Pulay to be an attractive version of the restarted DIIS method.

Some properties of range restricted GMRES methods

The GMRES method is one of the most popular iterative schemes for the solution of large linear systems of equations with a square nonsingular matrix. GMRES-type methods also have been applied to the solution of linear discrete ill-posed problems. Computational experience indicates that for the latter problems variants of the standard GMRES method, that require the solution to live in the range of a positive power of the matrix of the linear system of equations to be solved, generally yield more accurate approximations of the desired solution than standard GMRES. This paper investigates properties of these variants of GMRES.

An Efficient Parallel Gauss-Seidel Algorithm on a 3D Torus Network-on-Chip

Network-on-chip (NoC) multi-core architectures with a large number of processing elements are becoming a reality with the recent developments in technology. In these modern systems the processing elements are interconnected with regular NoC topologies such as meshes and tori. In this paper we propose a parallel Gauss-Seidel (GS) iterative algorithm for solving large systems of linear equations on a 3-dimensional torus NoC architecture. The proposed parallel algorithm is O(Nn2/k3) time complexity for solving a system with a matrix of order n on a k×k×k 3D torus NoC architecture with N iterations assuming n and N are large compared to k. We show that under these conditions the proposed parallel GS algorithm has near optimal speedup.

An effective method to suppress numerical dispersion in 2D acoustic and elastic modelling using a high-order Padé approximation

In this paper, we propose a method that can effectively reduce the numerical dispersion for solving the acoustic and elastic-wave equations. The method is a fourth-order Pade approximation scheme, in which the time difference operator is a rational function and a block tridiagonal system needs to be solved at each step. On the one hand, to efficiently solve a large linear system of equations we propose an explicit method for this implicit algrithom. On the other hand, to approximate the high-order spatial derivatives we use an eighth-order stereo-modelling method using wavefield displacements and their gradients simultaneously. For this new method, we investigate some mathematical properties including the stability, errors and the numerical dispersion relationship for 1D and 2D cases. We also present numerical results computed by the Pade approximation and compare them with the eighth-order Lax–Wendroff correction method and the eighth-order staggered-grid method. Numerical results show that the high-order Pade approximation scheme can effectively suppress the numerical dispersion caused by discretizing the wave equations when coarse spatial grids are used or models have strong velocity contrasts between adjacent grids. In contrast to other high-order finite-difference methods, the new method takes substantially less computational time and requires less memory because large spatial and time increments can be used. Thus the high-order Pade approximation method can potentially be used to solve large-scale wave propagation problems and seismic tomography based on the wave equations.

Tikhonov regularization via flexible Arnoldi reduction

Flexible GMRES, introduced by Saad, is a generalization of the standard GMRES method for the solution of large linear systems of equations. It is based on the flexible Arnoldi process for reducing a large square matrix to a small matrix. We describe how the flexible Arnoldi process can be applied to implement one-parameter and multi-parameter Tikhonov regularization of linear discrete ill-posed problems. The method proposed is well suited for large-scale problems. Moreover, computed examples show that our method can give approximate solutions of higher accuracy than available direct methods for small-scale problems.

Space-Time Regularization for Video Decompression

We consider the problem of reconstructing frames from a video which has been compressed using the video compressive sensing (VCS) method. In VCS data, each frame comes from first subsampling the original video data in space and then averaging the subsampled sequence in time. This results in a large linear system of equations whose inversion is ill-posed. We introduce a convex regularizer to invert the system, where the spatial component is regularized by the total variation seminorm, and the temporal component is regularized by enforcing sparsity on the difference between the spatial gradients of each frame. Since the regularizers are $L^1$-like norms, the model can be written in the form of an easy-to-solve saddle point problem. The saddle point problem is solved by the primal-dual algorithm, whose implementation calls for nearly pointwise operations (i.e., no direct linear inversion) and has a simple parallel version. Results show that our model decompresses videos more accurately than other popular mod...

Accurate mirror-based camera pose estimation with explicit geometric meanings

An accurate closed form solution is proposed to estimate camera pose by several mirrored reference object images acquired via a planar mirror under different unknown poses. Compared with state-of-the-art methods, our method is more accurate when there are more than three images and has explicit geometric meanings. This method also properly handles cases in which some of the mirror poses are parallel. The central idea is to minimize an error metric based on all reflections of rotation, which enables the camera rotation to be estimated directly by SVD of sum of mirrored camera rotations. After that, the camera translation is computed by solving a large system of linear equations to minimize object space collinearity error. Both synthesized data and real data experiments show the advantages of our approach.

On facilitating the use of HARDI in population studies by creating rotation-invariant markers.

We design and evaluate a novel method to compute rotationally invariant features using High Angular Resolution Diffusion Imaging (HARDI) data. These measures quantify the complexity of the angular diffusion profile modeled using a higher order model, thereby giving more information than classical diffusion tensor-derived parameters. The method is based on the spherical harmonic (SH) representation of the angular diffusion information, and is generalizable to a range of HARDI reconstruction models. These scalars are obtained as homogeneous polynomials of the SH representation of a HARDI reconstruction model. We show that finding such polynomials is equivalent to solving a large linear system of equations, and present a numerical method based on sparse matrices to efficiently solve this system. Among the solutions, we only keep a subset of algebraically independent polynomials, using an algorithm based on a numerical implementation of the Jacobian criterion. We compute a set of 12 or 25 rotationally invariant measures representative of the underlying white matter for the rank-4 or rank-6 spherical harmonics (SHs) representation of the apparent diffusion coefficient (ADC) profile, respectively. Synthetic data was used to investigate and quantify the difference in contrast. Real data acquired with multiple repetitions showed that within subject variation in the invariants was less than the difference across subjects - facilitating their use to study population differences. These results demonstrate that our measures are able to characterize white matter, especially complex white matter found in regions of fiber crossings and hence can be used to derive new biomarkers for HARDI and can be used for HARDI-based population analysis.

Large Scale Finite Element Analysis Via Assembly-Free Deflated Conjugate Gradient

Large-scale finite element analysis (FEA) with millions of degrees of freedom (DOF) is becoming commonplace in solid mechanics. The primary computational bottleneck in such problems is the solution of large linear systems of equations. In this paper, we propose an assembly-free version of the deflated conjugate gradient (DCG) for solving such equations, where neither the stiffness matrix nor the deflation matrix is assembled. While assembly-free FEA is a well-known concept, the novelty pursued in this paper is the use of assembly-free deflation. The resulting implementation is particularly well suited for large-scale problems and can be easily ported to multicore central processing unit (CPU) and graphics-programmable unit (GPU) architectures. For demonstration, we show that one can solve a 50 × 106 degree of freedom system on a single GPU card, equipped with 3 GB of memory. The second contribution is an extension of the “rigid-body agglomeration” concept used in DCG to a “curvature-sensitive agglomeration.” The latter exploits classic plate and beam theories for efficient deflation of highly ill-conditioned problems arising from thin structures.

Polynomial-based mean weighted residuals methods for elliptic problems with nonlocal boundary conditions in the rectangle

In this paper, polynomial-based mean weighted residuals methods for the solution of elliptic problems with nonlocal boundary conditions, in rectangular domains, are numerically studied. When these methods are employed, the nonclassical boundary conditions involve the solution of large systems of linear equations or least squares problems, hence some numerical techniques for these solvers are compared to show the importance of using efficient algorithms for this purpose. Different kinds of nodes are used to demonstrate how they can be employed to solve different numerical problems when large derivatives of the solution appear. We will also study how using extra precision and/or oversampling can often reduce the computational effort.These methods can also be combined with others (as, for example, finite difference or spline methods) to solve linear and nonlinear parabolic or hyperbolic partial differential equations. The numerical study of some techniques as those explained above can help to obtain significantly better numerical approximations with a smaller computational effort.

Hydraulic Analysis of Water Supply Networks Using a Modified Hardy Cross Method

There are different methods for the hydraulic analysis of water supply networks. In the solution process of most of these methods, a large system of linear equations is solved in each iteration. This usually requires a high computational effort. Hardy Cross method is one of the approaches that do not need such a process and may converge to the solution through scalar divisions. However, this method has two shortcomings: first, initial discharges should satisfy continuity equation at each node; second a large number of iterations are required to converge to solution. In this article an algorithm is suggested for initial discharges that are close to the final results while the continuity equations are automatically established. This algorithm may be directly implemented in the Hardy Cross method. To reduce the number of iterations the Hardy Cross method is combined with third-order and sixteenth-order methods. The results of some numerical examples demonstrate that the use of the combined approach with the suggested initial guess reduces the number of iterations and hydraulic analysis time and the solutions converge with a high accuracy.

A coarse space for heterogeneous Helmholtz problems based on the Dirichlet-to-Neumann operator

The Helmholtz equation governing wave propagation and scattering phenomena is difficult to solve numerically. Its discretization with piecewise linear finite elements results in typically large linear systems of equations. The inherently parallel domain decomposition methods constitute hence a promising class of preconditioners. An essential element of these methods is a good coarse space. Here, the Helmholtz equation presents a particular challenge, as even slight deviations from the optimal choice can be devastating.In this paper, we present a coarse space that is based on local eigenproblems involving the Dirichlet-to-Neumann operator. Our construction is completely automatic, ensuring good convergence rates without the need for parameter tuning. Moreover, it naturally respects local variations in the wave number and is hence suited also for heterogeneous Helmholtz problems. The resulting method is parallel by design and its efficiency is demonstrated on 2D homogeneous and heterogeneous numerical examples.

A Multiple-FPGA parallel computing architecture for real-time simulation of soft-object deformation

Hardware-based parallel computing is proposed for acceleration of finite-element (FE) analysis of linear elastic deformation models. An implementation of the Preconditioned Conjugate Gradient algorithm on N Field Programmable Gate Array (FPGA) devices solves the large linear system of equations arising from the FE discretization. The system employs a large number of customized fixed-point computing units with a high-throughput memory architecture. An implementation of this scalable architecture on four Altera EP3SE110 FPGA devices yields a peak performance of 604 Giga Operations per second. This enables haptic simulation of a 3-dimensional deformable object of 21000 elements at an update rate of 400Hz.

Vertex-Ball Spring Smoothing: An efficient method for unstructured dynamic hybrid meshes

Spring analogy approach is one of the most popular dynamic mesh deformation methods. In the ball-vertex method, perpendicular linear springs are introduced to deal with the element collapse problem. However, it is not very efficient as a large system of linear equations has to be resolved. In order to overcome this difficulty, the Vertex-Ball Spring Smoothing algorithm (VerBSS) is proposed in this paper. Following the mesh smoothing concept, a sub-spring system derived from the “ball-vertex” model is built and solved on a node by node basis using an LDLT solver. Interior nodes are smoothed layer by layer in an iterative manner to achieve the best result. Parallel scheme is also introduced in the smoothing process for further improvement of the efficiency. Numerical examples in two and three dimensions show that VerBSS is much more efficient than the ball-vertex method, and is capable of dealing with practical engineering objects with complex geometries subject to large deformations. VerBSS can be applied to complicated mesh topologies as well, not only to 2D/3D dynamic mesh, but also to the hybrid dynamic mesh.

Performance evaluation of kernel fusion BLAS routines on the GPU: iterative solvers as case study

Programmers usually implement iterative methods that solve partial differential equations by expressing them using a sequence of basic kernels from libraries optimized for the graphics processing unit (GPU). The global runtime of the resulting combination is often penalized by the smallest and most inefficient vector operations. To improve the GPU exploitation, we identify and analyze the potential kernels to be fused according to the data dependence, data type and size, and GPU resources. This paper provides an extensive analysis of the impact of fusing vector operations [level 1 of Basic Linear Algebra Subprograms (BLAS)] on the performance of the GPU. The experimental evaluation shows that this optimization provides noticeable improvement especially for kernels with lower memory requirements and on more modern GPUs. It is worth noting that the fused BLAS operations can be very useful to help programmers efficiently code iterative methods to solve large linear systems of equations for the GPU. Iterative methods such as biconjugate gradient method (BCG) are one of the examples that can benefit from this optimization strategy. Indeed, kernel fusion of vector routines makes the most efficient GPU implementation of BCG run between $$1.09\times $$ 1.09 × and $$1.27\times $$ 1.27 × faster on three GPUs of different characteristics.

Comment on “Plasmonic nanoparticle monomers and dimers: from nanoantennas to chiral metamaterials” by Chigrin et al. in Appl Phys B (2011) 105:81–97

Chigrin et al. (Appl Phys B 105:81, 2011) have recently presented a general direct solution of the integral equation for the electric field inside the scatterer, which does not require inversion of integral operator or, equivalently, solution of large system of linear equations. This comment points out that the rigorous applicability of the considered derivation is limited to ellipsoids much smaller than the wavelength.