Ill-conditioned Linear Systems Research Articles

Analysis and Imaging in Magnetic Induction Tomography using the Impedance Method

This article discusses the utilization of the impedance method in computation of the forward problem in magnetic induction tomography (MIT). The algorithms for the inverse problem were also developed. The new approach for solving the resulting ill-conditioned linear system of the inverse problem is proposed and the quality of images obtained is discussed based on a quality index proposed in the literature. The results show the prevalence of TAS in relation to linear correlation between the real image and obtained image. With respect to contrast TRT prevalece in relation TAS. The indices of average luminance presents similarity for the both methods. TAS prevails for smaller objects and TRT for larger objects, showing the greater robustness of TAS.

Accurate solution of dense linear systems, part I: Algorithms in rounding to nearest

We investigate how extra-precise accumulation of dot products can be used to solve ill-conditioned linear systems accurately. For a given p-bit working precision, extra-precise evaluation of a dot product means that the products and summation are executed in 2p-bit precision, and that the final result is rounded into the p-bit working precision. Denote by u=2−p the relative rounding error unit in a given working precision. We treat two types of matrices: first up to condition number u−1, and second up to condition number u−2. For both types of matrices we present two types of methods: first for calculating an approximate solution, and second for calculating rigorous error bounds for the solution together with the proof of non-singularity of the matrix of the linear system. In the first part of this paper we present algorithms using only rounding to nearest, in Part II we use directed rounding to obtain better results. All algorithms are given in executable Matlab code and are available from my homepage.

Accurate solution of dense linear systems, Part II: Algorithms using directed rounding

In Part I and this Part II of our paper we investigate how extra-precise evaluation of dot products can be used to solve ill-conditioned linear systems rigorously and accurately. In Part I only rounding to nearest is used. In this Part II we improve the results significantly by permitting directed rounding. Linear systems with tolerances in the data are treated, and a comfortable way is described to compute error bounds for extremely ill-conditioned linear systems with condition numbers up to about u−2/n, where u denotes the relative rounding error unit in a given working precision. We improve a method by Hansen/Bliek/Rohn/Ning/Kearfott/Neumaier. Of the known methods by Krawczyk, Rump, Hansen et al., Ogita and Nguyen we show that our presented Algorithm LssErrBnd seems the best compromise between accuracy and speed. Moreover, for input data with tolerances, a new method to compute componentwise inner bounds is presented. For not too wide input data they demonstrate that the computed inclusions are often almost optimal. All algorithms are given in executable Matlab code and are available from my homepage.

A numerical scheme for nonlinear Schrödinger equation by MQ quasi-interpolation

Quasi-interpolation is a very powerful tool in the field of approximation theory and its applications, which can avoid solving large scale ill-conditioned linear system arising in approximating an unknown function by means of radial basis functions. In this paper, we use an univariate multi-quadrics(MQ) quasi-interpolation scheme to solve one-dimensional nonlinear Schrödinger equation. In this novel numerical scheme, the spatial derivatives are approximated by using the derivative of the quasi-interpolation and the temporal derivative is approximated by finite difference method. The main advantage of this proposed scheme is its simplicity. Two numerical examples are given and compared with the finite difference method (FDM) to verify the good accuracy and easy implementation of this method.

Balancing Neumann–Neumann methods for the cardiac Bidomain model

Balancing Neumann---Neumann preconditioners are constructed, analyzed and numerically studied for the cardiac Bidomain model in three-dimensions. This reaction---diffusion system is discretized by low-order finite elements in space and implicit---explicit methods in time, yielding very ill-conditioned linear systems that must be solved at each time step. The proposed algorithm is based on decomposing the domain into nonoverlapping subdomains and on solving iteratively the Bidomain Schur complement obtained by implicitly eliminating the degrees of freedom interior to each subdomain. The iteration is preconditioned by a Balancing Neumann---Neumann method employing local Neumann solves on each subdomain and a coarse Bidomain solve. A novel approach for the estimation of the average operator of the nonoverlapping decomposition provides a framework for designing coarse spaces for Balancing Neumann---Neumann methods. The theoretical estimates obtained show that the proposed method is scalable, quasi-optimal and robust with respect to possible coefficient discontinuities of the Bidomain operator. The results of extensive parallel numerical tests in three dimensions confirm the convergence rates predicted by the theory.

A numerical solution of a two-dimensional IHCP

In this paper, we will first study the existence and uniqueness of the solution of a two-dimensional inverse heat conduction problem (IHCP) which is severely ill-posed, i.e., the solution does not depend continuously on the data. We propose a stable numerical approach based on the finite-difference method and the least-squares scheme to solve this problem in the presence of noisy data. We prove the convergence of the numerical solution, then to regularize the resultant ill-conditioned linear system of equations, we apply the Tikhonov regularization 0th, 1st and 2nd method to obtain the stable numerical approximation to the solution. The stability and accuracy of the scheme presented is evaluated by comparison with the Singular Value Decomposition (SVD) method.

An interior-point algorithm for the minimization arising from 3D contact problems with friction

The paper deals with a variant of the interior-point method for the minimization of strictly quadratic objective function subject to simple bounds and separable quadratic inequality constraints. Such minimizations arise from the finite element approximation of contact problems of linear elasticity with friction in three space dimensions. The main goal of the paper is the convergence analysis of the algorithm and its implementation. The optimal preconditioners for solving ill-conditioned inner linear systems are proposed. Numerical experiments illustrate the computational efficiency for large-scale problems.

A FAST LAGRANGE METHOD FOR LARGE-SCALE IMAGE RESTORATION PROBLEMS WITH REFLECTIVE BOUNDARY CONDITION

The goal of the image restoration is to find a good approximation of the original image for the degraded image, the blurring matrix, and the statistics of the noise vector given. Fast truncated Lagrange (FTL) method has been proposed by G. Landi as a image restoration method for large-scale ill-conditioned BTTB linear systems([3]). We implemented FTL method for the image restoration problem with reflective boundary condition which gives better reconstructions of the unknown, the true image.

A finite element multifrontal method for 3D CSEM modeling in the frequency domain

There has been a recent increase in the use of controlled-source electromagnetic (CSEM) surveys in the exploration for oil and gas. We developed a modeling scheme for 3D CSEM modeling in the frequency domain. The electric field was decomposed in primary and secondary components to eliminate the singularity originated by the source term. The primary field was calculated using a closed form solution, and the secondary field was computed discretizing a second-order partial differential equation for the electric field with the edge finite element. The solution to the linear system of equations was obtained using a massive parallel multifrontal solver, because such solvers are robust for indefinite and ill-conditioned linear systems. Recent trends in parallel computing were investigated for their use in mitigating the computational overburden associated with the use of a direct solver, and of its feasibility for 3D CSEM forward modeling with the edge finite element. The computation of the primary field was parallelized, over the computational domain and the number of sources, using a hybrid model of parallelism. When using a direct solver, the attainment of multisource solutions was only competitive if the same factors are used to achieve a solution for multi right-hand sides. This aspect was also investigated using the presented methodology. We tested our proposed approach using 1D and 3D synthetic models, and they demonstrated that it is robust and suitable for 3D CSEM modeling using a distributed memory system. The codes could thus be used to help design new surveys, as well to estimate subsurface conductivities through the implementation of an appropriate inversion scheme.

Inverse heat conduction problems in three-dimensional anisotropic functionally graded solids

A meshless method based on the local Petrov–Galerkin approach is applied to inverse transient heat conduction problems in three-dimensional solids with continuously inhomogeneous and anisotropic material properties. The Heaviside step function is used as a test function in the local weak form, leading to the derivation of local integral equations. Nodal points are randomly distributed in the domain analyzed, and each node is surrounded by a spherical subdomain in which a local integral equation is applied. A meshless approximation based on the moving least-squares method is employed in the implementation. After performing spatial integrations, we obtain a system of ordinary differential equations for certain nodal unknowns. A backward finite-difference method is used for the approximation of the diffusive term in the heat conduction equation. A truncated singular-value decomposition is used to solve the ill-conditioned linear system of algebraic equations at each time step. The effectiveness of the meshless local Petrov–Galerkin (MLPG) method for this inverse problem is demonstrated by numerical examples.

A Fast Direct Solver for Structured Linear Systems by Recursive Skeletonization

We present a fast direct solver for structured linear systems based on multilevel matrix compression. Using the recently developed interpolative decomposition of a low-rank matrix in a recursive manner, we embed an approximation of the original matrix into a larger, but highly structured sparse one that allows fast factorization and application of the inverse. The algorithm extends the Martinsson/Rokhlin method developed for 2D boundary integral equations and proceeds in two phases: a precomputation phase, consisting of matrix compression and factorization, followed by a solution phase to apply the matrix inverse. For boundary integral equations which are not too oscillatory, e.g., based on the Green's functions for the Laplace or low-frequency Helmholtz equations, both phases typically have complexity O(N) in two dimensions, where $N$ is the number of discretization points. In our current implementation, the corresponding costs in three dimensions are $O(N^{3/2})$ and $O(N \log N)$ for precomputation and solution, respectively. Extensive numerical experiments show a speedup of $\sim 100$ for the solution phase over modern fast multipole methods; however, the cost of precomputation remains high. Thus, the solver is particularly suited to problems where large numbers of iterations would be required. Such is the case with ill-conditioned linear systems or when the same system is to be solved with multiple right-hand sides. Our algorithm is implemented in Fortran and freely available.

Fast Algorithms for Image Reconstruction with Application to Partially Parallel MR Imaging

This paper presents two fast algorithms for total variation–based image reconstruction in a magnetic resonance imaging technique known as partially parallel imaging (PPI), where the inversion matrix is large and ill-conditioned. These algorithms utilize variable splitting techniques to decouple the original problem into more easily solved subproblems. The first method reduces the image reconstruction problem to an unconstrained minimization problem, which is solved by an alternating proximal minimization algorithm. One phase of the algorithm solves a total variation (TV) denoising problem, and the second phase solves an ill-conditioned linear system. Linear and sublinear convergence results are given, and an implementation based on a primal-dual hybrid gradient (PDHG) scheme for the TV problem and on a Barzilai–Borwein scheme for the linear inversion is proposed. The second algorithm exploits the special structure of the PPI reconstruction problem by decomposing it into one subproblem involving Fourier transforms and another subproblem that can be treated by the PDHG scheme. Numerical results and comparisons with recently developed methods indicate the efficiency of the proposed algorithms.

On the Importance of the Stokes-Brinkman Equations for Computing Effective Permeability in Karst Reservoirs

AbstractCavities and fractures significantly affect the flow paths in carbonate reservoirs and should be accurately accounted for in numerical models. Herein, we consider the problem of computing the effective permeability of rock samples based on high-resolution 3D CT scans containing millions of voxels. We use the Stokes-Brinkman equations in the entire domain, covering regions of free flow governed by the Stokes equations, porous Darcy flow, and transitions between them. The presence of different length scales and large (ten orders of magnitude) contrasts in permeability leads to highly ill-conditioned linear systems of equations, which are difficult to solve. To obtain a problem that is computationally tractable, we first analyze the relative importance of the Stokes and Darcy terms for a set of idealized 2D models. We find that, in terms of effective permeability, the Stokes-Brinkman equations are only applicable for a special parameter set where the effective free-flow permeability is less than four orders of magnitude different from the matrix permeability. All other cases can be accurately modeled with either the Stokes or the Darcy end-member flows, depending on if there do or do not exist percolating free-flow regions. The insights obtained are used to perform a direct computation of the effective permeability of a rock sample model with more than 8 million cells.

A rational Arnoldi approach for ill-conditioned linear systems

For the solution of full-rank ill-posed linear systems a new approach based on the Arnoldi algorithm is presented. Working with regularized systems, the method theoretically reconstructs the true solution by means of the computation of a suitable function of matrix. In this sense, the method can be referred to as an iterative refinement process. Numerical experiments arising from integral equations and interpolation theory are presented. Finally, the method is extended to work in connection with the standard Tikhonov regularization with the right-hand side contaminated by noise.

Numerical results for linear Fredholm integral equations of the first kind over surfaces in three dimensions

Linear Fredholm integral equations of the first kind over surfaces are less familiar than those of the second kind, although they arise in many applications like computer tomography, heat conduction and inverse scattering. This article emphasizes their numerical treatment, since discretization usually leads to ill-conditioned linear systems. Strictly speaking, the matrix is nearly singular and ordinary numerical methods fail. However, there exists a numerical regularization method – the Tikhonov method – to deal with this ill-conditioning and to obtain accurate numerical results.

A new iterative refinement of the solution of ill-conditioned linear system of equations

In this paper, a new iterative refinement of the solution of an ill-conditioned linear system of equations are given. The convergence properties of the method are studied. Some numerical experiments of the method are given and compared with that of two of the available methods.

Fast Structured AMG Preconditioning for the Bidomain Model in Electrocardiology

The electrical activity of the heart may be modeled by a nonlinear system of partial differential equations known as the bidomain model. Due to the rapid variations in the electrical field, accurate simulations require a fine-scale discretization of the equations and consequently the solution of large severely ill-conditioned linear systems at each time step. Solving these systems is a major bottleneck of the whole simulation. We propose a highly effective preconditioning strategy for a general and popular three-dimensional formulation of the problem. A theoretical analysis of the preconditioned matrix ensuring mesh independence of the spectrum is also described. Numerical comparisons with state-of-the-art approaches confirm the effectiveness of our preconditioning technique. Finally, we show that an equivalent but less exercised formulation provides the best performance, in terms of CPU time.

An inexact interior point method for L 1-regularized sparse covariance selection

Sparse covariance selection problems can be formulated as log-determinant (log-det) semidefinite programming (SDP) problems with large numbers of linear constraints. Standard primal–dual interior-point methods that are based on solving the Schur complement equation would encounter severe computational bottlenecks if they are applied to solve these SDPs. In this paper, we consider a customized inexact primal–dual path-following interior-point algorithm for solving large scale log-det SDP problems arising from sparse covariance selection problems. Our inexact algorithm solves the large and ill-conditioned linear system of equations in each iteration by a preconditioned iterative solver. By exploiting the structures in sparse covariance selection problems, we are able to design highly effective preconditioners to efficiently solve the large and ill-conditioned linear systems. Numerical experiments on both synthetic and real covariance selection problems show that our algorithm is highly efficient and outperforms other existing algorithms.

Iterative solvers for Tikhonov regularization of dense inverse problems

According to the special demands arising from the development of science and technology, in the last decades appeared a special class of problems that are inverse to the classical direct ones. Such an inverse problem is concerned with the opposite way, usually followed by a direct one: finding the cause of a given effect or finding the law of evolution given the cause and effect. Very frequently, such inverse problems are modelled by Fredholm first-kind integral equations that give rise after discretization to (very) ill-conditioned linear systems, in classical or least squares formulation. Then, an efficient numerical solution can be obtained by using the Tikhonov regularization technique. In this respect, in the present paper, we propose three Kovarik-like algorithms for numerical solution of the regularized problem. We prove convergence for all three methods and present numerical experiments on a mathematical model of an inverse problem concerned with the determination of charge distribution generating a given electric field.

Application of the shifted-Laplace preconditioner for iterative solution of a higher order finite element discretisation of the vector wave equation: First experiences

The analysis of electromagnetic scattering by electrically large one-side open–ended cavities remains a challenge. Finite element discretisation of the vector wave equation to solve for the electric field inside the cavity leads to ill-conditioned indefinite linear systems of large dimension, a result of the requirement of a fine nearly uniform discrete sampling in the computational domain. Direct methods based on frontal solution techniques to solve the resulting linear system have been used because efficient iterative methods were not available. With the arrival of the shifted-Laplace preconditioner of Erlangga, iterative solution of the indefinite system becomes tractable. This paper discusses the modifications required for application of the shifted-Laplace preconditioner to cavity scattering and some preliminary results of this approach. The shifted-Laplace preconditioner is shown to be very effective for improving the convergence rate of the iterative solution algorithm. However, to be able to handle problems with larger number of degrees of freedom, it is necessary to include a multigrid algorithm to solve the preconditioner system, as this will allow the use of short recurrence Krylov subspace methods, as opposed to the currently employed long-recurrence method, for which the storage requirements of the Krylov basis become unpractically large.