Row-action Methods Research Articles

On the randomized multiple row-action methods for solving linear least-squares problems

On the randomized multiple row-action methods for solving linear least-squares problems

Correction to: Survey of a class of iterative row-action methods: The Kaczmarz method

Correction to: Survey of a class of iterative row-action methods: The Kaczmarz method

Survey of a class of iterative row-action methods: The Kaczmarz method

AbstractThe Kaczmarz algorithm is an iterative method that solves linear systems of equations. It stands out among iterative algorithms when dealing with large systems for two reasons. First, at each iteration, the Kaczmarz algorithm uses a single equation, resulting in minimal computational work per iteration. Second, solving the entire system may only require the use of a small subset of the equations. These characteristics have attracted significant attention to the Kaczmarz algorithm. Researchers have observed that randomly choosing equations can improve the convergence rate of the algorithm. This insight led to the development of the Randomized Kaczmarz algorithm and, subsequently, several other variations emerged. In this paper, we extensively analyze the native Kaczmarz algorithm and many of its variations using large-scale systems as benchmarks. Through our investigation, we have verified that, for consistent systems, various row sampling schemes can outperform both the original and Randomized Kaczmarz method. Specifically, sampling without replacement and using quasirandom numbers are the fastest techniques. However, for inconsistent systems, the Conjugate Gradient method for Least-Squares problems overcomes all variations of the Kaczmarz method for these types of systems.

On the behaviour of the underrelaxed Hildreth’s row-action method for computing projections onto polyhedra

On the behaviour of the underrelaxed Hildreth’s row-action method for computing projections onto polyhedra

A wavelet-based sparse row-action method for image reconstruction in magnetic particle imaging.

Magnetic particle imaging (MPI) is a preclinical imaging technique capable of visualizing the spatio-temporal distribution of magnetic nanoparticles. The image reconstruction of this fast and dynamic process relies on efficiently solving an ill-posed inverse problem. Current approaches to reconstruct the tracer concentration from its measurements are either adapted to image characteristics of MPI but suffer from higher computational complexity and slower convergence or are fast but lack in the image quality of the reconstructed images. In this work we propose a novel MPI reconstruction method to combine the advantages of both approaches into a single algorithm. The underlying sparsity prior is based on an undecimated wavelet transform and is integrated into a fast row-action framework to solve the corresponding MPI minimization problem. Its performance is numerically evaluated against a classical FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) approach on simulated and real MPI data. The experimental results show that the proposed method increases image quality with significantly reduced computation times. In comparison to state-of-the-art MPI reconstruction methods, our approach shows better reconstruction results and at the same time accelerates the convergence rate of the underlying row-action algorithm.

Sampled limited memory methods for massive linear inverse problems

In many modern imaging applications the desire to reconstruct high resolution images, coupled with the abundance of data from acquisition using ultra-fast detectors, have led to new challenges in image reconstruction. A main challenge is that the resulting linear inverse problems are massive. The size of the forward model matrix exceeds the storage capabilities of computer memory, or the observational dataset is enormous and not available all at once. Row-action methods that iterate over samples of rows can be used to approximate the solution while avoiding memory and data availability constraints. However, their overall convergence can be slow. In this paper, we introduce a sampled limited memory row-action method for linear least squares problems, where an approximation of the global curvature of the underlying least squares problem is used to speed up the initial convergence and to improve the accuracy of iterates. We show that this limited memory method is a generalization of the damped block Kaczmarz method, and we prove linear convergence of the expectation of the iterates and of the error norm up to a convergence horizon. Numerical experiments demonstrate the benefits of these sampled limited memory row-action methods for massive 2D and 3D inverse problems in tomography applications.

An improved HASM method for dealing with large spatial data sets

Surface modeling with very large data sets is challenging. An efficient method for modeling massive data sets using the high accuracy surface modeling method (HASM) is proposed, and HASM_Big is developed to handle very large data sets. A large data set is defined here as a large spatial domain with high resolution leading to a linear equation with matrix dimensions of hundreds of thousands. An augmented system approach is employed to solve the equality-constrained least squares problem (LSE) produced in HASM_Big, and a block row action method is applied to solve the corresponding very large matrix equations. A matrix partitioning method is used to avoid information redundancy among each block and thereby accelerate the model. Experiments including numerical tests and real-world applications are used to compare the performances of HASM_Big with its previous version, HASM. Results show that the memory storage and computing speed of HASM_Big are better than those of HASM. It is found that the computational cost of HASM_Big is linearly scalable, even with massive data sets. In conclusion, HASM_Big provides a powerful tool for surface modeling, especially when there are millions or more computing grid cells.

Hildreth’s algorithm with applications to soft constraints for user interface layout

Hildreth’s algorithm is a row action method for solving large systems of inequalities. This algorithm is efficient for problems with sparse matrices, as opposed to direct methods such as Gaussian elimination or QR-factorization. We apply Hildreth’s algorithm, as well as a randomized version, along with prioritized selection of the inequalities, to efficiently detect the highest priority feasible subsystem of equations. We prove convergence results and feasibility criteria for both cyclic and randomized Hildreth’s algorithm, as well as a mixed algorithm which uses Hildreth’s algorithm for inequalities and Kaczmarz algorithm for equalities. These prioritized, sparse systems of inequalities commonly appear in constraint-based user interface (UI) layout specifications. The performance and convergence of these proposed algorithms are evaluated empirically using randomly generated UI layout specifications of various sizes. The results show that these methods offer improvements in performance over standard methods like Matlab’s LINPROG, a well-known efficient linear programming solver, and the recent developed Kaczmarz algorithm with prioritized IIS detection.

Generalized row-action methods for tomographic imaging

Row-action methods play an important role in tomographic image reconstruction. Many such methods can be viewed as incremental gradient methods for minimizing a sum of a large number of convex functions, and despite their relatively poor global rate of convergence, these methods often exhibit fast initial convergence which is desirable in applications where a low-accuracy solution is acceptable. In this paper, we propose relaxed variants of a class of incremental proximal gradient methods, and these variants generalize many existing row-action methods for tomographic imaging. Moreover, they allow us to derive new incremental algorithms for tomographic imaging that incorporate different types of prior information via regularization. We demonstrate the efficacy of the approach with some numerical examples.

Compressive sensing using the modified entropy functional

In most compressive sensing problems, ℓ1 norm is used during the signal reconstruction process. In this article, a modified version of the entropy functional is proposed to approximate the ℓ1 norm. The proposed modified version of the entropy functional is continuous, differentiable and convex. Therefore, it is possible to construct globally convergent iterative algorithms using Bregmanʼs row-action method for compressive sensing applications. Simulation examples with both 1D signals and images are presented.

A general extending and constraining procedure for linear iterative methods

Algebraic reconstruction techniques (ARTs), on both their successive and simultaneous formulations, have been developed since the early 1970s as efficient ‘row-action methods’ for solving the image-reconstruction problem in computerized tomography. In this respect, two important development directions were concerned with, first, their extension to the inconsistent case of the reconstruction problem and, second, their combination with constraining strategies, imposed by the particularities of the reconstructed image. In the first part of this paper, we introduce extending and constraining procedures for a general iterative method of an ART type and we propose a set of sufficient assumptions that ensure the convergence of the corresponding algorithms. As an application of this approach, we prove that Cimmino's simultaneous reflection method satisfies this set of assumptions, and we derive extended and constrained versions for it. Numerical experiments with all these versions are presented on a head phantom widely used in the image reconstruction literature. We also consider hard thresholding constraining used in sparse approximation problems and apply it successfully to a 3D particle image-reconstruction problem.

AIR Tools — A MATLAB package of algebraic iterative reconstruction methods

We present a MATLAB package with implementations of several algebraic iterative reconstruction methods for discretizations of inverse problems. These so-called row action methods rely on semi-convergence for achieving the necessary regularization of the problem. Two classes of methods are implemented: Algebraic Reconstruction Techniques (ART) and Simultaneous Iterative Reconstruction Techniques (SIRT). In addition we provide a few simplified test problems from medical and seismic tomography. For each iterative method, a number of strategies are available for choosing the relaxation parameter and the stopping rule. The relaxation parameter can be fixed, or chosen adaptively in each iteration; in the former case we provide a new “training” algorithm that finds the optimal parameter for a given test problem. The stopping rules provided are the discrepancy principle, the monotone error rule, and the NCP criterion; for the first two methods “training” can be used to find the optimal discrepancy parameter.

Extended and constrained diagonal weighting algorithm with application to inverse problems in image reconstruction

The algebraic reconstruction of images in computerized tomography gives rise to large, sparse and ill-conditioned inverse problems, in which the ‘effect’ (measurement of the attenuation of x-ray intensities after penetration of the analyzed body) is used to compute the ‘cause’ (the values of the attenuation function inside the body, i.e. the image). Algebraic reconstruction techniques (ART), on both their successive or simultaneous formulation, have been developed since the early 1970s as efficient ‘row action methods’ for solving the image reconstruction problem in computerized tomography. In this respect, two important development directions were concerned with their extension to the inconsistent case of the reconstruction problem as well as with their combination with constraining strategies, imposed by the particularities of the reconstructed image. In our paper we analyze, from these two points of view, the diagonal weighting (DW) algorithm proposed by Y Censor, D Gordon and R Gordon in 2001 as an improvement of the classical Cimmino's reflection method. In the first part of the paper we introduce general extended and constraining procedures for ART, based on a minimal set of sufficient assumptions that ensure the convergence of the corresponding algorithms. Starting from this general context we then design an extended form of the DW algorithm together with a constraining procedure for which we prove convergence under appropriate assumptions. Numerical experiments are presented on two phantoms widely used in the literature.

Delay-time tomography of the upper mantle below Europe, the Mediterranean, and Asia Minor

SUMMARY The method of linearized delay-time tomography is applied to retrieve the P velocity heterogeneity of the upper mantle in central Europe, the Mediterranean and Asia Minor. The upper mantle is parametrized in slowness cell-functions with lateral dimensions of 1. The model parametrization includes event mislocation vectors and station delay-time corrections. In all 20 000 parameters describe the model. Delay-time data, selected from the ISC data base (1964-1982), from regional events and from stations up to 90 of epicentral distance are used. About 18 000 events and 937 stations contributed nearly 500 000 delay times. We discuss the tomographic method, the data processing, the cell model employed, the computation of the ray geometry, the use of composite rays in tomography, and the cell hit count. The inversions are performed with two different iterative row-action methods, LSQR and SIRT. Results obtained-with the LSQR method will be shown only. The spatial resolution is investigated with sensitivity analyses using two different synthetic velocity models; a 'spike' model and a model with harmonically varying velocity anomalies. The spatial resolution is estimated to be high (on the order of the cell size) in the mantle below central and southeastern Europe. The mapped velocity heterogeneity exhibits primarily low velocities around a depth of 200 km and just above the transition from upper to lower mantle. This indicates a systematic deviation of the reference model (Jeffreys-Bullen) from the radially averaged earth structure, which may violate the basic assumptions underlying the linearized approach. A posteriori investigation of the data also indicates the presence of a low-velocity zone and the existence of strong velocity depth gradients.

Row-Action Inversion of the Barrick–Weber Equations

Abstract The Barrick–Weber equations describe the interaction of radar signals with the dynamic ocean surface, and so provide a mathematical basis for oceanic remote sensing. This report considers the inversion of these equations with several of the row-action methods commonly used to solve large linear systems with unstructured sparsity. It is found that the performance of the methods in inverting both synthetic and measured Doppler spectral data is comparable, with the method of Chahine–Twomey–Wyatt offering a slight advantage in the reliability of the recovery of the full directional wave spectrum and of parameters derived from its integration. Some remarks and open questions on the ill-posedness of the inversion problem conclude the paper.

An Open Question on Cyclic Relaxation

The problem discussed in this note is highly interesting. It is related to several dual iterative methods, such as the methods proposed by Kaczmarz, Hildreth, Agmon, Cryer, Mangasarian, Herman, Lent, Censor, and others. Cast as ‘row-action methods’ these algorithms have been proved as useful tools for solving large convex feasibility problems that arise in medical image reconstruction from projections, in inverse problems in radiation therapy, and in linear programming. The question that we want to answer is how these algorithms behave when the feasible region is empty. It is shown that under certain conditions the primal sequence still converges while the dual sequence {y k } obeys the rule y k =u k +k v, where {u k } is a converging sequence and v is a fixed vector that satisfies A T v=0,v≥0,and,b T v>0. It is conjectured that these properties hold whenever the feasible region is empty. However, the validity of this claim remains an open question.

Maximum Entropy Reconstruction Methods in Electron Paramagnetic Resonance Imaging

Electron Paramagnetic Resonance (EPR) is a spectroscopic technique that detects and characterizes molecules with unpaired electrons (i.e., free radicals). Unlike the closely related nuclear magnetic resonance (NMR) spectroscopy, EPR is still under development as an imaging modality. Athough a number of physical factors have hindered its development, EPR's potential is quite promising in a number of important application areas, including in vivo oximetry. EPR images are generally reconstructed using a tomographic imaging technique, of which filtered backprojection (FBP) is the most commonly used. We apply two iterative methods for maximum-entropy image reconstruction in EPR. The first is the multiplicative algebraic reconstruction technique (MART), a well-known row-action method. We propose a second method, known as LSEnt (least-squares entropy), that maximizes entropy and performs regularization by maintaining a desired distance from the measurements. LSEnt is in part motivated by the barrier method of interior-point programming. We present studies in which images of two physical phantoms, reconstructed using FBP, MART, and LSEnt, are compared. The images reconstructed using MART and LSEnt have lower variance, better contrast recovery, subjectively better resolution, and reduced streaking artifact than those reconstructed using FBP. These results suggest that maximum-entropy reconstruction methods (particularly the more flexible LSEnt) may be critical in overcoming some of the physical challenges of EPR imaging.

The adventures of a simple algorithm

This paper presents, analyzes, and tests a row-action method for solving the regularized linear programming problem: minimize 1/2∥ x ∥ 2 2+α c T x subject to A x ⩾ b , where α is a given positive constant. The proposed method is based on the observation that the dual of this problem has the form: maximize b T y −1/2∥A T y −α c ∥ 2 2 subject to y ⩾ 0 , and if y α solves the dual then the point x α=A T y α−α c provides the unique solution of the primal. Maximizing the dual objective function by changing one variable at a time results in an effective row-action method which is suitable for solving large sparse problems. One aim of this paper is to clarify the convergence properties of the proposed scheme. Let y k denote the estimated dual solution at the end of the kth iteration, and let x k=A T y k−α c denote the corresponding primal estimate. It is proved that the sequence { x k} converges to x α , while the sequence { y k} converges to a point y α that solves the dual. The only assumption which is needed in order to establish these claims is that the feasible region is not empty. Yet perhaps the more striking features of the algorithm are related to temporary situations in which it attempts to solve an inconsistent system. In such cases the sequence { y k} obeys the rule: y k= u k+k v , where { u k} is a fastly converging sequence and v is a fixed vector that satisfies A T v = 0 and b T v >0 . So the sequence { x k} is almost unchanged for many iterations. The paper ends with numerical experiments that illustrate the effects of this phenomenon and the close links with Kaczmarz's SOR method.

Approximation of linear programs by Bregman's DF projections

The motivation of this paper is twofold: to contribute to the theory of Bregman's D F projections and to point out its link to finding ϵ-optimal solutions to linear programming problems. A symmetric primal–dual pair of the D F projection problem is presented, we call it Young programming problem, because the usual inequality that relates the primal and dual objective functions reduces to the well-known Young inequality. Some basic properties of Young programs are derived and an associated linear program is defined in a natural way. It will be shown that an ϵ-optimal solution of the associated linear program is obtained by solving an ϵ-parametric Young program and the optimal solutions of ϵ-parametric Young programs converge to an optimal solution of the associated linear program when ϵ approaches 0. It may also be interesting that the explicit formulation of the dual program enables us to give a dual interpretation of row-action methods, a powerful tool for solving D F projection problems in general, and allows us to find new functions satisfying sufficient conditions for the convergence of row-action methods.

Asymptotically optimal row-action methods for generalized least squares problems

Recently row-action methods have been frequently used in actual applications, for instance, image reconstruction. The asymptotic convergence behavior of row-action methods has been studied by lusem and De Pierro (1990) and Mandel (1984). Here asymptotically optimal behavior of row-action methods are investigated. Some asymptotically optimal row-action methods, for example, asymptotically optimal Kaczmarz and Hildreth methods, are proposed. Several numerical algorithms of asymptotically optimal row-action methods for generalized least squares problems are given.