Large Systems Of Linear Equations Research Articles

Parameter estimation with model order reduction for elliptic differential equations

In this paper a new method for parameter estimation for elliptic partial differential equations is introduced. Parameter estimation includes minimizing an objective function, which is a measure for the difference between the parameter-dependent solution of the differential equation and some given data. We assume, that the given data results in a good approximation of the state of the system.In order to evaluate the objective function the solution of a differential equation has to be computed and hence, a large system of linear equations has to be solved. Minimization methods involve many evaluations of the objective function and therefore, the differential equation has to be solved multiple times. Thus, the computing time for parameter estimation can be large. Model order reduction was developed in order to reduce the computational effort of solving these differential equations multiple times. We use the given approximation of the state of the system as reduced basis and omit computing any snapshots. Therefore, our approach decreases the effort of the offline phase drastically. Furthermore, the dimension of the reduced system is one and thus, is much smaller than the dimension of other approaches. However, the obtained reduced model is a good approximation only close to the given data. Hence, the reduced system can lead to large errors for parameter sets, which correspond to solutions far away from the given approximation of the state of the system. In order to prevent convergence of the parameter estimator to such a local minimizer we penalise the approximation error between the full and the reduced system.

Global Golub–Kahan bidiagonalization applied to large discrete ill-posed problems

We consider the solution of large linear systems of equations that arise from the discretization of ill-posed problems. The matrix has a Kronecker product structure and the right-hand side is contaminated by measurement error. Problems of this kind arise, for instance, from the discretization of Fredholm integral equations of the first kind in two space-dimensions with a separable kernel and in image restoration problems. Regularization methods, such as Tikhonov regularization, have to be employed to reduce the propagation of the error in the right-hand side into the computed solution. We investigate the use of the global Golub–Kahan bidiagonalization method to reduce the given large problem to a small one. The small problem is solved by employing Tikhonov regularization. A regularization parameter determines the amount of regularization. The connection between global Golub–Kahan bidiagonalization and Gauss-type quadrature rules is exploited to inexpensively compute bounds that are useful for determining the regularization parameter by the discrepancy principle.

Slope tomography based on eikonal solvers and the adjoint-state method

Velocity macromodel building is a crucial step in the seismic imaging workflow as it provides the necessary background model for migration or full waveform inversion. In this study, we present a new formulation of stereotomography that can handle more efficiently long-offset acquisition, complex geological structures and large-scale data sets. Stereotomography is a slope tomographic method based upon a semi-automatic picking of local coherent events. Each local coherent event, characterized by its two-way traveltime and two slopes in common-shot and common-receiver gathers, is tied to a scatterer or a reflector segment in the subsurface. Ray tracing provides a natural forward engine to compute traveltime and slopes but can suffer from non-uniform ray sampling in presence of complex media and long-offset acquisitions. Moreover, most implementations of stereotomography explicitly build a sensitivity matrix, leading to the resolution of large systems of linear equations, which can be cumbersome when large-scale data sets are considered. Overcoming these issues comes with a new matrix-free formulation of stereotomography: a factored eikonal solver based on the fast sweeping method to compute first-arrival traveltimes and an adjoint-state formulation to compute the gradient of the misfit function. By solving eikonal equation from sources and receivers, we make the computational cost proportional to the number of sources and receivers while it is independent of picked events density in each shot and receiver gather. The model space involves the subsurface velocities and the scatterer coordinates, while the dips of the reflector segments are implicitly represented by the spatial support of the adjoint sources and are updated through the joint localization of nearby scatterers. We present an application on the complex Marmousi model for a towed-streamer acquisition and a realistic distribution of local events. We show that the estimated model, built without any prior knowledge of the velocities, provides a reliable initial model for frequency-domain FWI of long-offset data for a starting frequency of 4 Hz, although some artefacts at the reservoir level result from a deficit of illumination. This formulation of slope tomography provides a computationally efficient alternative to waveform inversion method such as reflection waveform inversion or differential-semblance optimization to build an initial model for pre-stack depth migration and conventional FWI.

Deconvolution based photoacoustic reconstruction with sparsity regularization.

In most photoacoustic tomography (PAT) reconstruction approaches, it is assumed that the receiving transducers have omnidirectional response and can fully surround the region of interest. These assumptions are not satisfied in practice. To deal with these limitations, we present a novel deconvolution based photoacoustic reconstruction with sparsity regularization (DPARS) technique. The DPARS algorithm is a semi-analytical reconstruction approach in which the projections of the absorber distribution derived from a deconvolution-based method are computed and used to generate a large linear system of equations. In these projections, computed over limited viewing angles, the directivity effect of the transducer is taken into account. The distribution of absorbers is computed using a sparse representation of absorber coefficients obtained from the discrete cosine transform. This sparse representation helps improve the numerical conditioning of the system of equations and reduces the computation time of the deconvolution-based approach by one order of magnitude relative to Tikhonov regularization. The algorithm has been tested in simulations, and using two-dimensional and three-dimensional experimental data obtained with a conventional ultrasound transducer. The results show that DPARS, when evaluated using contrast-to-noise ratio and root-mean-square errors, outperforms the conventional delay-and-sum (DAS) reconstruction method.

On the stability of ADI methods

When solving parabolic equations in two space dimensions implicit methods are preferred to the explicit method because of their better stability properties. Straightforward implementation of implicit methods require time-consuming solution of large systems of linear equations, and ADI methods are preferred instead. We expect the ADI methods to inherit the stability properties of the implicit methods they are derived from, and we demonstrate that this is partly true. The Douglas-Rachford and Peaceman-Rachford methods are absolutely stable in the sense that their growth factors are ≤ 1 in absolute value. Near jump discontinuities, however, there are differences w.r.t. how the ADI methods react to the situation: do they produce oscillations and how effectively do they damp them. We demonstrate the behaviour on two simple examples.

A Preconditioned Low-Rank Projection Method with a Rank-Reduction Scheme for Stochastic Partial Differential Equations

In this study, we consider the numerical solution of large systems of linear equations obtained from the stochastic Galerkin formulation of stochastic partial differential equations (PDEs). We propose an iterative algorithm that exploits the Kronecker product structure of the linear systems. The proposed algorithm efficiently approximates the solutions in low-rank tensor format. Using standard Krylov subspace methods for the data in tensor format is computationally prohibitive due to the rapid growth of tensor ranks during the iterations. To keep tensor ranks low over the entire iteration process, we devise a rank-reduction scheme that can be combined with the iterative algorithm. The proposed rank-reduction scheme identifies an important subspace in the stochastic domain and compresses tensors of high rank on-the-fly during the iterations. The proposed reduction scheme is a coarse-grid method in that the important subspace can be identified inexpensively in a coarse spatial grid setting. The efficiency of the proposed method is illustrated by numerical experiments on benchmark problems.

Approximate model for geometrical complex structures

Many engineering structures are nowadays made of composite materials or metal foam. These modern engineering materials contain very complex inner geometry. To simulate the deformational behaviour of these structures often requires a high number of discretisation elements. This in turn yields a very large system of linear equations that are extremely time and memory consuming or practically impossible to solve. It is therefore desirable to find an approach to overcome this obstacle. In this paper a numerical method is proposed to find an approximate substitute model for geometrical complex structures.

Parallelization of the Rational Arnoldi Algorithm

Rational Krylov methods are applicable to a wide range of scientific computing problems, and the rational Arnoldi algorithm is a commonly used procedure for computing an orthonormal basis of a rational Krylov space. Typically, the computationally most expensive component of this algorithm is the solution of a large linear system of equations at each iteration. We explore the option of solving several linear systems simultaneously, thus constructing the rational Krylov basis in parallel. If this is not done carefully, the basis being orthogonalized may become badly conditioned, leading to numerical instabilities in the orthogonalization process. We introduce the new concept of continuation pairs, which gives rise to a near-optimal parallelization strategy that allows us to control the growth of the condition number of this nonorthogonal basis. As a consequence we obtain a significantly more accurate and reliable parallel rational Arnoldi algorithm. The computational benefits are illustrated using several numerical examples from different application areas.

Hydrodynamics of suspensions of passive and active rigid particles: a rigid multiblob approach

We develop a rigid multiblob method for numerically solving the mobility problem for suspensions of passive and active rigid particles of complex shape in Stokes flow in unconfined, partially confined, and fully confined geometries. As in a number of existing methods, we discretize rigid bodies using a collection of minimally-resolved spherical blobs constrained to move as a rigid body, to arrive at a potentially large linear system of equations for the unknown Lagrange multipliers and rigid-body motions. Here we develop a block-diagonal preconditioner for this linear system and show that a standard Krylov solver converges in a modest number of iterations that is essentially independent of the number of particles. For unbounded suspensions and suspensions sedimented against a single no-slip boundary, we rely on existing analytical expressions for the Rotne-Prager tensor combined with a fast multipole method or a direct summation on a Graphical Processing Unit to obtain an simple yet efficient and scalable implementation. For fully confined domains, such as periodic suspensions or suspensions confined in slit and square channels, we extend a recently-developed rigid-body immersed boundary method to suspensions of freely-moving passive or active rigid particles at zero Reynolds number. We demonstrate that the iterative solver for the coupled fluid and rigid body equations converges in a bounded number of iterations regardless of the system size. We optimize a number of parameters in the iterative solvers and apply our method to a variety of benchmark problems to carefully assess the accuracy of the rigid multiblob approach as a function of the resolution. We also model the dynamics of colloidal particles studied in recent experiments, such as passive boomerangs in a slit channel, as well as a pair of non-Brownian active nanorods sedimented against a wall.

Analyzing Barrettes as Large-Section Supports by CCT

Most of soil structure interaction methods for analyzing large-section supports such as barrette foundation modeling the barrette and surrounding soil using 3D FE model. In which, the model leads to a large finite element mesh of a large system of linear equations to be solved. In this paper, a Composed Coefficient Technique (CCT) is adapted for analyzing barrette. The technique takes into account the 3D full interactions between barrette and the surrounding soil. Due to the high rigidity of the barrette relative to the surrounding soil, a uniform or variable settlement along the barrette height can be considered. This enables to compose the stiffness coefficients of the soil matrix into composed coefficients, which consequently leads to a significant reduction in the soil stiffness matrix. An application for analyzing barrette by CCT is carried on the soil of the new area of East Port Said, in where the typical soil stratification is very week and structures in this area need to be supported by deep foundations such as barrettes. The application presents guidelines and diagrams for barrettes that may be used in East Port Said.

SEBSM-Based Iterative Method for Solving Large Systems of Linear Equations and Its Applications in Engineering Computation

In this paper, a new iterative method, for solving sparse nonsymmetrical systems of linear equations is proposed based on the Simultaneous Elimination and Back-Substitution Method (SEBSM), and the method is applied to solve systems resulted in engineering problems solved using Finite Element Method (FEM). First, SEBSM is introduced for solving general linear systems using the direct method. And, then an iterative method based on SEBSM is presented. In the method, the coefficient matrix [Formula: see text] is split into lower, diagonally banded and upper matrices. The iterative convergence can be controlled by selecting a suitable bandwidth of the diagonally banded matrix. And the size of the working array needing to be stored in iteration is as small as the bandwidth of the diagonally banded matrix. Finally, an accelerating strategy for this iterative method is proposed by introducing a relaxation factor, which can speed up the convergence effectively if an optimal relaxation factor is chosen. Two numerical examples are given to demonstrate the behavior of the proposed method.

Adaptive FEM-based nonrigid image registration using truncated hierarchical B-splines

We present an efficient approach of Finite Element Method (FEM)-based nonrigid image registration, in which the spatial transformation is constructed using truncated hierarchical B-splines (THB-splines). The image registration framework minimizes an energy functional using an FEM-based method and thus involves solving a large system of linear equations. This framework is carried out on a set of successively refined grids. However, due to the increased number of control points during subdivision, large linear systems are generated which are generally demanding to solve. Instead of using uniform subdivision, an adaptive local refinement scheme is carried out, only refining the areas of large change in deformation of the image. By incorporating the key advantages of THB-spline basis functions such as linear independence, partition of unity and reduced overlap into the FEM-based framework, we improve the matrix sparsity and computational efficiency. The performance of the proposed method is demonstrated on 2D synthetic and medical images.

Two pass AGA iterative methods to solve large systems of linear equations

Two pass AGA iterative methods to solve large systems of linear equations

Quantum algorithms and the finite element method

The finite element method is used to approximately solve boundary value problems for differential equations. The method discretises the parameter space and finds an approximate solution by solving a large system of linear equations. Here we investigate the extent to which the finite element method can be accelerated using an efficient quantum algorithm for solving linear equations. We consider the representative general question of approximately computing a linear functional of the solution to a boundary value problem, and compare the quantum algorithm's theoretical performance with that of a standard classical algorithm -- the conjugate gradient method. Prior work had claimed that the quantum algorithm could be exponentially faster, but did not determine the overall classical and quantum runtimes required to achieve a predetermined solution accuracy. Taking this into account, we find that the quantum algorithm can achieve a polynomial speedup, the extent of which grows with the dimension of the partial differential equation. In addition, we give evidence that no improvement of the quantum algorithm could lead to a super-polynomial speedup when the dimension is fixed and the solution satisfies certain smoothness properties.

Preconditioning techniques for an image deblurring problem

SummaryIn this paper, we consider the solution of a large linear system of equations, which is obtained from discretizing the Euler–Lagrange equations associated with the image deblurring problem. The coefficient matrix of this system is of the generalized saddle point form with high condition number. One of the blocks of this matrix has the block Toeplitz with Toeplitz block structure. This system can be efficiently solved using the minimal residual iteration method with preconditioners based on the fast Fourier transform. Eigenvalue bounds for the preconditioner matrix are obtained. Numerical results are presented. Copyright © 2016 John Wiley & Sons, Ltd.

Filtering Chromatic Aberration for Wide Acceptance Angle Electrostatic Lenses II--Experimental Evaluation and Software-Based Imaging Energy Analyzer.

Here, the experimental results of the method of filtering the effect of chromatic aberration for wide acceptance angle electrostatic lens-based system are described. This method can eliminate the effect of chromatic aberration from the images of a measured spectral image sequence by determining and removing the effect of higher and lower kinetic energy electrons on each different energy image, which leads to significant improvement of image and spectral quality. The method is based on the numerical solution of a large system of linear equations and equivalent with a multivariate strongly nonlinear deconvolution method. A matrix whose elements describe the strongly nonlinear chromatic aberration-related transmission function of the lens system acts on the vector of the ordered pixels of the distortion free spectral image sequence, and produces the vector of the ordered pixels of the measured spectral image sequence. Since the method can be applied not only on 2D real- and $k$ -space diffraction images, but also along a third dimension of the image sequence that is along the optical or in the 3D parameter space, the energy axis, it functions as a software-based imaging energy analyzer (SBIEA). It can also be applied in cases of light or other type of optics for different optical aberrations and distortions. In case of electron optics, the SBIEA method makes possible the spectral imaging without the application of any other energy filter. It is notable that this method also eliminates the disturbing background significantly in the present investigated case of reflection electron energy loss spectra. It eliminates the instrumental effects and makes possible to measure the real physical processes better.

RGIsearch: A C++ program for the determination of renormalization group invariants

RGIsearch is a C++ program that searches for invariants of a user-defined set of renormalization group equations. Based on the general shape of the β-functions of quantum field theories, RGIsearch searches for several types of invariants that require different methods. Additionally, it supports the computation of invariants up to two-loop level. A manual for the program is given, including the settings and set-up of the program, as well as a test case. Program summaryProgram title: RGIsearchCatalogue identifier: AEZQ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEZQ_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 129911No. of bytes in distributed program, including test data, etc.: 8364946Distribution format: tar.gzProgramming language: C++.Computer: Desktop PC.Operating system: Linux.RAM: Ranging from several MB to several GBClassification: 2.9, 4.8, 5.Nature of problem: The determination of renormalization group invariants, which can be used as a probe for the high-energy behavior of theories in particle physics.Solution method: Several types of invariants are considered based on the general shape of the renormalization group equations. The problem of computing these invariants is then reduced to creating and solving very large systems of linear equations or generalized eigenvalue problems. Using the specific sparsity structure of these systems, the systems can be solved.Restrictions: Since the algorithms that solve the linear systems and generalized eigenvalue problems are very efficient, the main restriction is the amount of available RAM, where the systems are stored. The program supports polynomial renormalization group equations, which are typical particle physics, up to two-loop order.Running time: Ranging from seconds to minutes

Sampling method based projection approach for the reconstruction of 3D acoustically penetrable scatterers

We present a projection based regularization parameter choice approach within the framework of the linear sampling method for the reconstruction of acoustically penetrable objects. Using the Golub–Kahan bidiagonalization algorithm and the Lanczos tridiagonalization process we form appropriate subspaces which generate a sequence of regularized solutions. As a result two new and efficient methods are developed and used for the solution of problems that involve large linear systems of equations. The effectiveness of our approach is illustrated with reconstructions of three dimensional objects.

Quantum algorithms: an overview

AbstractQuantum computers are designed to outperform standard computers by running quantum algorithms. Areas in which quantum algorithms can be applied include cryptography, search and optimisation, simulation of quantum systems and solving large systems of linear equations. Here we briefly survey some known quantum algorithms, with an emphasis on a broad overview of their applications rather than their technical details. We include a discussion of recent developments and near-term applications of quantum algorithms.

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.