Smoothing and Regularization with Modified Sparse Approximate Inverses

Abstract

Sparse approximate inverses 𝑀 which satisfy min_𝑀‖𝐴𝑀−𝐼‖_𝐹 have shown to be an attractive alternative to classical smoothers like Jacobi or Gauss-Seidel (Tang and Wan; 2000). The static and dynamic computation of a SAI and a SPAI (Grote and Huckle; 1997), respectively, comes along with advantages like inherent parallelism and robustness with equal smoothing properties (Bröker et al.; 2001). Here, we are interested in developing preconditioners that can incorporate probing conditions for improving the approximation relative to high- or low-frequency subspaces. We present analytically derived optimal smoothers for the discretization of the constant-coefficient Laplace operator. On this basis, we introduce probing conditions in the generalized Modified SPAI (MSPAI) approach (Huckle and Kallischko; 2007) which yields efficient smoothers for multigrid. In the second part, we transfer our approach to the domain of ill-posed problems to recover original information from blurred signals. Using the probing facility of MSPAI, we impose the preconditioner to act as approximately zero on the noise subspace. In combination with an iterative regularization method, it thus becomes possible to reconstruct the original information more accurately in many cases. A variety of numerical results demonstrate the usefulness of this approach.

1. Introduction

For applying an iterative solution method to an ill-conditioned system of linear equations 𝐴𝑥=𝑏 with sparse matrix 𝐴∈ℝ^𝑛×𝑛, 𝑥,𝑏∈ℝ^𝑛, it is often crucial to include an efficient preconditioner. Here, the original problem 𝐴𝑥=𝑏 is replaced by the preconditioned system 𝑀𝐴𝑥=𝑀𝑏 or 𝐴𝑥=𝐴(𝑀𝑦)=𝑏. Often, used preconditioners as Jacobi, Gauss-Seidel, or Incomplete LU (ILU) decomposition are of unsatisfactory quality or strongly sequential. In a parallel environment, both the computation of the preconditioner 𝑀 as well as the application of the preconditioner on any given vector 𝑣 should be efficient. Furthermore, an iterative solver applied on 𝐴𝑀𝑦=𝑏 or 𝑀𝐴𝑥=𝑀𝑏 should converge much faster than for 𝐴𝑥=𝑏 (e.g., it holds cond(𝑀𝐴)≪cond(𝐴)).

The first two conditions can be easily satisfied by using a sparse matrix 𝑀 as approximation to 𝐴⁻¹. Note that the inverse of a sparse system 𝐴 is nearly dense, but in many cases, the entries of 𝐴⁻¹ are rapidly decaying, so most of the entries are very small (see Demko et al. [1]).

Sparse approximate inverses can be computed by minimizing 𝐴𝑀−𝐼 in the Frobenius norm, where 𝐼 denotes the identity matrix. In this Frobenius norm minimization, we can include further approximation conditions, described by the Modified SPAI (MSPAI) [2, 3] method. This additional feature allows us to control the approximation property of the preconditioner. So, by means of probing vectors, we can choose subspaces for which the preconditioner satisfies certain conditions. This is especially important for iterative solution methods that differ between high-frequency and low-frequency components, for example, smoothing in multigrid or regularization techniques based on iterative solvers.

The outline of the paper is the following. In Section 1, we will give a survey of SPAI and MSPAI and a short description of multigrid methods and iterative solvers for regularization problems. In Section 2, we show that for multigrid methods the smoothing property can be greatly improved by using MSPAI in comparison to SAI or SPAI smoothers. With a different subspace approach, we focus on the reconstruction of signals associated to ill-posed problems in Section 3. We present numerical results throughout the paper to demonstrate the impact of our approach at the corresponding parts. A conclusion with a short outlook closes the discussion.

Benson and Frederickson [4] were the first to propose the computation of an explicit approximation 𝑀≈𝐴⁻¹ to the inverse of a system matrix 𝐴∈ℝ^𝑛×𝑛. For an a priori prescribed sparsity pattern 𝒫, this can be done in a static way by solvingmin_{𝒫(𝑀)=𝒫}‖𝐴𝑀−𝐼‖²_𝐹=^𝑛_𝑘=1min_{𝒫𝑚𝑘=𝒫𝑘}𝐴𝑚_𝑘−𝑒_𝑘²₂, (1) with 𝑚_𝑘 the 𝑘th column of the preconditioner 𝑀 and 𝑒_𝑘 the 𝑘th column of the identity matrix 𝐼. 𝒫(𝐴) denotes the pattern 𝒫 of 𝐴 and 𝒫_𝑘 the pattern of the 𝑘th column of 𝒫. This well-known approach of providing a SAI preconditioner, naturally leads to inherent parallelism which is one of its main advantages. Each of the 𝑛 small least squares problems regarding one column can be computed independently of one another.

Using both the index set 𝒥_𝑘, which is implicitly given by 𝒫_𝑘 and contains the indices 𝑗 such that 𝑚_𝑘(𝑗)≠0, and its corresponding so-called shadow ℐ_𝑘, that is the indices of nonzero rows in 𝐴(∶,𝒥_𝑘), each subproblem in (1) is related to a small matrix 𝐴_𝑘∶=𝐴(ℐ_𝑘,𝒥_𝑘) if 𝒫_𝑘 is sparse. We refer to the reduced sparse column vectors as 𝑚_𝑘∶=𝑚_𝑘(𝒥_𝑘) and 𝑒_𝑘̂∶=𝑒_𝑘(ℐ_𝑘), respectively. The solution of every reduced problem min𝑚_𝑘‖𝐴𝑚_𝑘−𝑒_𝑘̂‖₂ can be obtained, for example, by using QR decomposition using Householder or the modified Gram-Schmidt algorithm.

1.1. The SPAI Algorithm

The SPAI algorithm is an additional feature in this Frobenius norm minimization that introduces different strategies for choosing new profitable indices in 𝑚_𝑘 to improve on an already computed approximation. We assume that by solving (1) for a given index set 𝒥_𝑘, we already have determined an optimal solution 𝑚_𝑘(𝒥_𝑘) inducing the sparse vector 𝑚_𝑘 with residual 𝑟_𝑘. Dynamically, we want to augment new entries in 𝑚_𝑘 and solve (1) for this enlarged index set 𝒥_𝑘 such that we derive a reduction in the norm of the new residual 𝑟_𝑘̃=𝐴(ℐ_𝑘,𝒥_𝑘)𝑚_𝑘(𝒥_𝑘)−𝑒_𝑘(ℐ_𝑘).

Following Cosgrove et al. [5] and Grote and Huckle [6], in SPAI, we test one possible new index 𝑗∈𝒥_new out of a given set of possible new indices 𝒥_new to improve 𝑚_𝑘. Therefore, we can consider the reduced 1D minimization problemmin_𝜆𝑗𝐴𝑚_𝑘+𝜆_𝑗𝑒_𝑗−𝑒_𝑘2=min_𝜆𝑗𝜆_𝑗𝐴_𝑗+𝑟_𝑘2. (2) The solution of (2) is given by 𝜆_𝑗𝑟=−^𝑇𝑘𝐴𝑒_𝑗𝐴e_𝑗²₂, (3) and leads to an improved squared residual norm 𝜌²_𝑗=𝑟_𝑘²−𝑟^𝑇_𝑘𝐴𝑒_𝑗²𝐴𝑒_𝑗²₂. (4) For improving 𝑚_𝑘, we only have to consider indices 𝑗 in rows of 𝐴 that are related to nonzero entries in the old residual 𝑟_𝑘; otherwise, they do not lead to a reduction of the residual norm. Thus, we have to determine those column indices 𝑗, which satisfy 𝑟^𝑇_𝑘𝐴𝑒_𝑗≠0. Let us denote the index set of nonzero entries in 𝑟_𝑘 by ℒ. By 𝒥_𝑖, we denote the set of new indices that are related to the nonzero elements in the 𝑖th row of 𝐴, and by 𝒥_new∶=⋃_𝑖∈ℒ𝒥_𝑖 the set of all possible new indices that will lead to a reduction of the residual norm. The one or more newly added indices 𝒥_𝑐 are chosen to be a subset of 𝒥_new that corresponds to the maximal reduction in 𝑟_𝑘. For the enlarged index set 𝒥_𝑘∪𝒥_𝑐, we have to update the QR decomposition of the related least squares submatrix and solve the new column 𝑚_𝑘.

It is possible to influence the sparsity and approximation quality during the computation of the SPAI indirectly by different parameters, for example, how many entries are to be added in one step, how many pattern updates are to be done, which residual norm should be reached, or which initial pattern 𝒫 is to be used. Note that SPAI can also be applied on dense systems to compute a sparse preconditioner.

1.2. Modified SPAI

Holland et al. [7] have generalized the SPAI ansatz allowing a sparse target matrix 𝐵 on the right hand side in the form min_{𝒫(𝑀)=𝒫}‖𝐴𝑀−𝐵‖_𝐹. This approach is useful in connection with some kind of two-level preconditioning: first, compute a standard sparse preconditioner 𝐵 for 𝐴, and then improve this preconditioner by an additional Frobenius norm minimization with target 𝐵. From an algorithmic point of view, the minimization with sparse target matrix 𝐵, instead of 𝐼, introduces no additional difficulties. Simply, 𝒫(𝑀) should be chosen more carefully with respect to 𝐴 and 𝐵.

In [2], we combine this approach with classical probing techniques [8–10], which are, for example, applied to preconditioning Schur complements. In contrast to classical probing, our basic formulationmin_{𝒫(𝑀)=𝒫}‖𝐶𝑀−𝐵‖_𝐹=min_{𝒫(𝑀)=𝒫}𝐶₀𝜌𝑒^𝑇𝐵𝑀−⁰𝜌𝑓^𝑇_F, (5) with sparse matrices 𝐶₀ and 𝐵₀, is not restricted to special probing subspaces as it allows any choice of 𝑒 and 𝑓. The resulting preconditioner 𝑀 satisfies both 𝐶₀𝑀≈𝐵₀ and 𝑒^𝑇𝑀≈𝑓^𝑇. We refer to the first 𝑛 rows of (5), that is, 𝐶₀𝑀−𝐵₀, as full approximation part and to the additional rows as probing part. The weight 𝜌≥0 enables us to control how much emphasis is put on the probing constraints, and the matrices 𝑒,𝑓∈ℝ^𝑛×𝑘 represent the 𝑘-dimensional subspace on which the preconditioner should be optimal. Choosing 𝜌=0, 𝐶₀=𝐴, and 𝐵₀=𝐼 in (5) leads to the classical SPAI formulation. Setting 𝐶₀=𝐼 and 𝐵₀=𝐴, we end up with a formulation computing explicit sparse approximations to 𝐴. In this case, the derived approximation on 𝐴 can have a considerably fewer number of nonzeros (nnz) than 𝐴 but choosing 𝑓^𝑇=𝑒^𝑇𝐴, the preconditioner 𝑀 will have a similar action on 𝑒^𝑇 as 𝐴.

Furthermore, it is also possible to include individual probing conditions for each column 𝑚_𝑘. As a new approach, we use probing masks defined by sparse row vectors 𝑠_𝑘, 𝑘∈{1,…,𝑛}, containing the same pattern as 𝑚_𝑘, that is, 𝒫(𝑠_𝑘)=𝒫(𝑚_𝑘), and add the condition min𝑚_𝑘|𝑠_𝑘̂𝑚_𝑘−𝑓_𝑘| to the Frobenius norm minimization of the 𝑘th column. Corresponding to 𝑚_𝑘, 𝑠_𝑘̂ denotes the reduced form of 𝑠_𝑘. The masks for each column of 𝑀 can be stored in a sparse rectangular matrix 𝑆, whereby the individual probing conditions result from diag(𝑆𝑀)≈𝑓. Compared to MSPAI using global probing vectors 𝑒,𝑓, this approach gives considerably more freedom for choosing probing conditions individually for each column of the preconditioner. Note that, for example, for a tridiagonal pattern, the mask 𝑠_𝑘̂=(1,0,−1) can be used to enforce a quasisymmetry in the column vector 𝑚_𝑘 such that 𝑚_{𝑘−1,𝑘}≈𝑚_𝑘+1,𝑘.

The field of applications using MSPAI is versatile: we can improve preconditioners resulting from ILU, IC, FSAI, FSPAI, or AINV (see [11] for an overview) by adding probing information. We also overcome the main drawbacks of MILU and classical probing such as the restriction to certain vectors like (1,1,…,1)^𝑇 as probing subspace and the rather difficult efficient implementation on parallel computers. The numerical examples in [2, 3] demonstrate MSPAI’s effectiveness for preconditioning various PDE matrices and preconditioning Schur complements arising from domain decompositions.

1.3. Multigrid

The crucial observation leading to multigrid (MG) methods is the following: applying a stationary iterative solver like Gauss-Seidel iteration gives a satisfactory reduction of the error on the subspace related to high-frequency components. Therefore, Gauss-Seidel iteration is considered as a smoother in a first step of the MG algorithm. The error mainly contains smooth components and can be projected to a smaller linear system. This reduced system can be tackled recursively by the same approach based on smoothing steps and projection on an again reduced system. On the coarsest level, the small linear system can be solved explicitly. Afterwards, the coarse solutions have to be prolongated step by step back to the finer levels including postsmoothing steps.

Here, we are mainly interested in the smoother. To derive a convergent method the smoother has to reduce the error in the high-frequency subspace. For a given matrix 𝐴 and an approximate inverse smoother 𝑀, the iteration 𝑘 is described by 𝑥^(𝑘+1)=𝑥^(𝑘)+𝑀(𝑏−𝐴𝑥^(𝑘)), and thus the error is given by 𝐼−𝑀𝐴. As the eigenvalues and eigenvectors are analytically well-known for a discretization of the constant coefficient Laplace operator, it is possible to fully discuss the convergence behavior of MG in this special case. For analyzing the smoothing property, the eigenvalues are separated into high- and low-frequency eigenvectors. The projection 𝑃=𝑈^𝑇(𝐼−𝑀𝐴)𝑈 on the high-frequency subspace 𝑈 gives the smoothing factor defined as the spectral radius of 𝑃. In the constant coefficient case, the smoothing factor can also be described as 1−𝑚𝑎, where 𝑎 and 𝑚 are generating functions representing 𝐴 and 𝑀 [12, 13]. Note that the technique of generating functions is similar to Local Fourier Analysis (LFA) used in the standard multigrid literature [14]. In 2D the functions are defined in 𝑥,𝑦∈[0,𝜋], where the high-frequency domain is given by the difference 𝐺 between the two squares [0,𝜋]² and [0,𝜋/2]² with corners at 𝜋0,2,𝜋2,𝜋2,𝜋2,0,(𝜋,0),(𝜋,𝜋),(0,𝜋). (6) Thus, the smoothing factor is given by 𝜇∶=max_{𝑥,𝑦∈𝐺}1−𝑚𝑥,𝑦𝑎𝑥,𝑦. (7) In this paper, we want to use MSPAI as a smoother. SAI was already considered by Tang and Wan in [15] and SPAI by Bröker et al. in [16].

1.4. Iterative Regularization

For ill-posed problems, as they arise in image restoration, regularization techniques are important in order to recover the original information. Let us consider the model problem𝑥⟶^blur𝐻𝑥⟶^noise𝐻𝑥+𝜂=𝑏, (8) where 𝑥 is the original image, 𝐻 is the blur operator, 𝜂 is a vector representing the noise, and 𝑏 is the observed image. We want to recover 𝑥 as good as possible and as fast as possible. Because 𝐻 may be extremely ill-conditioned or even singular, and, because of the presence of noise, (8) cannot be solved directly. Consequently, to solve 𝐻𝑥=𝑏 on the signal subspace, a regularization technique has to be applied. One of the classical methods is the Tikhonov regularization [17] which solves min_𝑥‖𝐻𝑥−𝑏‖²₂+𝛾‖𝑏‖²₂⟺𝐻^𝑇𝐻+𝛾𝐼𝑥=𝐻^𝑇𝑏, (9) instead of (8) for a fixed regularization parameter 𝛾≥0.

Another regularization method is based on an iterative solver such as the Conjugate Gradient (CG) method [8, 18] for spd matrices or CG on the normal equations in the general unsymmetric case. The usual observation which coincides with the CG convergence analysis is that in the first iterations, the error is reduced relative to large eigenvalues. In later steps, the eigenspectrum related to noise and small eigenvalues dominates the evolution of the approximate solution. Therefore, the restoration has to stop after a few iterations before the method starts to reduce the error relative to the noise space.

Preconditioning usually should accelerate the convergence without destroying the quality of the reconstruction [19–21]. Structured preconditioners like Toeplitz or circulant matrices are considered typically whenever spatially invariant blur operators are treated. For general 𝐻, a preconditioner like ILU will lead to faster convergence, but the quality of the reconstruction will deteriorate, because the preconditioner also improves the solution relative to the unwanted noise subspace. Therefore, it is even more demanding to develop preconditioners for general 𝐻, for example, for spatially variant blur.

The application of the preconditioner can have three positive effects: (i)reduce the necessary number of iterations, (ii)result in a better reconstruction of the original vector, and (iii)result in a flat convergence curve such that it is easier to find the best reconstruction.

In general, we have to expect that not all three conditions can be reached simultaneously. Therefore, we have to present different preconditioners depending on the application.

Using preconditioners within iterative methods to restore original data has been successfully proposed by Nagy, for example, in [22, 23], mostly in connection with nearly structured problems. Furthermore, for the analysis and solution of discrete ill-posed systems there are several MATLAB packages such as Regularization Tools developed by Hansen [24] and RestoreTools developed by Nagy et al. [25].

2. MSPAI Smoothing

2.1. A 1D Model Problem

To derive approximate inverse smoothers, we consider the standard 1D discretized Laplace operator with constant coefficients which is of the form 𝐴₁𝑥=𝑏, with 𝐴₁∶=tridiag(−1/2,1,−1/2). The matrix 𝐴₁ is related to the generating function or symbol𝑎₁(𝑥)𝑒=1−^𝑖𝑥+𝑒^−𝑖𝑥2=1−cos(𝑥). (10) The high-frequency part of 𝐴₁ is represented by (10) for 𝑥∈[𝜋/2,𝜋]=∶𝐼₁. Hence, the optimal smoothing parameter 𝜔 in the Jacobi smoother 1−𝜔𝑎₁(𝑥) is found by solving the problemmin_𝜔max_𝑥∈𝐼11−𝜔𝑎₁(𝑥)=min_𝜔max_𝑥∈𝐼1|𝑠(𝑥)|=^𝜔=2/313. (11) The solution can be found by replacing (11) with the maximum over the two boundary values min_𝜔max{|𝑠(𝜋/2)|,|𝑠(𝜋)|}.

2.1.1. Analytical Derivation of the Optimal Smoother

As approximate inverse smoother, we choose a trigonometric polynomial of the same degree𝑚₁(𝑥)=𝑎+2𝑏cos(𝑥), (12) and a tridiagonal Toeplitz matrix 𝑀=tridiag(𝑏,𝑎,𝑏), respectively. The smoothing condition for (12) can be written as min_𝑎,𝑏max_𝑥∈𝐼11−𝑚₁(𝑥)𝑎₁(𝑥)=min_𝑎,𝑏max_𝑥∈𝐼1|𝑠(𝑥)|, (13) or using the boundary values of 𝐼₁ asmin_𝑎,𝑏max1−𝑚₁𝜋2𝑎₁𝜋2,1−𝑚₁(𝜋)𝑎₁(𝜋),1−𝑚₁(𝑢)𝑎₁(𝑢)=min_𝑎,𝑏max|𝑎−1|,|2(𝑎−2𝑏)−1|,1−𝑚₁(𝑢)𝑎₁(𝑢), (14) in which 𝑢∈𝐼₁ is the local extreme point of the function 1−𝑚₁(𝑥)𝑎₁(𝑥) with derivative equal to zero. We thus obtain for 𝑢=1/2−𝑎/4𝑏 the quadratic condition𝑠(𝑢)=(𝑎+2𝑏)²8𝑏−1, (15) and the overall solution (𝑎,𝑏)=(16/17,4/17) which leads to𝑀3_opt4∶=tridiag,1716,417.17 (16) The related smoothing factor for this optimal preconditioner 𝑚(𝑥) is 1/17=0.0588, which is significantly smaller than the optimal smoothing factor of 1/3=0.333 related to 𝜔=2/3 for the Jacobi smoother.

2.1.2. Individual Probing Masks

The minimization conditions (14) on 𝑎 and 𝑏 can be also seen as conditions for the entries of the matrix directly. A column 𝑘 of 𝑀 is described by the three values (𝑚_{𝑘−1,𝑘},𝑚_𝑘,𝑘,𝑚_𝑘+1,𝑘)^𝑇=∶𝑚_𝑘, where the main diagonal entry is related to 𝑎 and the upper and lower entries to 𝑏. Therefore, we can translate the two linear minimization conditions of (14) into individual probing masks 𝑠_𝑘̂ on entries of 𝑀 directly 𝑆3_𝐿1𝑚∶min_𝑘(0,1,0)𝑚_𝑘−1,𝑆3_𝐿2𝑚∶min_𝑘(−1,1,−1)𝑚_𝑘−12. (17) Notice that we refer to individual probing conditions as local probing conditions. There are several possibilities to eliminate and replace the quadratic condition (15) by linear conditions that can be related to masks itself. In one possible method, we assume that the linear conditions of (14) are satisfied exactly. Inserting 𝑎=1 into the second condition yields 𝑏=1/4. We replace the denominator of the quadratic condition with this value and get the new linear condition min_𝑎,𝑏(𝑎+2𝑏)²8⋅1/4−1=min_𝑎,𝑏(𝑎+2𝑏)²−2⟶min_𝑎,𝑏𝑎+2𝑏−√2, (18) which is related to the probing condition with isotropic probing mask 𝑆3_𝐿3𝑚∶min_𝑘(1,1,1)𝑚_𝑘−√2. (19) We call a mask 𝑠_𝑘̂ to be isotropic if the minimization min𝑚_𝑘|𝑠_𝑘̂𝑚_𝑘−𝑓_𝑘| can be described by one probing vector 𝑒∈ℝ^𝑛 with 𝑒(𝒥_𝑘)=𝛼𝑠^𝑇_𝑘̂, 𝛼∈ℝ and 𝑘∈{1,…,𝑛}, in the form min_{𝒫(𝑀)=𝒫}‖𝑒^𝑇𝑀−𝑓^𝑇‖₂. Here, for instance, 𝑠_𝑘̂=(−1,1,−1) of 𝑆3_𝐿2 is related to 𝑒=(1,−1,1,−1,…)^𝑇 and 𝑠_𝑘̂=(1,1,1) of 𝑆3_𝐿3 to 𝑒=(1,1,…,1)^𝑇.

In our following numerical experiments, we are interested in a comparison between MSPAI with probing conditions, the optimal smoother 𝑀3_opt (16), and the tridiagonal preconditioner given by SPAI. The SPAI matrix is derived by the minimization over the Frobenius norm in which the reduced Least Squares problem for a typical column of 𝑀 is given by𝑚min_𝑘𝑚−0.5001−0.50−0.51−0.50−0.5100−0.5^𝑘−00100₂. (20) The inner columns of the solution of (20) have Toeplitz structure described by the matrix tridiag(2/5,6/5,2/5) related to the generating function 𝑚(𝑥)=6/5+(4/5)cos(𝑥). Preconditioning with the SPAI matrix results in the smoothing factor 0.250. Indeed, this is a degradation compared to the smoothing factor 0.0588 for the optimal-derived smoother 𝑀3_opt but still an improvement compared to Jacobi. Note that MSPAI minimizes the 2-norm of a combination of conditions while 𝑀3_opt is derived by minimizing the 1-norm. Therefore, to derive efficient smoothing factors by MSPAI, it is important to find a good combination and weighting of probing conditions.

Let us consider a numerical experiment for the matrix 𝐴₁ of order 𝑛=10³. Using an approximate inverse preconditioner satisfying the individual probing conditions, we achieve highly reduced smoothing factors in comparison to SPAI and Jacobi. Table 1 shows that depending on the probing mask and weight 𝜌, reductions down to 𝜇=0.075 are possible for 𝑆3_𝐿3 weighed with 𝜌=10². Moreover, the choice of the subspace weight is stable; that is, increasing values lead to a saturation of the achievable smoothing factors.

Table 1 

Smoothing factors by using 𝑀3_opt (16), SPAI, and MSPAI with local probing conditions for the constant coefficient system 𝐴₁ of order 𝑛=10³.

2.1.3. Global Probing Vectors

Considering the conditions 𝑆3_𝐿2 and 𝑆3_𝐿3 with isotropic masks it is possible to derive conditions with low- and high-frequency global probing vectors 𝑒_𝑆∶=(1,1,…,1)^𝑇 and 𝑒_𝑁1∶=(1,−1,1,−1,…)^𝑇 for the generalized Frobenius norm minimization (5) of MSPAI min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑆1𝑀−2𝑒^𝑇_𝑆2,min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁1𝑀−√2𝑒^𝑇_𝑁12, (21) respectively. In the following, we use a different approach and derive global probing conditions with respect to the optimal derived preconditioner 𝑀3_opt (16). This can also be applied to general probing vectors like 𝑒_𝑁2∶=(1,0,−1,0,1,…)^𝑇 and 𝑒_𝑁3∶=(0,1,0,−1,0,1,…)^𝑇, representing additional high-frequency subspaces. For a given probing vector 𝑒, we observe that 𝑒^𝑇𝑀3_opt=𝛼𝑒^𝑇 with 𝛼∈ℝ. As we want to find the smoother 𝑀, based on MSPAI with the same action on 𝑒^𝑇, we define probing conditions 𝑒^𝑇𝑀≈𝛼𝑒^𝑇 and obtain 𝑆3_𝐺1∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑆𝑀−24𝑒17^𝑇𝑆₂,𝑆3_𝐺2∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁18𝑀−𝑒17^𝑇𝑁1₂,𝑆3_𝐺3∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁2𝑀−16𝑒17^𝑇𝑁2₂,𝑆3_𝐺4∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁3𝑀−16𝑒17^𝑇𝑁3₂. (22)

Figure 1(a) shows the impact of using these global probing conditions for matrix 𝐴₁. Inducing MSPAI to satisfy 𝑆3_𝐺1 leads to a strong reduction of the smoothing factor in comparison to SPAI, close to the optimal value of 𝑀3_opt. A combination of global conditions can lead to an efficient smoother as well. Again the smoothing factor stays stable for increasing values of 𝜌.

(a) MSPAI with global probing conditions for 𝐴1

(b) MSPAI using global probing vectors with action on 𝐵1

Figure 1 

Smoothing factor 𝜇 against subspace weight 𝜌. (a) considers MSPAI using global probing conditions for the constant coefficient system 𝐴1 of size 𝑛=103. (b) shows MSPAI using global probing vectors with action on 𝐵1 with varying coefficients of size 𝑛=103. MSPAI is compared to Jacobi, SPAI, and 𝑀3opt (16).

In a similar way, it is possible to derive global probing conditions with action on 𝐴, that is, min_{𝒫(𝑀)=𝒫}‖𝑒^𝑇𝐴𝑀−𝛼𝑒^𝑇‖₂. For probing vectors, for example, 𝑒_𝑁1,𝑒_𝑁2, and 𝑒_𝑁3, the approximation 𝑒^𝑇𝐴≈const⋅𝑒^𝑇 holds, and therefore 𝑒^𝑇𝐴𝑀≈const⋅𝑒^𝑇𝑀 up to some boundary perturbations. Note that for an approximate inverse preconditioner which should nearly be exact on 𝑒 the condition 𝑒^𝑇𝐴𝑀≈𝑒^𝑇 should be satisfied.

As a generalization, we are interested in the impact of using global probing conditions when considering the 1D tridiagonal matrix 𝐵₁ with varying coefficients and 𝑘th row defined by 𝐵_1𝑘,∶∶=0,…,0,−𝑏_𝑘−1,𝑏_𝑘−1+𝑏_𝑘,−𝑏_𝑘,−1,0,…,0with𝑏∶=2−𝑗2𝑛_{𝑗=0,…,𝑛}. (23) We denote the exact conditions min_{𝒫(𝑀)=𝒫}‖𝑒^𝑇𝐵₁𝑀−(16/17)𝑒^𝑇‖₂ for the probing vectors 𝑒_𝑁1, 𝑒_𝑁2, and 𝑒_𝑁3 as 𝑆3_𝐺5, 𝑆3_𝐺6, and 𝑆3_𝐺7. The approximate conditions for 𝛼=1 are indicated by 𝑆3_𝐺5, 𝑆3_𝐺6, and 𝑆3_𝐺7.

Figure 1(b) shows that both the exact and approximate probing on the high-frequency part of the system lead to significant improvement compared to SPAI and Jacobi, similar to the 1D problem 𝐴₁.

2.2. A 2D Model Problem

We consider the block tridiagonal matrix 𝐴₂𝐴∶=¹⊗𝐼+𝐼⊗𝐴₁21=blocktridiag0,−4−1,041,1,−410,−4,0, (24) with 𝐴₂∈ℝ^𝑛2×𝑛2 related to the generating function 𝑎₂(𝑥)𝑦𝑥,𝑦=1−cos+cos2. (25) As before the smoothing corresponds to the rectangle 𝐼₂∶={(𝑥,𝑦)∣𝑥∈[𝜋/2,𝜋]∧𝑦∈[0,𝜋]}. Hence, the solution of min_𝜔max_{𝑥,𝑦∈𝐼2}|1−𝜔𝑎₂(𝑥,𝑦)| yields the optimal Jacobi smoother with smoothing factor 3/5=0.6 for 𝜔=4/5.

2.2.1. Analytical Derivation of the Optimal Smoother

As approximate inverse smoother, we use the trigonometric polynomial 𝑚₂(𝑥,𝑦)=𝑎+2𝑏(cos(𝑥)+cos(𝑦)). Via the minimization 1−𝑚₂(𝑥,𝑦)𝑎₂(𝑥,𝑦) over 𝐼₂, the solution for the optimal smoothing preconditioner is given by𝑀5_opt8∶=blocktridiag0,841,0,4148,8418410,,41,0 (26) with smoothing factor 9/41=0.2195.

2.2.2. Individual Probing Masks

Analogously to the 1D case, the minimization can be derived by considering min_𝑎,𝑏ℎ1−𝑎+2𝑏ℎ𝑥,𝑦1−𝑥,𝑦2,withcos(𝑥)+cos𝑦=∶ℎ𝑥,𝑦∈[−2,1]. (27) The corners of 𝐼₂ result in the minimization conditionsmin_𝑎,𝑏|,𝑎|2𝑎−8𝑏−12,𝑎+𝑏−12𝑎+𝑏+²16𝑏−1. (28) To get rid of the quadratic condition, we make use of the fact that the linear conditions in (28) have the same absolute value for the optimal 𝑎 and 𝑏. We obtain 𝑎=6𝑏 and the quadratic term can be replaced either by 6𝑏𝑎/𝑏=6𝑎 or (6𝑏)²/𝑏=36𝑏. Consequently, we end up with two additional linear conditionsmin_𝑎,𝑏7𝑎8,𝑎+𝑏−12+13𝑏4−1. (29)

Once again, we can see the minimizations (28) and (29) on 𝑎 and 𝑏 as conditions for the entries of the smoother directly. A column 𝑘 of 𝑀 is now described by the five degrees of freedom (0,𝑚_{𝑘−𝑛,𝑘},0|𝑚_{𝑘−1,𝑘},𝑚_𝑘,𝑘,𝑚_𝑘+1,𝑘|0,𝑚_{𝑘+𝑛,𝑘},0)^𝑇=∶𝑚_𝑘. Thus, we gain the following local probing conditions for the 5-point stencil 𝐴₂: 𝑆5_𝐿1𝑚∶min_𝑘(0,−2,0|−2,2,−2|0,−2,0)𝑚_𝑘−1,𝑆5_𝐿2𝑚∶min_𝑘10,41,04,12,1410,4,0𝑚_𝑘−1,𝑆5_𝐿3𝑚∶min_𝑘10,41,04,78,1410,4,0𝑚_𝑘−1,𝑆5_𝐿4𝑚∶min_𝑘0,1316,013,1162,13160,13𝑚16,0_𝑘−1. (30) It is possible to combine the conditions 𝑆5_𝐿2 and 𝑆5_𝐿4 to obtain an isotropic mask to which the low-frequency probing vector 𝑒_𝑆 corresponds to. To derive this, we add the diagonal and subdiagonal values of the weighed conditions 𝑟⋅𝑆5_𝐿2 and 𝑠⋅𝑆5_𝐿4 and set them to be equal. The isotropy condition leads to the equation (𝑟+𝑠)/2=𝑟/4+13𝑠/16 with the solution 𝑟=5𝑠/4. We obtain 𝑎=𝑏=(𝑟+𝑠)/2=9𝑠/8 and 𝑟+𝑠=5𝑠/4+𝑠=9𝑠/4 for the right hand side. This leads to the additional individual isotropic probing condition 𝑆5_𝐿5𝑚∶min_𝑘(0,1,0|1,1,1|0,1,0)𝑚_𝑘−2. (31)

Let us consider 𝐴₂ of size 𝑛=1024. Following Figure 2(a) we can see that MSPAI satisfying individual probing conditions reduces the smoothing factor in comparison to SPAI with factor 0.339. Utilizing more degrees of freedom with a combination of probing masks, the smoothing factor can be reduced further towards the optimal value of 𝑀5_opt.

(a) MSPAI using individual probing conditions

(b) MSPAI using global probing conditions

Figure 2 

Smoothing factor 𝜇 against subspace weight 𝜌. MSPAI using individual (a) and global (b) probing conditions compared to Jacobi, SPAI, and 𝑀5opt (26) for the 2D constant-coefficient model problem 𝐴2 of size 𝑛=1024.

2.2.3. Global Probing Vectors

Similar to the 1D model problem, we observe that (𝑒⊗𝑓)^𝑇𝑀5_opt=𝛼(𝑒⊗𝑓)^𝑇 for certain 𝑒,𝑓, and α. Therefore, we define global probing conditions by using (𝑒⊗𝑓)^𝑇𝑀=𝛼(𝑒⊗𝑓)^𝑇. Regarding the 2D case, we use global vectors resulting from the Kronecker products 𝑒_𝑆∶=𝑒_𝑆⊗𝑒_𝑆, 𝑒_𝑁1∶=𝑒_𝑁1⊗𝑒_𝑁1, 𝑒_𝑁2∶=𝑒_𝑁2⊗𝑒_𝑁2, and 𝑒_𝑁3∶=𝑒_𝑁3⊗𝑒_𝑁3. We obtain 𝑆5_𝐺1∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑆𝑀−8041𝑒^𝑇_𝑆2,𝑆5_𝐺2∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁1𝑀−1641𝑒^𝑇_𝑁12,𝑆5_𝐺3∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁2𝑀−4841𝑒^𝑇_𝑁22,𝑆5_𝐺4∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁3𝑀−4841𝑒^𝑇_𝑁32. (32)

Like MSPAI with individual probing masks an approximation on the subspace spanned by the global probing vectors reduces the smoothing factor compared to SPAI and Jacobi (see Figure 2(b)). The global condition 𝑆5_𝐺1 leads to a stable smoothing factor of 0.252. Similar to the 1D model problem it is also possible to derive global probing conditions with action on 𝐴₂, not considered here.

2.2.4. 9-Point Stencil of the 2D Model Problem

As a second 2D example, we consider the 9-point stencil of 𝐴₂, which is the block tridiagonal matrix 𝐴₃−1∶=blocktridiag81,−81,−8−181,1,−8−181,−81,−8, (33) related to the generating function𝑎₃ℎ𝑥,𝑦=1−𝑥,𝑦4, (34) with ℎ(𝑥,𝑦)∶=cos(𝑥)+cos(𝑦)+cos(𝑥−𝑦)+cos(𝑥+𝑦).

We use the pattern of 𝑎₃(𝑥,𝑦) for our approximate inverse smoother 𝑚₃𝑥,𝑦=𝑎+2𝑏cos(𝑥)+cos𝑦+cos𝑥−𝑦+cos𝑥+𝑦. (35) With ℎ(𝑥,𝑦)∈[−1,2], this yields the minimization problemmin_𝑎,𝑏ℎ1−𝑎+2𝑏ℎ𝑥,𝑦1−𝑥,𝑦4. (36) The optimal solution is given by (𝑎,𝑏)=(160/153,16/153) which leads to𝑀9_opt=blocktridiag16,15316,1531615316,153160,1531615316,15316,15316,153 (37) with smoothing factor 1/17=0.0588.

We consider the boundary values ℎ(𝑥,𝑦)=−1 and ℎ(𝑥,𝑦)=2 to get linear conditions. A third, quadratic condition can be deduced by setting the derivative of (36) to zero and inserting its solution ℎ^′(𝑥,𝑦)=2−(𝑎/4𝑏) into (36). We obtain the following minimization conditions at high-frequency values: min_𝑎,𝑏3(𝑎−4𝑏)2,3()−1𝑎+2𝑏4,𝑎−12𝑎+2𝑏+².32𝑏−1 (38) The direct translation into linear probing masks and their combination reveals the individual probing conditions 𝑆9_𝐿1𝑚∶min_𝑘−343,−43,−4−34,323,−4−343,−43,−4⋅𝑚_𝑘−1,𝑆9_𝐿2𝑚∶min_𝑘3,316,316316,3164,3316,316,316𝑚16^𝑘−1,𝑆9_𝐿3𝑚∶min_𝑘14,14,1414,13,116414,14,14𝑚_𝑘−1,𝑆9_𝐿4𝑚∶min_𝑘41,6441,64416441,1642,416441,6441,6441𝑚64^𝑘−1,𝑆9_𝐿5𝑚∶min_𝑘(1,1,1|1,1,1|1,1,1)𝑚_𝑘−20.11 (39) In Table 2, we can see that in comparison to SPAI, all individual probing conditions lead to a reduced smoothing factor and are stable for increasing values of 𝜌. It is almost feasible to reach the optimal value 0.0588. Again, it is possible to derive global probing conditions, for example, with action on 𝐴₃, by using Kronecker products of the 1D probing vectors. By observing (𝑒⊗𝑓)^𝑇𝐴𝑀9_opt=𝛼(𝑒⊗𝑓)^𝑇 similar to Section 2.2.3, we can define conditions with (𝑒⊗𝑓)^𝑇𝐴𝑀=𝛼(𝑒⊗𝑓)^𝑇 for some 𝑒,𝑓, and α.

Table 2 

Smoothing factors by using 𝑀9_opt (37), SPAI, and MSPAI with local probing conditions for the 2D 9-point stencil 𝐴₃ of order 𝑛=1024.

3. MSPAI in Regularization

An optimal preconditioner for iterative regularization methods should treat the large eigenvalues and have no effect on the smaller eigenvalues not amplifying the noise. Since SPAI is an efficient smoother and satisfies the first condition, we propose MSPAI probing as regularizing preconditioner for general 𝐻 to suppress a reconstruction on the noise space. Following [20], such a preconditioner 𝑀 should have the following properties: (i)𝑀≈|𝐻|⁻¹ on the signal subspace with |𝐻|=𝐻^𝑇𝐻, and (ii)𝑀≈𝐼 or 𝑀≈0 on the noise subspace.

For circulant matrices, the eigendecomposition is known, and, therefore, these conditions can be satisfied by manipulating the eigenvalues. Most of the preconditioners make use of properties of structured matrices. For general matrices, this is usually not possible, and we thus use the probing facility of MSPAI in order to derive a different approximation quality on the signal or noise subspace, respectively. (i)For the signal space, we could use the vector 𝑒_𝑆∶=(1,1,…,1)^𝑇 representing smooth components, and therefore, the important part of the signal subspace. In this paper, we expect that the signal subspace is already taken into account by SPAI itself, and thus we omit this probing possibility. (ii)For the noise space, we use 𝑒_𝑁1∶=(1,−1,1,−1,…)^𝑇, 𝑒_𝑁2∶=(1,0,−1,0,1,…)^𝑇, and 𝑒_𝑁3∶=(0,1,0,−1,0,…)^𝑇 as typical vectors related to fast oscillations. In case of using probing masks, 𝑠_𝑘̂=(−1,1,−1) and 𝑠_𝑘̂=(1,0,−1) are related to 𝑒_𝑁1 and 𝑒_𝑁2,  𝑒_𝑁3, respectively.

For higher dimensional problems, probing vectors typically result from a Kronecker product of 1D probing vectors. The conditions in MSPAI are given by 𝐴𝑀≈𝐼 in order to derive a good preconditioner and fast convergence on the signal subspace and 𝜌𝑒^𝑇_𝑁𝑀≈0 with 𝜌>0 for the noise subspace in order to avoid a deterioration of the reconstruction by the preconditioner.

3.1. A 1D Model Problem

Let us consider the matrix 𝐻₁∶=tridiag(1/2,1,1/2) related to the symbol ℎ₁(𝑥)𝑒=1+^𝑖𝑥+𝑒^−𝑖𝑥2=1+cos(𝑥). (40) As approximate inverse preconditioner, we choose the trigonometric polynomial of symbol 𝑚₄(𝑥)=𝑎+2𝑏cos(𝑥) related to a Toeplitz matrix 𝑀=tridiag(𝑏,𝑎,𝑏). The entries 𝑎 and 𝑏 should be determined such that the preconditioner 𝑀 acts both nearly as the inverse on the signal subspace and as zero on the noise subspace. For 𝑥∈{0,𝜋/2,𝜋}, this leads to the minimization conditionsmin_𝑎,𝑏𝑚₄(0)ℎ₁(0)−1,𝑚₄𝜋2ℎ₁𝜋2−1,𝜌𝑚₄(𝜋)=min_𝑎,𝑏|2(𝑎+2𝑏)−1|,|𝑎−1|,𝜌|𝑎−2𝑏|. (41) By using the factor 𝜌, we introduce a weighting between the signal and the noise conditions. We obtain the optimal solution of (41) by solving the equality (2𝑎+4𝑏−1)=(𝑎−1)=−𝜌(𝑎−2𝑏). (42) Thus, we are able to derive the optimal regularizing preconditioner𝑀3_𝜌∶=tridiag𝜌−0.5,2+5𝜌2+2𝜌,2+5𝜌𝜌−0.5.2+5𝜌 (43) In the limiting case, for only preconditioning on the signal space (𝜌=0), we get 𝑀3₀=tridiag(−1/4,1,−1/4), and for only regularizing on the noise subspace (𝜌=∞) we get 𝑀3_∞=tridiag(1/5,2/5,1/5)=(2/5)𝐻₁=(2/5)𝐻^𝑇₁. Therefore, the preconditioner for 𝜌=∞ is equivalent to the normal equations. With the parameter 𝜌, we can choose between preconditioning on the signal subspace or suppression of noise.

Again, the conditions (41) can be seen as conditions for the entries of the preconditioner directly, where a column 𝑘 is described by the three degrees of freedom 𝑚_𝑘∶=(𝑚_{𝑘−1,𝑘},𝑚_𝑘,𝑘,𝑚_𝑘+1,𝑘)^𝑇. Both the first and the second condition are covered by SPAI and are not used in MSPAI, as we can think of SPAI approximation as acting mainly on the signal subspace. The translation of the last condition of (41) into an individual probing condition with probing mask 𝑠_𝑘̂ yields 𝑅3_𝐿1𝑚∶min_𝑘(−1,1,−1)𝑚_𝑘. (44) It is possible to transfer this individual regularization condition to the global condition 𝑒^𝑇𝑀≈0 with probing vector 𝑒_𝑁1=(1,−1,1,−1,…)^𝑇, representing high-frequency noisy components. Thus, we force 𝑀 to act as approximately zero on the noise subspace. We refer to the high-frequency global probing condition as 𝑅3_𝐺1∶min_{𝒫(𝑀)=𝒫}𝑒^𝑇_𝑁1𝑀₂. (45)

As a first numerical example, we examine the given 1D blur operator 𝐻₁^2𝑘 with 𝑘>0 for a given vector 𝑥 representing the original data. We use the MSPAI regularization property to reconstruct two different signals perturbed by random noise. By using powers of 𝐻, we make the problem more ill-conditioned and thus enforce the difference between the reconstruction results of all methods. In order to be able to use CG as iterative method, we ensure to have a symmetric and positive definite (spd) preconditioner via corresponding powers of 𝑀𝑀^𝑇 after the computation of 𝑀; that is, the preconditioner is (𝑀𝑀^𝑇)^𝑘. We implement the CG algorithm without a stopping criterion to iterate to the maximum number of specified iterations, that is, to observe the semiconvergent behavior. We keep track of the reconstruction error between the original signal 𝑥 and the reconstruction 𝑥 as well as of the residual in each iteration. If not mentioned otherwise, we refer to CG as using CG without preconditioner and to PCG as CG using a preconditioner, for example, MSPAI with properties on a certain subspace. Notice that in general, an appropriate stopping criterion should be employed [26] when using an iterative regularization method.

Let us consider in a first example the ill-posed problem with operator 𝐻⁸₁ and original data 𝑥₁ of size 𝑛=10³ affected by random noise of order 0.01%. Here, 𝑥₁ denotes the smooth signal 𝑥₁∶=sin4𝜋𝑗𝑛_{𝑗=1,…,𝑛}. (46) The measurement of the error ‖𝑥₁−𝑥₁‖₂ within each iteration shows that for increasing values of 𝜌 the reconstruction is stable, and the region of optimal reconstruction quality is broader and smoother compared to the unpreconditioned CG (see Figure 3). Moreover, by putting more weight to the regularization property of 𝑀3_𝜌 it is possible to approximate the distribution of CG using normal equations. However, we achieve almost similar reconstruction quality in fewer iterations. Refer to Table 3 for details. Hence, with this family of preconditioners, we can steer the iteration to faster convergence (small 𝜌) or slightly better reconstruction quality (𝜌>1).

Table 3 

Optimal reconstruction error ‖𝑥₁−𝑥₁‖₂ for the 1D operator 𝐻⁸₁ and original data 𝑥₁ of size 𝑛=10³ affected by random noise of order 0.01%.

Figure 3 

Reconstruction error against the iteration for the 1D operator 𝐻81 and original data 𝑥1 of size 𝑛=103 affected by random noise of order 0.01%. CG is compared to PCG using the optimal preconditioner 𝑀3𝜌 (43) for different weights 𝜌.

In our second 1D example, we observe the behavior for the blur operator 𝐻⁴₁ and a signal which has strongly increasing entries near the boundaries and an additional nonanalytic point in the middle 𝑥₂1∶=𝑗−1+1𝑛+1−𝑗𝑛/2+1/4−𝑗+sin4𝑗𝜋𝑛_{𝑗=1,…,𝑛}. (47) The problem has size 𝑛=10³ and is perturbed by random noise of order 0.1%. Having a closer look at the error between the original and the reconstructed signal vector shows that the reconstruction is very good for interior components, different from 𝑛/2, but deteriorates near the boundary. Therefore, we introduce some correction at the boundary by changing the components (𝑒_𝑁1)₁=𝜏 and (𝑒_𝑁1)_𝑛=±𝜏, in which 𝜏 is a weight factor of heuristic choice, and in most cases 𝜏∈[0,1]. Additionally, we correct our probing subspace in the middle by (𝑒_𝑁1)_𝑛/2=0 and (𝑒_𝑁1)_𝑛/2±1=1/2(𝑒_𝑁1)_𝑛/2±1 again with heuristic weight 1/2. Following Figure 4, we can see the typical CG behavior of reaching the optimal value after a few iterations but afterwards deteriorating the reconstruction very fast by reducing the error relative to the noise subspace. Using a MSPAI preconditioner the optimum is reached after more iterations but in a stable and smooth manner with slightly smaller value. Table 4 shows the optimal reached reconstruction error with its corresponding number of iterations. Compared to the unpreconditioned CG, the reconstruction 𝑥₂ is by a factor of 2.86 times more accurate when using boundary corrections and 4.72 times when using both corrections within MSPAI.

Table 4 

Optimal reconstruction error ‖𝑥₂−𝑥₂‖₂ for the 1D blur operator 𝐻⁴₁ and the original data 𝑥₂. The problem has size 𝑛=10³ and is affected by random noise of order 0.1%. 𝜌=10².

Figure 4 

Reconstruction error against the iteration for the 1D blur operator 𝐻41 and the original data 𝑥2 of size 𝑛=103 affected by random noise of order 0.1%. CG is compared to CG for the normal equations and to PCG using MSPAI satisfying 𝑅3𝐺1 with 𝜌=102.

Likewise, it is possible to use subspace corrections within mask probing. We change the individual probing condition 𝑅3_𝐿1 for the first and last column of the preconditioner, in order to take into account the missing value −1 which lies outside the vector. This lost information can be incorporated by using the probing mask (1−𝜏,−1+(−1+𝜏))=(1−𝜏,𝜏−2) for the first and the last column. Note that similarly, interior discontinuities can also be treated by modifying the related masks. Furthermore, also for 𝑀3_𝜌 it is possible to build in similar corrections, for example, by weighting the nondiagonal entries. Various experiments revealed that using such corrections at discontinuities of the data vector 𝑥 within the probing subspaces lead to better reconstruction, that is, smaller reconstruction errors compared to the unpreconditioned case. We will use similar techniques in Section 3.3.

As a more general problem, we consider the reconstruction of 𝑥₁ based on the blur operator 𝐵₂ of size 𝑛=10³ affected by random noise of order 1% and 0.1%, respectively. 𝐵₂ denotes the tridiagonal system with varying coefficients whose 𝑘th row is given by 𝐵_2𝑘,∶∶=0,…,0,𝑤_𝑘−1,𝑤_𝑘−1+𝑤_𝑘,𝑤_𝑘,0,…,0,for𝑤∶=2𝑗²𝑛²+1_{𝑗=0,…,𝑛}, (48) with 𝑘∈{1,…,𝑛}. Following Figure 5, CG using normal equations yields both better reconstruction after slightly more iterations and smoother convergence compared to the unpreconditioned CG. In case of applying MSPAI satisfying the global probing condition 𝑅3_𝐺1 to CG(𝐵₂), we are able to reconstruct 𝑥₁ with smaller error and in a smooth way. These positive effects can be enforced by using it for CG with normal equations. The optimal reconstruction is achieved in fewer iterations, is much better, and the convergence curve much smoother.

(a) Noise level 0.1%

(b) Noise level 1%

Figure 5 

Reconstruction error against the iteration for the 1D blur operator 𝐵2 and the original data 𝑥1 of size 𝑛=103 affected by random noise of order 0.1% (a) and 1% (b), respectively. CG is compared to CG using the normal equations, PCG for 𝐵2 and 𝐵𝑇2𝐵2 using MSPAI with 𝑅3𝐺1.

3.2. A 2D Model Problem

In 2D, we consider the matrix𝐻₂𝐻∶=¹⊗𝐼+𝐼⊗𝐻₁21=blocktridiag0,41,041,1,410,4,0, (49) with its corresponding symbol ℎ₂(𝑥,𝑦)=1+cos(𝑥)/2+cos(𝑦)/2 and preconditioner 𝑚₅(𝑥,𝑦)=𝑎+2𝑏(cos(𝑥)+cos(𝑦)). With cos(𝑥)+cos(𝑦)=∶ℎ(𝑥,𝑦)∈[−1,2]=∶𝐼₃, we have to minimizemin_𝑎,𝑏max_{𝑥,𝑦∈𝐼3}ℎ1−𝑎+2𝑏ℎ𝑥,𝑦1−𝑥,𝑦2, (50) relative to the signal subspace and min_𝑎,𝑏𝜌|𝑎−4𝑏| for ℎ(𝑥,𝑦)=−2 relative to the noise subspace. This yields the optimal block tridiagonal preconditioner𝑀5_𝜌3=blocktridiag0,2𝜌−23,02𝜌−23,9+8𝜌,2𝜌−230,2𝜌−2,0. (51)

The individual probing masks for a reduced inner column of the preconditioner 𝑚_𝑘∶=(0,𝑚_{𝑘−𝑛,𝑘},0|𝑚_{𝑘−1,𝑘},𝑚_𝑘,𝑘,𝑚_𝑘+1,𝑘|0,𝑚_{𝑘+𝑛,𝑘},0)^𝑇 as well as the high-frequency probing vector for the 2D case can be derived via Kronecker products of the 1D probing vectors. Thus, we obtain the individual and global probing condition 𝑅5_𝐿1𝑚∶min_𝑘((−1,1,−1)⊗(−1,1,−1))𝑚_𝑘,𝑅5_𝐺1∶min_{𝒫(𝑀)=𝒫}𝑒_𝑁1⊗𝑒_𝑁1^𝑇𝑀₂, (52) respectively.

We observe the behavior for the 2D problem with blur operator 𝐻⁴₂ and original data 𝑥₁ of size 𝑛=50² affected by random noise of order 0.1%. Again, for 𝑀5_𝜌, with increasing values of 𝜌, we approximate the convergence of CG using normal equations (refer to Figure 6). Applying MSPAI with 𝑅5_𝐺1 to PCG it is almost possible to reach the value of CG(𝐻⁴₂^𝑇𝐻⁴₂) but in less iterations. Higher values of 𝜌 lead to smooth and broad convergence curves similar to 𝑀5_𝜌. CG has its optimal value after 9 iterations with error 2.219⋅10⁻³, MSPAI using 𝑅5_𝐺1, 𝜌=1 after 12 iterations with value 1.496⋅10⁻³, and CG using normal equations reaches the error 1.282⋅10⁻³ after 44 iterations.

(a) PCG using the optimal preconditioner 𝑀5𝜌(51)

(b) PCG using MSPAI satisfying 𝜌𝑅5𝐺1

Figure 6 

Reconstruction error against the iteration for the 2D problem with blur operator 𝐻42 and original data 𝑥1 of size 𝑛=502 affected by random noise of order 0.1%. CG is compared to CG for the normal equations and to 𝑀5𝜌 (a) and MSPAI satisfying 𝜌𝑅5𝐺1 (b), respectively.

3.3. Examples from the Regularization Toolbox

We are interested in the impact of using sparse approximate inverse preconditioners on some problems of the MATLAB package Regularization Tools Version 4.1 for analysis and solution of discrete ill-posed problems, developed by Hansen [24].

In our first example, we consider the deriv2 example which is a discretization of a first kind Fredholm integral equation. We choose deriv2(n,case) with case=2. Our problem has size 𝑛=2⋅10³ and is affected with random noise of order 0.001% and 0.01%, respectively. Note that the system matrix 𝐾 is dense. We apply both CG and PCG to 𝐾^𝑇𝐾 and force the MSPAI to act as approximately zero on the noise subspace by using the high-frequency probing conditions 𝑅3_𝐿1 and 𝑅3_𝐿2∶min𝑚_𝑘|(1,0,−1)𝑚_𝑘| simultaneously. We weigh the subspace with 𝜌=10³and apply the symmetric preconditioner 𝑀+𝑀^𝑇 to the normal equations. To avoid deterioration at the boundary, we adjust 𝑀 by resetting the values 𝑀_1,1=𝑀_2,2, 𝑀_2,1=𝑀_3,2, 𝑀_𝑛,𝑛=𝑀_{𝑛−1,𝑛−1}, and 𝑀_{𝑛−1,𝑛}=𝑀_{𝑛−2,𝑛−1}.

Following Figure 7, we are able to achieve better reconstruction 𝑔 of the original data 𝑔 and in fewer iterations when applying MSPAI. For other noise levels and for deriv(n,1), we observed similar behavior.

(a) Noise level 0.001%

(b) Noise level 0.01%

Figure 7 

Reconstruction error against the iteration for the deriv2 problem of [24] invoked with deriv2(2000,2) and affected by random noise of order 0.001% (a) and 0.01% (b), respectively. CG is compared to PCG using MSPAI with 𝑅3𝐿1, 𝑅3𝐿2, and 𝜌=103 for the normal equations. MSPAI is symmetrized via 𝑀+𝑀𝑇.

Let us focus on the blur [24] test problem as a second example which is deblurring images degraded by atmospheric turbulence blur. The matrix 𝐺 is an 𝑛²-by-𝑛² symmetric, doubly block Toeplitz matrix that models blurring of an 𝑛-by-𝑛 image by a Gaussian point spread function. The parameter 𝜎 controls the width of 𝐺 and thus the amount of smoothing and ill-posedness. 𝐺=𝐴₁⊗𝐴₁ is symmetric block banded and possibly positive definite depending on 𝑛 and 𝜎. We choose the problem to be of size 𝑛=150 with bandwidth 4 and 𝜎=1, that is, we invoke blur(150,4,1) and 𝐺 is of size 150²×150². The original data vector is denoted by 𝑝.

We compare the unpreconditioned CG to PCG both for 𝐺 and for the normal equations 𝐺^𝑇𝐺. In case of preconditioning 𝐺, MSPAI is symmetrized via 𝑀+𝑀^𝑇 and for 𝐺^𝑇𝐺 the preconditioner 𝑀^𝑇𝑀 is applied. For MSPAI, we impose the blocktridiagonal pattern 𝑚_𝑘∶=(0,𝑚_{𝑘−𝑛,𝑘},0|𝑚_{𝑘−1,𝑘},𝑚_𝑘,𝑘,𝑚_𝑘+1,𝑘|0,𝑚_{𝑘+𝑛,𝑘},0)^𝑇. In view of the structure of 𝐺, we build the high-frequency subspace by Kronecker products of oscillatory probing vectors in the regularizing global conditions 𝑅5_𝐺1 and the new ones 𝑅5_𝐺2∶min_{𝒫(𝑀)=𝒫}‖(𝑒_𝑁1⊗𝑒_𝑁2)^𝑇𝑀‖₂ and 𝑅5_𝐺3∶min_{𝒫(𝑀)=𝒫}‖(𝑒_𝑁2⊗𝑒_𝑁1)^𝑇𝑀‖₂, all weighed with 𝜌=1.

Following Figure 8 and Table 5, we obtain better reconstruction 𝑝 when applying MSPAI in contrast to the unpreconditioned CG for 𝐺 or 𝐺^𝑇𝐺. We observed similar behavior for other values of 𝑏𝑎𝑛𝑑 and 𝜎.

Table 5 

Optimal reconstruction error ‖𝑝−𝑝‖₂ at given iteration for the blur problem of [24] invoked with blur(150,4,1) for noise of order 0.01%, 0.1%, 1%, and 10%. Subspace weight is 𝜌=1.

(a) Noise level 0.01%

(b) Noise level 1%

Figure 8 

Reconstruction error against the iteration for the blur problem of [24] invoked with blur(150,4,1) and affected by random noise of order 0.01% (a) and 1% (b), respectively. CG is compared to CG using normal equations, PCG using MSPAI with 𝑅5𝐺1, 𝑅5𝐺2, 𝑅5𝐺3, and PCG using MSPAI for normal equations with 𝜌=1. We use 𝑀+𝑀𝑇 as preconditioner within PCG(𝐺) and 𝑀𝑇𝑀 for PCG(𝐺𝑇𝐺).

As a last example we reduce the 2D blur example of Hansen to the 1D case. For the blur operator we take the 1D analogon of 𝐺 and reduce the 2D right hand side 𝑝 to appropriate size. The consideration of example 𝐻⁴₁ with original data 𝑥₂ (Table 4 and Figure 4) shows that the preconditioner should take also into account the behavior of the original or blurred data vector. Smoothing, for example, with 𝑀=tridiag(1/2,1,1/2) makes sense to remove noisy components only as long as the data is continuous. At discontinuities, smoothing would cause additional errors. Therefore, we use a modified tridiagonal smoothing preconditioner with 𝑗th row (0,…,0,𝑟_𝑗−1,1,𝑟_𝑗,0,…,0). Here, 𝑟_𝑗≈1/2 near continuous components 𝑥_𝑗, but 𝑟_𝑗≈0 near discontinuities. In case that we know the original data 𝑥 we define 𝑀_𝑟 with 𝑗th row 𝑀_𝑟𝑗,∶∶=0,…,0,𝑟_𝑗−1,1,𝑟_𝑗,1,0,…,0for𝑟∶=2⋅11+𝜌𝑥_𝑗−𝑥_{𝑗+1𝑗=1,…,𝑛−1}. (53) Otherwise we define the preconditioner 𝑀_̃𝑟 with 𝑗th row 𝑀_{̃𝑟𝑗,∶}∶=0,…,0,𝑟̃_𝑗−1,1,𝑟_𝑗̃,̃1,0,…,0via𝑟∶=2⋅11+𝜌𝑏_𝑗−𝑏_𝑗+1^𝑘_{𝑗=1,…,𝑛−1}. (54) The parameters 𝜌 and 𝑘 have to be chosen in such a way that discontinuities are revealed as good as possible.

Following Figure 9, the reconstruction is strongly improved if the preconditioner is adjusted relative to the discontinuities of 𝑥. Also, using only the observed data improves the reconstruction. Therefore, we could also consider an iterative process, where a first approximation 𝑥₁ on 𝑥 is used to define the tridiagonal preconditioner 𝑀₁ delivering an improved approximation 𝑥₂ which again gives a new preconditioner 𝑀₂, and so on.

(a) 𝐺 with 𝜎=1

(b) 𝐺 with 𝜎=10

Figure 9 

Reconstruction error against iterations for the 1D blur problem of size 𝑛=103 and band=4 affected by random noise of order 0.1%. 𝐺 has 𝜎=1 (a) and 𝜎=10 (b), respectively. CG is compared to PCG using 𝑀𝑟 and 𝑀̃𝑟, respectively, for 𝐺 and 𝐺𝑇𝐺.

4. Conclusion

We have considered the derivation of preconditioners with special behavior on certain subspaces. For this purpose, analytic minimization problems in functions have been translated into MSPAI minimizations for vectors based on masks. Such mask-based probing conditions can be different for each column 𝑚_𝑘 of the preconditioner 𝑀 and can be written in the form min𝑚_𝑘|𝑠_𝑘̂𝑚_𝑘−𝑓_𝑘| with reduced vectors 𝑠_𝑘̂, 𝑚_𝑘, and scalar 𝑓_𝑘∈ℝ. Mask probing has the advantage that for each sparse column 𝑚_𝑘 we can use a different sparse probing vector 𝑠_𝑘. Furthermore, we have introduced probing conditions based on probing vectors that are global for the whole matrix 𝑀 in the form min_{𝒫(𝑀)=𝒫}‖𝑒^𝑇𝑀−𝑓^𝑇‖₂ and min_{𝒫(𝑀)=𝒫}‖𝑒^𝑇𝐴𝑀−𝑓^𝑇‖₂, respectively, with 𝑒,𝑓∈ℝ^𝑛. The probing vectors are related to the low-frequency or smooth subspace, represented by 𝑒_𝑆=(1,1,…,1)^𝑇, or to the high-frequency oscillatory subspace, represented, for example, by 𝑒_𝑁1=(1,−1,1,−1,…)^𝑇, 𝑒_𝑁2=(1,0,−1,0,1,…)^𝑇, or 𝑒_𝑁3=(0,1,0,−1,0,…)^𝑇.

For multigrid methods, we have shown that the smoothing property of approximate inverses like SAI or SPAI can be improved significantly by using MSPAI with appropriate probing masks or probing vectors. In special cases, we can analytically determine the optimal approximate inverse smoother and its corresponding smoothing factor. SPAI is far from being optimal in this class, but including convenient individual or global probing conditions, we can nearly reach the optimal smoothing factor with MSPAI. Moreover, an increasing weighting of the subspace leads to stable behavior. Our tests on systems with varying coefficients and for the 2D case demonstrate that even the usage of global probing conditions with action on 𝐴, only satisfying the approximation 𝑒^𝑇𝐴𝑀≈𝑒^𝑇, reduces the smoothing factor in comparison to SPAI and the damped Jacobi method.

A different subspace approach becomes necessary during the recovery process of blurred signals. Here, our main focus is to allow a better and more stable reconstruction of the original signal. Applying a preconditioner in iterative regularization can easily lead to a deterioration of the reconstruction by approximating the inverse also in the noise subspace, or by removing high-frequency components in the original signal. Therefore, preconditioners have to be developed very carefully. We derive approximate inverse preconditioners analytically based on generating functions, or by applying MSPAI with probing masks or probing vectors. These preconditioners allow to incorporate filtering for noise reduction, and they can be adjusted both to the system matrix, for example, the blur operator, and to the data vector 𝑥. So, the deterioration of the reconstruction at discontinuities of 𝑥 can be reduced by modifying the probing conditions relative to the variation of the signal data. We show that these preconditioners can be used for faster convergence or better reconstruction. The application to more general problems is an interesting and important task which will be investigated in the future.