# Smoothing and Regularization with Modified Sparse Approximate Inverses

Authors: T. Huckle; M. Sedlacek
Journal: Journal of Electrical and Computer Engineering
 (2010)

---

## Abstract

Sparse approximate inversesM which satisfy minM∥AM−I∥F
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.


---

## Body

## 1. Introduction

For applying an iterative solution method to an ill-conditioned system of linear equationsAx=b with sparse matrix A∈ℝn×n, x,b∈ℝn, it is often crucial to include an efficient preconditioner. Here, the original problem Ax=b is replaced by the preconditioned system MAx=Mb or Ax=A(My)=b. 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 M as well as the application of the preconditioner on any given vector v should be efficient. Furthermore, an iterative solver applied on AMy=b or MAx=Mb should converge much faster than for Ax=b (e.g., it holds cond(MA)≪cond(A)).The first two conditions can be easily satisfied by using a sparse matrixM as approximation to A-1. Note that the inverse of a sparse system A is nearly dense, but in many cases, the entries of A-1 are rapidly decaying, so most of the entries are very small (see Demko et al. [1]).Sparse approximate inverses can be computed by minimizingAM-I in the Frobenius norm, where I 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 Section1, 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 M≈A-1 to the inverse of a system matrix A∈ℝn×n. For an a priori prescribed sparsity pattern 𝒫, this can be done in a static way by solving(1)min⁡P(M)=P⁡‖AM-I‖F2=∑k=1nmin⁡P(mk)=Pk⁡‖Amk-ek‖22,
with mk the kth column of the preconditioner M and ek the kth column of the identity matrix I. 𝒫(A) denotes the pattern 𝒫 of A and 𝒫k the pattern of the kth 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 n small least squares problems regarding one column can be computed independently of one another.Using both the index set𝒥k, which is implicitly given by 𝒫k and contains the indices j such that mk(j)≠0, and its corresponding so-called shadow ℐk, that is the indices of nonzero rows in A(:,𝒥k), each subproblem in (1) is related to a small matrix Âk:=A(ℐk,𝒥k) if 𝒫k is sparse. We refer to the reduced sparse column vectors as m̂k:=mk(𝒥k) and êk:=ek(ℐk), respectively. The solution of every reduced problem minm̂k∥Âm̂k-êk∥2 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 inmk to improve on an already computed approximation. We assume that by solving (1) for a given index set 𝒥k, we already have determined an optimal solution mk(𝒥k) inducing the sparse vector mk with residual rk. Dynamically, we want to augment new entries in mk and solve (1) for this enlarged index set 𝒥̃k such that we derive a reduction in the norm of the new residual r̃k=A(ℐ̃k,𝒥̃k)mk(𝒥̃k)-ek(ℐ̃k).Following Cosgrove et al. [5] and Grote and Huckle [6], in SPAI, we test one possible new index j∈𝒥new out of a given set of possible new indices 𝒥new to improve mk. Therefore, we can consider the reduced 1D minimization problem(2)min⁡λj⁡‖A(mk+λjej)-ek‖2=min⁡λj⁡‖λjAj+rk‖2.
The solution of (2) is given by (3)λj=-rkTAej‖Aej‖22,
and leads to an improved squared residual norm (4)ρj2=‖rk‖2-(rkTAej)2‖Aej‖22.
For improving mk, we only have to consider indices j in rows of A that are related to nonzero entries in the old residual rk; otherwise, they do not lead to a reduction of the residual norm. Thus, we have to determine those column indices j, which satisfy rkTAej≠0. Let us denote the index set of nonzero entries in rk by ℒ. By 𝒥̃i, we denote the set of new indices that are related to the nonzero elements in the ith row of A, and by 𝒥new:=⋃i∈ℒ𝒥̃i the set of all possible new indices that will lead to a reduction of the residual norm. The one or more newly added indices 𝒥c are chosen to be a subset of 𝒥new that corresponds to the maximal reduction in rk. For the enlarged index set 𝒥k∪𝒥c, we have to update the QR decomposition of the related least squares submatrix and solve the new column mk.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 B on the right hand side in the form min𝒫(M)=𝒫∥AM-B∥F. This approach is useful in connection with some kind of two-level preconditioning: first, compute a standard sparse preconditioner B for A, and then improve this preconditioner by an additional Frobenius norm minimization with target B. From an algorithmic point of view, the minimization with sparse target matrix B, instead of I, introduces no additional difficulties. Simply, 𝒫(M) should be chosen more carefully with respect to A and B.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 formulation(5)min⁡P(M)=P⁡‖CM-B‖F=min⁡P(M)=P⁡‖(C0ρeT)M-(B0ρfT)‖F,
with sparse matrices C0 and B0, is not restricted to special probing subspaces as it allows any choice of e and f. The resulting preconditioner M satisfies both C0M≈B0 and eTM≈fT. We refer to the first n rows of (5), that is, C0M-B0, 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 e,f∈ℝn×k represent the k-dimensional subspace on which the preconditioner should be optimal. Choosing ρ=0, C0=A, and B0=I in (5) leads to the classical SPAI formulation. Setting C0=I and B0=A, we end up with a formulation computing explicit sparse approximations to A. In this case, the derived approximation on A can have a considerably fewer number of nonzeros (nnz) than A but choosing fT=eTA, the preconditioner M will have a similar action on eT as A.Furthermore, it is also possible to include individual probing conditions for each columnmk. As a new approach, we use probing masks defined by sparse row vectors sk, k∈{1,…,n}, containing the same pattern as mk, that is, 𝒫(sk)=𝒫(mk), and add the condition minm̂k|ŝkm̂k-fk| to the Frobenius norm minimization of the kth column. Corresponding to m̂k, ŝk denotes the reduced form of sk. The masks for each column of M can be stored in a sparse rectangular matrix S, whereby the individual probing conditions result from diag(SM)≈f. Compared to MSPAI using global probing vectors e,f, 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 ŝk=(1,0,-1) can be used to enforce a quasisymmetry in the column vector mk such that mk-1,k≈mk+1,k.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)T 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 matrixA and an approximate inverse smoother M, the iteration k is described by x(k+1)=x(k)+M(b-Ax(k)), and thus the error is given by I-MA. 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 P=UT(I-MA)U on the high-frequency subspace U gives the smoothing factor defined as the spectral radius of P. In the constant coefficient case, the smoothing factor can also be described as 1-ma, where a and m are generating functions representing A and M [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 x,y∈[0,π], where the high-frequency domain is given by the difference G between the two squares [0,π]2 and [0,π/2]2 with corners at (6)(0,π2),(π2,π2),(π2,0),(π,0),(π,π),(0,π).
Thus, the smoothing factor is given by (7)μ:=max⁡x,y∈G⁡{|1-m(x,y)a(x,y)|}.
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(8)x⟶blurHx⟶noiseHx+η=b,
where x is the original image, H is the blur operator, η is a vector representing the noise, and b is the observed image. We want to recover x as good as possible and as fast as possible. Because H may be extremely ill-conditioned or even singular, and, because of the presence of noise, (8) cannot be solved directly. Consequently, to solve Hx=b on the signal subspace, a regularization technique has to be applied. One of the classical methods is the Tikhonov regularization [17] which solves (9)min⁡x⁡{‖Hx-b‖22+γ‖b‖22}⟺(HTH+γI)x=HTb,
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 H, 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 H, 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].


## 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 inmk to improve on an already computed approximation. We assume that by solving (1) for a given index set 𝒥k, we already have determined an optimal solution mk(𝒥k) inducing the sparse vector mk with residual rk. Dynamically, we want to augment new entries in mk and solve (1) for this enlarged index set 𝒥̃k such that we derive a reduction in the norm of the new residual r̃k=A(ℐ̃k,𝒥̃k)mk(𝒥̃k)-ek(ℐ̃k).Following Cosgrove et al. [5] and Grote and Huckle [6], in SPAI, we test one possible new index j∈𝒥new out of a given set of possible new indices 𝒥new to improve mk. Therefore, we can consider the reduced 1D minimization problem(2)min⁡λj⁡‖A(mk+λjej)-ek‖2=min⁡λj⁡‖λjAj+rk‖2.
The solution of (2) is given by (3)λj=-rkTAej‖Aej‖22,
and leads to an improved squared residual norm (4)ρj2=‖rk‖2-(rkTAej)2‖Aej‖22.
For improving mk, we only have to consider indices j in rows of A that are related to nonzero entries in the old residual rk; otherwise, they do not lead to a reduction of the residual norm. Thus, we have to determine those column indices j, which satisfy rkTAej≠0. Let us denote the index set of nonzero entries in rk by ℒ. By 𝒥̃i, we denote the set of new indices that are related to the nonzero elements in the ith row of A, and by 𝒥new:=⋃i∈ℒ𝒥̃i the set of all possible new indices that will lead to a reduction of the residual norm. The one or more newly added indices 𝒥c are chosen to be a subset of 𝒥new that corresponds to the maximal reduction in rk. For the enlarged index set 𝒥k∪𝒥c, we have to update the QR decomposition of the related least squares submatrix and solve the new column mk.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 B on the right hand side in the form min𝒫(M)=𝒫∥AM-B∥F. This approach is useful in connection with some kind of two-level preconditioning: first, compute a standard sparse preconditioner B for A, and then improve this preconditioner by an additional Frobenius norm minimization with target B. From an algorithmic point of view, the minimization with sparse target matrix B, instead of I, introduces no additional difficulties. Simply, 𝒫(M) should be chosen more carefully with respect to A and B.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 formulation(5)min⁡P(M)=P⁡‖CM-B‖F=min⁡P(M)=P⁡‖(C0ρeT)M-(B0ρfT)‖F,
with sparse matrices C0 and B0, is not restricted to special probing subspaces as it allows any choice of e and f. The resulting preconditioner M satisfies both C0M≈B0 and eTM≈fT. We refer to the first n rows of (5), that is, C0M-B0, 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 e,f∈ℝn×k represent the k-dimensional subspace on which the preconditioner should be optimal. Choosing ρ=0, C0=A, and B0=I in (5) leads to the classical SPAI formulation. Setting C0=I and B0=A, we end up with a formulation computing explicit sparse approximations to A. In this case, the derived approximation on A can have a considerably fewer number of nonzeros (nnz) than A but choosing fT=eTA, the preconditioner M will have a similar action on eT as A.Furthermore, it is also possible to include individual probing conditions for each columnmk. As a new approach, we use probing masks defined by sparse row vectors sk, k∈{1,…,n}, containing the same pattern as mk, that is, 𝒫(sk)=𝒫(mk), and add the condition minm̂k|ŝkm̂k-fk| to the Frobenius norm minimization of the kth column. Corresponding to m̂k, ŝk denotes the reduced form of sk. The masks for each column of M can be stored in a sparse rectangular matrix S, whereby the individual probing conditions result from diag(SM)≈f. Compared to MSPAI using global probing vectors e,f, 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 ŝk=(1,0,-1) can be used to enforce a quasisymmetry in the column vector mk such that mk-1,k≈mk+1,k.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)T 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 matrixA and an approximate inverse smoother M, the iteration k is described by x(k+1)=x(k)+M(b-Ax(k)), and thus the error is given by I-MA. 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 P=UT(I-MA)U on the high-frequency subspace U gives the smoothing factor defined as the spectral radius of P. In the constant coefficient case, the smoothing factor can also be described as 1-ma, where a and m are generating functions representing A and M [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 x,y∈[0,π], where the high-frequency domain is given by the difference G between the two squares [0,π]2 and [0,π/2]2 with corners at (6)(0,π2),(π2,π2),(π2,0),(π,0),(π,π),(0,π).
Thus, the smoothing factor is given by (7)μ:=max⁡x,y∈G⁡{|1-m(x,y)a(x,y)|}.
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(8)x⟶blurHx⟶noiseHx+η=b,
where x is the original image, H is the blur operator, η is a vector representing the noise, and b is the observed image. We want to recover x as good as possible and as fast as possible. Because H may be extremely ill-conditioned or even singular, and, because of the presence of noise, (8) cannot be solved directly. Consequently, to solve Hx=b on the signal subspace, a regularization technique has to be applied. One of the classical methods is the Tikhonov regularization [17] which solves (9)min⁡x⁡{‖Hx-b‖22+γ‖b‖22}⟺(HTH+γI)x=HTb,
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 H, 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 H, 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 formA1x=b, with A1:=tridiag(-1/2,1,-1/2). The matrix A1 is related to the generating function or symbol(10)a1(x)=1-eix+e-ix2=1-cos⁡⁡(x).
The high-frequency part of A1 is represented by (10) for x∈[π/2,π]=:I1. Hence, the optimal smoothing parameter ω in the Jacobi smoother 1-ωa1(x) is found by solving the problem(11)min⁡ω⁡max⁡x∈I1⁡|1-ωa1(x)|=min⁡ω⁡max⁡x∈I1⁡|s(x)|=ω=2/313.
The solution can be found by replacing (11) with the maximum over the two boundary values minωmax{|s(π/2)|,|s(π)|}.


#### 2.1.1. Analytical Derivation of the Optimal Smoother

As approximate inverse smoother, we choose a trigonometric polynomial of the same degree(12)m1(x)=a+2bcos⁡⁡(x),
and a tridiagonal Toeplitz matrix M=tridiag(b,a,b), respectively. The smoothing condition for (12) can be written as (13)min⁡a,b⁡max⁡x∈I1⁡|1-m1(x)a1(x)|=min⁡a,b⁡max⁡x∈I1⁡|s(x)|,
or using the boundary values of I1 as(14)min⁡a,b⁡max⁡⁡{|1-m1(π2)a1(π2)|,|1-m1(π)a1(π)|,|1-m1(u)a1(u)|(π2)}=min⁡a,b⁡max⁡⁡{|a-1|,|2(a-2b)-1|,|1-m1(u)a1(u)|},
in which u∈I1 is the local extreme point of the function 1-m1(x)a1(x) with derivative equal to zero. We thus obtain for u=1/2-a/4b the quadratic condition(15)s(u)=(a+2b)28b-1,
and the overall solution (a,b)=(16/17,4/17) which leads to(16)M3opt:=tridiag(417,1617,417).
The related smoothing factor for this optimal preconditioner m(x) 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 a and b can be also seen as conditions for the entries of the matrix directly. A column k of M is described by the three values (mk-1,k,mk,k,mk+1,k)T=:m̂k, where the main diagonal entry is related to a and the upper and lower entries to b. Therefore, we can translate the two linear minimization conditions of (14) into individual probing masks ŝk on entries of M directly (17)S3L1:min⁡m̂k⁡|(0,1,0)m̂k-1|,S3L2:min⁡m̂k⁡|(-1,1,-1)m̂k-12|.
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 a=1 into the second condition yields b=1/4. We replace the denominator of the quadratic condition with this value and get the new linear condition (18)min⁡a,b⁡|(a+2b)28⋅1/4-1|=min⁡a,b⁡|(a+2b)2-2|⟶min⁡a,b⁡|a+2b-2|,
which is related to the probing condition with isotropic probing mask (19)S3L3:min⁡m̂k⁡|(1,1,1)m̂k-2|.
We call a mask ŝk to be isotropic if the minimization minm̂k|ŝkm̂k-fk| can be described by one probing vector e∈ℝn with e(𝒥k)=αŝkT, α∈ℝ and k∈{1,…,n}, in the form min𝒫(M)=𝒫∥eTM-fT∥2. Here, for instance, ŝk=(-1,1,-1) of S3L2 is related to e=(1,-1,1,-1,…)T and ŝk=(1,1,1) of S3L3 to e=(1,1,…,1)T.In our following numerical experiments, we are interested in a comparison between MSPAI with probing conditions, the optimal smootherM3opt (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 M is given by(20)min⁡m̂k⁡‖(-0.5001-0.50-0.51-0.50-0.5100-0.5)m̂k-(00100)‖2.
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 m(x)=6/5+(4/5)cos(x). 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 M3opt but still an improvement compared to Jacobi. Note that MSPAI minimizes the 2-norm of a combination of conditions while M3opt 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 matrixA1 of order n=103. 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 S3L3 weighed with ρ=102. 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 usingM3opt (16), SPAI, and MSPAI with local probing conditions for the constant coefficient system A1 of order n=103.
M3optSPAIWeightLocal probing conditionsρρS3L1∧0.7S3L2ρS3L30.05880.2501.00.1340.0832.00.1070.0771010.0950.075102


#### 2.1.3. Global Probing Vectors

Considering the conditionsS3L2 and S3L3 with isotropic masks it is possible to derive conditions with low- and high-frequency global probing vectors eS:=(1,1,…,1)T and eN1:=(1,-1,1,-1,…)T for the generalized Frobenius norm minimization (5) of MSPAI (21)min⁡P(M)=P⁡‖eSTM-12eST‖2,min⁡P(M)=P⁡‖eN1TM-2eN1T‖2,
respectively. In the following, we use a different approach and derive global probing conditions with respect to the optimal derived preconditioner M3opt (16). This can also be applied to general probing vectors like eN2:=(1,0,-1,0,1,…)T and eN3:=(0,1,0,-1,0,1,…)T, representing additional high-frequency subspaces. For a given probing vector e, we observe that eTM3opt=αeT with α∈ℝ. As we want to find the smoother M, based on MSPAI with the same action on eT, we define probing conditions eTM≈αeT and obtain (22)S3G1:min⁡P(M)=P⁡‖eSTM-2417eST‖2,S3G2:min⁡P(M)=P⁡‖eN1TM-817eN1T‖2,S3G3:min⁡P(M)=P⁡‖eN2TM-1617eN2T‖2,S3G4:min⁡P(M)=P⁡‖eN3TM-1617eN3T‖2.Figure1(a) shows the impact of using these global probing conditions for matrix A1. Inducing MSPAI to satisfy S3G1 leads to a strong reduction of the smoothing factor in comparison to SPAI, close to the optimal value of M3opt. A combination of global conditions can lead to an efficient smoother as well. Again the smoothing factor stays stable for increasing values of ρ.Smoothing factorμ against subspace weight ρ. (a) considers MSPAI using global probing conditions for the constant coefficient system A1 of size n=103. (b) shows MSPAI using global probing vectors with action on B1 with varying coefficients of size n=103. MSPAI is compared to Jacobi, SPAI, and M3opt (16).
(a)
MSPAI with global probing conditions forA1(b)
MSPAI using global probing vectors with action onB1In a similar way, it is possible to derive global probing conditions with action onA, that is, min𝒫(M)=𝒫∥eTAM-αeT∥2. For probing vectors, for example, eN1,eN2, and eN3, the approximation eTA≈const·eT holds, and therefore eTAM≈const·eTM up to some boundary perturbations. Note that for an approximate inverse preconditioner which should nearly be exact on e the condition eTAM≈eT should be satisfied.As a generalization, we are interested in the impact of using global probing conditions when considering the 1D tridiagonal matrixB1 with varying coefficients and kth row defined by (23)(B1)k,::=(0,…,0,-bk-1,bk-1+bk,-bk,0,…,0),withb:=(-12-j2n)j=0,…,n.
We denote the exact conditions min𝒫(M)=𝒫∥eTB1M-(16/17)eT∥2 for the probing vectors eN1, eN2, and eN3 as S3G5, S3G6, and S3G7. The approximate conditions for α=1 are indicated by S3G̃5, S3G̃6, and S3G̃7.Figure1(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 A1.


### 2.2. A 2D Model Problem

We consider the block tridiagonal matrix(24)A2:=A1⊗I+I⊗A12=blocktridiag(0,-14,0|-14,1,-14|0,-14,0),
with A2∈ℝn2×n2 related to the generating function (25)a2(x,y)=1-cos⁡⁡(x)+cos⁡⁡(y)2.
As before the smoothing corresponds to the rectangle I2:={(x,y)∣x∈[π/2,π]∧y∈[0,π]}. Hence, the solution of minωmaxx,y∈I2|1-ωa2(x,y)| 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 polynomialm2(x,y)=a+2b(cos(x)+cos(y)). Via the minimization 1-m2(x,y)a2(x,y) over I2, the solution for the optimal smoothing preconditioner is given by(26)M5opt:=blocktridiag(0,841,0|841,4841,841|0,841,0),
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(27)min⁡a,b⁡|1-(a+2bh(x,y))(1-h(x,y)2)|,withcos⁡⁡(x)+cos⁡⁡(y)=:h(x,y)∈[-2,1].
The corners of I2 result in the minimization conditions(28)min⁡a,b⁡{|2a-8b-1|,|a2+b-1|,|a2+b+a216b-1|}.
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 a and b. We obtain a=6b and the quadratic term can be replaced either by 6ba/b=6a or (6b)2/b=36b. Consequently, we end up with two additional linear conditions(29)min⁡a,b⁡{|7a8+b-1|,|a2+13b4-1|}.Once again, we can see the minimizations (28) and (29) on a and b as conditions for the entries of the smoother directly. A column k of M is now described by the five degrees of freedom (0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T=:m̂k. Thus, we gain the following local probing conditions for the 5-point stencil A2: (30)S5L1:min⁡m̂k⁡|(0,-2,0|-2,2,-2|0,-2,0)m̂k-1|,S5L2:min⁡m̂k⁡|(0,14,0|14,12,14|0,14,0)m̂k-1|,S5L3:min⁡m̂k⁡|(0,14,0|14,78,14|0,14,0)m̂k-1|,S5L4:min⁡m̂k⁡|(0,1316,0|1316,12,1316|0,1316,0)m̂k-1|.
It is possible to combine the conditions S5L2 and S5L4 to obtain an isotropic mask to which the low-frequency probing vector eS corresponds to. To derive this, we add the diagonal and subdiagonal values of the weighed conditions r·S5L2 and s·S5L4 and set them to be equal. The isotropy condition leads to the equation (r+s)/2=r/4+13s/16 with the solution r=5s/4. We obtain a=b=(r+s)/2=9s/8 and r+s=5s/4+s=9s/4 for the right hand side. This leads to the additional individual isotropic probing condition (31)S5L5:min⁡m̂k⁡|(0,1,0|1,1,1|0,1,0)m̂k-2|.Let us considerA2 of size n=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 M5opt.Smoothing factorμ against subspace weight ρ. MSPAI using individual (a) and global (b) probing conditions compared to Jacobi, SPAI, and M5opt (26) for the 2D constant-coefficient model problem A2 of size n=1024.
(a)
MSPAI using individual probing conditions(b)
MSPAI using global probing conditions


#### 2.2.3. Global Probing Vectors

Similar to the 1D model problem, we observe that(e⊗f)TM5opt=α(e⊗f)T for certain e,f, and α. Therefore, we define global probing conditions by using (e⊗f)TM=α(e⊗f)T. Regarding the 2D case, we use global vectors resulting from the Kronecker products e¯S:=eS⊗eS, e¯N1:=eN1⊗eN1, e¯N2:=eN2⊗eN2, and e¯N3:=eN3⊗eN3. We obtain (32)S5G1:min⁡P(M)=P⁡‖e¯STM-8041e¯ST‖2,S5G2:min⁡P(M)=P⁡‖e¯N1TM-1641e¯N1T‖2,S5G3:min⁡P(M)=P⁡‖e¯N2TM-4841e¯N2T‖2,S5G4:min⁡P(M)=P⁡‖e¯N3TM-4841e¯N3T‖2.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 Figure2(b)). The global condition S5G1 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 A2, 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 ofA2, which is the block tridiagonal matrix (33)A3:=blocktridiag(-18,-18,-18|-18,1,-18|-18,-18,-18),
related to the generating function(34)a3(x,y)=1-h(x,y)4,
with h(x,y):=cos(x)+cos(y)+cos(x-y)+cos(x+y).We use the pattern ofa3(x,y) for our approximate inverse smoother (35)m3(x,y)=a+2b(cos⁡⁡(x)+cos⁡⁡(y)+cos⁡⁡(x-y)+cos⁡⁡(x+y)).
With h(x,y)∈[-1,2], this yields the minimization problem(36)min⁡a,b⁡|1-(a+2bh(x,y))(1-h(x,y)4)|.
The optimal solution is given by (a,b)=(160/153,16/153) which leads to(37)M9opt=blocktridiag(16153,16153,16153|16153,160153,16153|16153,16153,16153),
with smoothing factor 1/17=0.0588.We consider the boundary valuesh(x,y)=-1 and h(x,y)=2 to get linear conditions. A third, quadratic condition can be deduced by setting the derivative of (36) to zero and inserting its solution h′(x,y)=2-(a/4b) into (36). We obtain the following minimization conditions at high-frequency values: (38)min⁡a,b⁡{|3(a-4b)2-1|,|3(a+2b)4-1|,|a2+2b+a232b-1|}.
The direct translation into linear probing masks and their combination reveals the individual probing conditions (39)S9L1:min⁡m̂k⁡|(-34,-34,-34|-34,32,-34|-34,-34,-34)⋅m̂k-134|,S9L2:min⁡m̂k⁡|(316,316,316|316,34,316|316,316,316)m̂k-1|,S9L3:min⁡m̂k⁡|(14,14,14|14,1316,14|14,14,14)m̂k-1|,S9L4:min⁡m̂k⁡|(4164,4164,4164|4164,12,4164|4164,4164,4164)m̂k-1|,S9L5:min⁡m̂k⁡|(1,1,1|1,1,1|1,1,1)m̂k-2011|.
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 A3, by using Kronecker products of the 1D probing vectors. By observing (e⊗f)TAM9opt=α(e⊗f)T similar to Section 2.2.3, we can define conditions with (e⊗f)TAM=α(e⊗f)T for some e,f, and α.Table 2
Smoothing factors by usingM9opt (37), SPAI, and MSPAI with local probing conditions for the 2D 9-point stencil A3 of order n=1024.
M9optSPAIweightlocal probing conditionsρρS9L1ρS9L2ρS9L3ρS9L4ρS9L50.05880.16570.20.1470.1610.1500.1160.0951.00.1290.1360.0770.0660.0692.00.1300.1290.0630.0631010.1260.1260.0700.0620.0681020.1300.1270.071


## 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 formA1x=b, with A1:=tridiag(-1/2,1,-1/2). The matrix A1 is related to the generating function or symbol(10)a1(x)=1-eix+e-ix2=1-cos⁡⁡(x).
The high-frequency part of A1 is represented by (10) for x∈[π/2,π]=:I1. Hence, the optimal smoothing parameter ω in the Jacobi smoother 1-ωa1(x) is found by solving the problem(11)min⁡ω⁡max⁡x∈I1⁡|1-ωa1(x)|=min⁡ω⁡max⁡x∈I1⁡|s(x)|=ω=2/313.
The solution can be found by replacing (11) with the maximum over the two boundary values minωmax{|s(π/2)|,|s(π)|}.


### 2.1.1. Analytical Derivation of the Optimal Smoother

As approximate inverse smoother, we choose a trigonometric polynomial of the same degree(12)m1(x)=a+2bcos⁡⁡(x),
and a tridiagonal Toeplitz matrix M=tridiag(b,a,b), respectively. The smoothing condition for (12) can be written as (13)min⁡a,b⁡max⁡x∈I1⁡|1-m1(x)a1(x)|=min⁡a,b⁡max⁡x∈I1⁡|s(x)|,
or using the boundary values of I1 as(14)min⁡a,b⁡max⁡⁡{|1-m1(π2)a1(π2)|,|1-m1(π)a1(π)|,|1-m1(u)a1(u)|(π2)}=min⁡a,b⁡max⁡⁡{|a-1|,|2(a-2b)-1|,|1-m1(u)a1(u)|},
in which u∈I1 is the local extreme point of the function 1-m1(x)a1(x) with derivative equal to zero. We thus obtain for u=1/2-a/4b the quadratic condition(15)s(u)=(a+2b)28b-1,
and the overall solution (a,b)=(16/17,4/17) which leads to(16)M3opt:=tridiag(417,1617,417).
The related smoothing factor for this optimal preconditioner m(x) 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 a and b can be also seen as conditions for the entries of the matrix directly. A column k of M is described by the three values (mk-1,k,mk,k,mk+1,k)T=:m̂k, where the main diagonal entry is related to a and the upper and lower entries to b. Therefore, we can translate the two linear minimization conditions of (14) into individual probing masks ŝk on entries of M directly (17)S3L1:min⁡m̂k⁡|(0,1,0)m̂k-1|,S3L2:min⁡m̂k⁡|(-1,1,-1)m̂k-12|.
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 a=1 into the second condition yields b=1/4. We replace the denominator of the quadratic condition with this value and get the new linear condition (18)min⁡a,b⁡|(a+2b)28⋅1/4-1|=min⁡a,b⁡|(a+2b)2-2|⟶min⁡a,b⁡|a+2b-2|,
which is related to the probing condition with isotropic probing mask (19)S3L3:min⁡m̂k⁡|(1,1,1)m̂k-2|.
We call a mask ŝk to be isotropic if the minimization minm̂k|ŝkm̂k-fk| can be described by one probing vector e∈ℝn with e(𝒥k)=αŝkT, α∈ℝ and k∈{1,…,n}, in the form min𝒫(M)=𝒫∥eTM-fT∥2. Here, for instance, ŝk=(-1,1,-1) of S3L2 is related to e=(1,-1,1,-1,…)T and ŝk=(1,1,1) of S3L3 to e=(1,1,…,1)T.In our following numerical experiments, we are interested in a comparison between MSPAI with probing conditions, the optimal smootherM3opt (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 M is given by(20)min⁡m̂k⁡‖(-0.5001-0.50-0.51-0.50-0.5100-0.5)m̂k-(00100)‖2.
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 m(x)=6/5+(4/5)cos(x). 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 M3opt but still an improvement compared to Jacobi. Note that MSPAI minimizes the 2-norm of a combination of conditions while M3opt 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 matrixA1 of order n=103. 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 S3L3 weighed with ρ=102. 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 usingM3opt (16), SPAI, and MSPAI with local probing conditions for the constant coefficient system A1 of order n=103.
M3optSPAIWeightLocal probing conditionsρρS3L1∧0.7S3L2ρS3L30.05880.2501.00.1340.0832.00.1070.0771010.0950.075102


### 2.1.3. Global Probing Vectors

Considering the conditionsS3L2 and S3L3 with isotropic masks it is possible to derive conditions with low- and high-frequency global probing vectors eS:=(1,1,…,1)T and eN1:=(1,-1,1,-1,…)T for the generalized Frobenius norm minimization (5) of MSPAI (21)min⁡P(M)=P⁡‖eSTM-12eST‖2,min⁡P(M)=P⁡‖eN1TM-2eN1T‖2,
respectively. In the following, we use a different approach and derive global probing conditions with respect to the optimal derived preconditioner M3opt (16). This can also be applied to general probing vectors like eN2:=(1,0,-1,0,1,…)T and eN3:=(0,1,0,-1,0,1,…)T, representing additional high-frequency subspaces. For a given probing vector e, we observe that eTM3opt=αeT with α∈ℝ. As we want to find the smoother M, based on MSPAI with the same action on eT, we define probing conditions eTM≈αeT and obtain (22)S3G1:min⁡P(M)=P⁡‖eSTM-2417eST‖2,S3G2:min⁡P(M)=P⁡‖eN1TM-817eN1T‖2,S3G3:min⁡P(M)=P⁡‖eN2TM-1617eN2T‖2,S3G4:min⁡P(M)=P⁡‖eN3TM-1617eN3T‖2.Figure1(a) shows the impact of using these global probing conditions for matrix A1. Inducing MSPAI to satisfy S3G1 leads to a strong reduction of the smoothing factor in comparison to SPAI, close to the optimal value of M3opt. A combination of global conditions can lead to an efficient smoother as well. Again the smoothing factor stays stable for increasing values of ρ.Smoothing factorμ against subspace weight ρ. (a) considers MSPAI using global probing conditions for the constant coefficient system A1 of size n=103. (b) shows MSPAI using global probing vectors with action on B1 with varying coefficients of size n=103. MSPAI is compared to Jacobi, SPAI, and M3opt (16).
(a)
MSPAI with global probing conditions forA1(b)
MSPAI using global probing vectors with action onB1In a similar way, it is possible to derive global probing conditions with action onA, that is, min𝒫(M)=𝒫∥eTAM-αeT∥2. For probing vectors, for example, eN1,eN2, and eN3, the approximation eTA≈const·eT holds, and therefore eTAM≈const·eTM up to some boundary perturbations. Note that for an approximate inverse preconditioner which should nearly be exact on e the condition eTAM≈eT should be satisfied.As a generalization, we are interested in the impact of using global probing conditions when considering the 1D tridiagonal matrixB1 with varying coefficients and kth row defined by (23)(B1)k,::=(0,…,0,-bk-1,bk-1+bk,-bk,0,…,0),withb:=(-12-j2n)j=0,…,n.
We denote the exact conditions min𝒫(M)=𝒫∥eTB1M-(16/17)eT∥2 for the probing vectors eN1, eN2, and eN3 as S3G5, S3G6, and S3G7. The approximate conditions for α=1 are indicated by S3G̃5, S3G̃6, and S3G̃7.Figure1(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 A1.


## 2.1.1. Analytical Derivation of the Optimal Smoother

As approximate inverse smoother, we choose a trigonometric polynomial of the same degree(12)m1(x)=a+2bcos⁡⁡(x),
and a tridiagonal Toeplitz matrix M=tridiag(b,a,b), respectively. The smoothing condition for (12) can be written as (13)min⁡a,b⁡max⁡x∈I1⁡|1-m1(x)a1(x)|=min⁡a,b⁡max⁡x∈I1⁡|s(x)|,
or using the boundary values of I1 as(14)min⁡a,b⁡max⁡⁡{|1-m1(π2)a1(π2)|,|1-m1(π)a1(π)|,|1-m1(u)a1(u)|(π2)}=min⁡a,b⁡max⁡⁡{|a-1|,|2(a-2b)-1|,|1-m1(u)a1(u)|},
in which u∈I1 is the local extreme point of the function 1-m1(x)a1(x) with derivative equal to zero. We thus obtain for u=1/2-a/4b the quadratic condition(15)s(u)=(a+2b)28b-1,
and the overall solution (a,b)=(16/17,4/17) which leads to(16)M3opt:=tridiag(417,1617,417).
The related smoothing factor for this optimal preconditioner m(x) 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 a and b can be also seen as conditions for the entries of the matrix directly. A column k of M is described by the three values (mk-1,k,mk,k,mk+1,k)T=:m̂k, where the main diagonal entry is related to a and the upper and lower entries to b. Therefore, we can translate the two linear minimization conditions of (14) into individual probing masks ŝk on entries of M directly (17)S3L1:min⁡m̂k⁡|(0,1,0)m̂k-1|,S3L2:min⁡m̂k⁡|(-1,1,-1)m̂k-12|.
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 a=1 into the second condition yields b=1/4. We replace the denominator of the quadratic condition with this value and get the new linear condition (18)min⁡a,b⁡|(a+2b)28⋅1/4-1|=min⁡a,b⁡|(a+2b)2-2|⟶min⁡a,b⁡|a+2b-2|,
which is related to the probing condition with isotropic probing mask (19)S3L3:min⁡m̂k⁡|(1,1,1)m̂k-2|.
We call a mask ŝk to be isotropic if the minimization minm̂k|ŝkm̂k-fk| can be described by one probing vector e∈ℝn with e(𝒥k)=αŝkT, α∈ℝ and k∈{1,…,n}, in the form min𝒫(M)=𝒫∥eTM-fT∥2. Here, for instance, ŝk=(-1,1,-1) of S3L2 is related to e=(1,-1,1,-1,…)T and ŝk=(1,1,1) of S3L3 to e=(1,1,…,1)T.In our following numerical experiments, we are interested in a comparison between MSPAI with probing conditions, the optimal smootherM3opt (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 M is given by(20)min⁡m̂k⁡‖(-0.5001-0.50-0.51-0.50-0.5100-0.5)m̂k-(00100)‖2.
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 m(x)=6/5+(4/5)cos(x). 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 M3opt but still an improvement compared to Jacobi. Note that MSPAI minimizes the 2-norm of a combination of conditions while M3opt 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 matrixA1 of order n=103. 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 S3L3 weighed with ρ=102. 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 usingM3opt (16), SPAI, and MSPAI with local probing conditions for the constant coefficient system A1 of order n=103.
M3optSPAIWeightLocal probing conditionsρρS3L1∧0.7S3L2ρS3L30.05880.2501.00.1340.0832.00.1070.0771010.0950.075102


## 2.1.3. Global Probing Vectors

Considering the conditionsS3L2 and S3L3 with isotropic masks it is possible to derive conditions with low- and high-frequency global probing vectors eS:=(1,1,…,1)T and eN1:=(1,-1,1,-1,…)T for the generalized Frobenius norm minimization (5) of MSPAI (21)min⁡P(M)=P⁡‖eSTM-12eST‖2,min⁡P(M)=P⁡‖eN1TM-2eN1T‖2,
respectively. In the following, we use a different approach and derive global probing conditions with respect to the optimal derived preconditioner M3opt (16). This can also be applied to general probing vectors like eN2:=(1,0,-1,0,1,…)T and eN3:=(0,1,0,-1,0,1,…)T, representing additional high-frequency subspaces. For a given probing vector e, we observe that eTM3opt=αeT with α∈ℝ. As we want to find the smoother M, based on MSPAI with the same action on eT, we define probing conditions eTM≈αeT and obtain (22)S3G1:min⁡P(M)=P⁡‖eSTM-2417eST‖2,S3G2:min⁡P(M)=P⁡‖eN1TM-817eN1T‖2,S3G3:min⁡P(M)=P⁡‖eN2TM-1617eN2T‖2,S3G4:min⁡P(M)=P⁡‖eN3TM-1617eN3T‖2.Figure1(a) shows the impact of using these global probing conditions for matrix A1. Inducing MSPAI to satisfy S3G1 leads to a strong reduction of the smoothing factor in comparison to SPAI, close to the optimal value of M3opt. A combination of global conditions can lead to an efficient smoother as well. Again the smoothing factor stays stable for increasing values of ρ.Smoothing factorμ against subspace weight ρ. (a) considers MSPAI using global probing conditions for the constant coefficient system A1 of size n=103. (b) shows MSPAI using global probing vectors with action on B1 with varying coefficients of size n=103. MSPAI is compared to Jacobi, SPAI, and M3opt (16).
(a)
MSPAI with global probing conditions forA1(b)
MSPAI using global probing vectors with action onB1In a similar way, it is possible to derive global probing conditions with action onA, that is, min𝒫(M)=𝒫∥eTAM-αeT∥2. For probing vectors, for example, eN1,eN2, and eN3, the approximation eTA≈const·eT holds, and therefore eTAM≈const·eTM up to some boundary perturbations. Note that for an approximate inverse preconditioner which should nearly be exact on e the condition eTAM≈eT should be satisfied.As a generalization, we are interested in the impact of using global probing conditions when considering the 1D tridiagonal matrixB1 with varying coefficients and kth row defined by (23)(B1)k,::=(0,…,0,-bk-1,bk-1+bk,-bk,0,…,0),withb:=(-12-j2n)j=0,…,n.
We denote the exact conditions min𝒫(M)=𝒫∥eTB1M-(16/17)eT∥2 for the probing vectors eN1, eN2, and eN3 as S3G5, S3G6, and S3G7. The approximate conditions for α=1 are indicated by S3G̃5, S3G̃6, and S3G̃7.Figure1(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 A1.


## 2.2. A 2D Model Problem

We consider the block tridiagonal matrix(24)A2:=A1⊗I+I⊗A12=blocktridiag(0,-14,0|-14,1,-14|0,-14,0),
with A2∈ℝn2×n2 related to the generating function (25)a2(x,y)=1-cos⁡⁡(x)+cos⁡⁡(y)2.
As before the smoothing corresponds to the rectangle I2:={(x,y)∣x∈[π/2,π]∧y∈[0,π]}. Hence, the solution of minωmaxx,y∈I2|1-ωa2(x,y)| 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 polynomialm2(x,y)=a+2b(cos(x)+cos(y)). Via the minimization 1-m2(x,y)a2(x,y) over I2, the solution for the optimal smoothing preconditioner is given by(26)M5opt:=blocktridiag(0,841,0|841,4841,841|0,841,0),
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(27)min⁡a,b⁡|1-(a+2bh(x,y))(1-h(x,y)2)|,withcos⁡⁡(x)+cos⁡⁡(y)=:h(x,y)∈[-2,1].
The corners of I2 result in the minimization conditions(28)min⁡a,b⁡{|2a-8b-1|,|a2+b-1|,|a2+b+a216b-1|}.
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 a and b. We obtain a=6b and the quadratic term can be replaced either by 6ba/b=6a or (6b)2/b=36b. Consequently, we end up with two additional linear conditions(29)min⁡a,b⁡{|7a8+b-1|,|a2+13b4-1|}.Once again, we can see the minimizations (28) and (29) on a and b as conditions for the entries of the smoother directly. A column k of M is now described by the five degrees of freedom (0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T=:m̂k. Thus, we gain the following local probing conditions for the 5-point stencil A2: (30)S5L1:min⁡m̂k⁡|(0,-2,0|-2,2,-2|0,-2,0)m̂k-1|,S5L2:min⁡m̂k⁡|(0,14,0|14,12,14|0,14,0)m̂k-1|,S5L3:min⁡m̂k⁡|(0,14,0|14,78,14|0,14,0)m̂k-1|,S5L4:min⁡m̂k⁡|(0,1316,0|1316,12,1316|0,1316,0)m̂k-1|.
It is possible to combine the conditions S5L2 and S5L4 to obtain an isotropic mask to which the low-frequency probing vector eS corresponds to. To derive this, we add the diagonal and subdiagonal values of the weighed conditions r·S5L2 and s·S5L4 and set them to be equal. The isotropy condition leads to the equation (r+s)/2=r/4+13s/16 with the solution r=5s/4. We obtain a=b=(r+s)/2=9s/8 and r+s=5s/4+s=9s/4 for the right hand side. This leads to the additional individual isotropic probing condition (31)S5L5:min⁡m̂k⁡|(0,1,0|1,1,1|0,1,0)m̂k-2|.Let us considerA2 of size n=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 M5opt.Smoothing factorμ against subspace weight ρ. MSPAI using individual (a) and global (b) probing conditions compared to Jacobi, SPAI, and M5opt (26) for the 2D constant-coefficient model problem A2 of size n=1024.
(a)
MSPAI using individual probing conditions(b)
MSPAI using global probing conditions


### 2.2.3. Global Probing Vectors

Similar to the 1D model problem, we observe that(e⊗f)TM5opt=α(e⊗f)T for certain e,f, and α. Therefore, we define global probing conditions by using (e⊗f)TM=α(e⊗f)T. Regarding the 2D case, we use global vectors resulting from the Kronecker products e¯S:=eS⊗eS, e¯N1:=eN1⊗eN1, e¯N2:=eN2⊗eN2, and e¯N3:=eN3⊗eN3. We obtain (32)S5G1:min⁡P(M)=P⁡‖e¯STM-8041e¯ST‖2,S5G2:min⁡P(M)=P⁡‖e¯N1TM-1641e¯N1T‖2,S5G3:min⁡P(M)=P⁡‖e¯N2TM-4841e¯N2T‖2,S5G4:min⁡P(M)=P⁡‖e¯N3TM-4841e¯N3T‖2.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 Figure2(b)). The global condition S5G1 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 A2, 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 ofA2, which is the block tridiagonal matrix (33)A3:=blocktridiag(-18,-18,-18|-18,1,-18|-18,-18,-18),
related to the generating function(34)a3(x,y)=1-h(x,y)4,
with h(x,y):=cos(x)+cos(y)+cos(x-y)+cos(x+y).We use the pattern ofa3(x,y) for our approximate inverse smoother (35)m3(x,y)=a+2b(cos⁡⁡(x)+cos⁡⁡(y)+cos⁡⁡(x-y)+cos⁡⁡(x+y)).
With h(x,y)∈[-1,2], this yields the minimization problem(36)min⁡a,b⁡|1-(a+2bh(x,y))(1-h(x,y)4)|.
The optimal solution is given by (a,b)=(160/153,16/153) which leads to(37)M9opt=blocktridiag(16153,16153,16153|16153,160153,16153|16153,16153,16153),
with smoothing factor 1/17=0.0588.We consider the boundary valuesh(x,y)=-1 and h(x,y)=2 to get linear conditions. A third, quadratic condition can be deduced by setting the derivative of (36) to zero and inserting its solution h′(x,y)=2-(a/4b) into (36). We obtain the following minimization conditions at high-frequency values: (38)min⁡a,b⁡{|3(a-4b)2-1|,|3(a+2b)4-1|,|a2+2b+a232b-1|}.
The direct translation into linear probing masks and their combination reveals the individual probing conditions (39)S9L1:min⁡m̂k⁡|(-34,-34,-34|-34,32,-34|-34,-34,-34)⋅m̂k-134|,S9L2:min⁡m̂k⁡|(316,316,316|316,34,316|316,316,316)m̂k-1|,S9L3:min⁡m̂k⁡|(14,14,14|14,1316,14|14,14,14)m̂k-1|,S9L4:min⁡m̂k⁡|(4164,4164,4164|4164,12,4164|4164,4164,4164)m̂k-1|,S9L5:min⁡m̂k⁡|(1,1,1|1,1,1|1,1,1)m̂k-2011|.
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 A3, by using Kronecker products of the 1D probing vectors. By observing (e⊗f)TAM9opt=α(e⊗f)T similar to Section 2.2.3, we can define conditions with (e⊗f)TAM=α(e⊗f)T for some e,f, and α.Table 2
Smoothing factors by usingM9opt (37), SPAI, and MSPAI with local probing conditions for the 2D 9-point stencil A3 of order n=1024.
M9optSPAIweightlocal probing conditionsρρS9L1ρS9L2ρS9L3ρS9L4ρS9L50.05880.16570.20.1470.1610.1500.1160.0951.00.1290.1360.0770.0660.0692.00.1300.1290.0630.0631010.1260.1260.0700.0620.0681020.1300.1270.071


## 2.2.1. Analytical Derivation of the Optimal Smoother

As approximate inverse smoother, we use the trigonometric polynomialm2(x,y)=a+2b(cos(x)+cos(y)). Via the minimization 1-m2(x,y)a2(x,y) over I2, the solution for the optimal smoothing preconditioner is given by(26)M5opt:=blocktridiag(0,841,0|841,4841,841|0,841,0),
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(27)min⁡a,b⁡|1-(a+2bh(x,y))(1-h(x,y)2)|,withcos⁡⁡(x)+cos⁡⁡(y)=:h(x,y)∈[-2,1].
The corners of I2 result in the minimization conditions(28)min⁡a,b⁡{|2a-8b-1|,|a2+b-1|,|a2+b+a216b-1|}.
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 a and b. We obtain a=6b and the quadratic term can be replaced either by 6ba/b=6a or (6b)2/b=36b. Consequently, we end up with two additional linear conditions(29)min⁡a,b⁡{|7a8+b-1|,|a2+13b4-1|}.Once again, we can see the minimizations (28) and (29) on a and b as conditions for the entries of the smoother directly. A column k of M is now described by the five degrees of freedom (0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T=:m̂k. Thus, we gain the following local probing conditions for the 5-point stencil A2: (30)S5L1:min⁡m̂k⁡|(0,-2,0|-2,2,-2|0,-2,0)m̂k-1|,S5L2:min⁡m̂k⁡|(0,14,0|14,12,14|0,14,0)m̂k-1|,S5L3:min⁡m̂k⁡|(0,14,0|14,78,14|0,14,0)m̂k-1|,S5L4:min⁡m̂k⁡|(0,1316,0|1316,12,1316|0,1316,0)m̂k-1|.
It is possible to combine the conditions S5L2 and S5L4 to obtain an isotropic mask to which the low-frequency probing vector eS corresponds to. To derive this, we add the diagonal and subdiagonal values of the weighed conditions r·S5L2 and s·S5L4 and set them to be equal. The isotropy condition leads to the equation (r+s)/2=r/4+13s/16 with the solution r=5s/4. We obtain a=b=(r+s)/2=9s/8 and r+s=5s/4+s=9s/4 for the right hand side. This leads to the additional individual isotropic probing condition (31)S5L5:min⁡m̂k⁡|(0,1,0|1,1,1|0,1,0)m̂k-2|.Let us considerA2 of size n=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 M5opt.Smoothing factorμ against subspace weight ρ. MSPAI using individual (a) and global (b) probing conditions compared to Jacobi, SPAI, and M5opt (26) for the 2D constant-coefficient model problem A2 of size n=1024.
(a)
MSPAI using individual probing conditions(b)
MSPAI using global probing conditions


## 2.2.3. Global Probing Vectors

Similar to the 1D model problem, we observe that(e⊗f)TM5opt=α(e⊗f)T for certain e,f, and α. Therefore, we define global probing conditions by using (e⊗f)TM=α(e⊗f)T. Regarding the 2D case, we use global vectors resulting from the Kronecker products e¯S:=eS⊗eS, e¯N1:=eN1⊗eN1, e¯N2:=eN2⊗eN2, and e¯N3:=eN3⊗eN3. We obtain (32)S5G1:min⁡P(M)=P⁡‖e¯STM-8041e¯ST‖2,S5G2:min⁡P(M)=P⁡‖e¯N1TM-1641e¯N1T‖2,S5G3:min⁡P(M)=P⁡‖e¯N2TM-4841e¯N2T‖2,S5G4:min⁡P(M)=P⁡‖e¯N3TM-4841e¯N3T‖2.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 Figure2(b)). The global condition S5G1 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 A2, 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 ofA2, which is the block tridiagonal matrix (33)A3:=blocktridiag(-18,-18,-18|-18,1,-18|-18,-18,-18),
related to the generating function(34)a3(x,y)=1-h(x,y)4,
with h(x,y):=cos(x)+cos(y)+cos(x-y)+cos(x+y).We use the pattern ofa3(x,y) for our approximate inverse smoother (35)m3(x,y)=a+2b(cos⁡⁡(x)+cos⁡⁡(y)+cos⁡⁡(x-y)+cos⁡⁡(x+y)).
With h(x,y)∈[-1,2], this yields the minimization problem(36)min⁡a,b⁡|1-(a+2bh(x,y))(1-h(x,y)4)|.
The optimal solution is given by (a,b)=(160/153,16/153) which leads to(37)M9opt=blocktridiag(16153,16153,16153|16153,160153,16153|16153,16153,16153),
with smoothing factor 1/17=0.0588.We consider the boundary valuesh(x,y)=-1 and h(x,y)=2 to get linear conditions. A third, quadratic condition can be deduced by setting the derivative of (36) to zero and inserting its solution h′(x,y)=2-(a/4b) into (36). We obtain the following minimization conditions at high-frequency values: (38)min⁡a,b⁡{|3(a-4b)2-1|,|3(a+2b)4-1|,|a2+2b+a232b-1|}.
The direct translation into linear probing masks and their combination reveals the individual probing conditions (39)S9L1:min⁡m̂k⁡|(-34,-34,-34|-34,32,-34|-34,-34,-34)⋅m̂k-134|,S9L2:min⁡m̂k⁡|(316,316,316|316,34,316|316,316,316)m̂k-1|,S9L3:min⁡m̂k⁡|(14,14,14|14,1316,14|14,14,14)m̂k-1|,S9L4:min⁡m̂k⁡|(4164,4164,4164|4164,12,4164|4164,4164,4164)m̂k-1|,S9L5:min⁡m̂k⁡|(1,1,1|1,1,1|1,1,1)m̂k-2011|.
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 A3, by using Kronecker products of the 1D probing vectors. By observing (e⊗f)TAM9opt=α(e⊗f)T similar to Section 2.2.3, we can define conditions with (e⊗f)TAM=α(e⊗f)T for some e,f, and α.Table 2
Smoothing factors by usingM9opt (37), SPAI, and MSPAI with local probing conditions for the 2D 9-point stencil A3 of order n=1024.
M9optSPAIweightlocal probing conditionsρρS9L1ρS9L2ρS9L3ρS9L4ρS9L50.05880.16570.20.1470.1610.1500.1160.0951.00.1290.1360.0770.0660.0692.00.1300.1290.0630.0631010.1260.1260.0700.0620.0681020.1300.1270.071


## 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 generalH to suppress a reconstruction on the noise space. Following [20], such a preconditioner M should have the following properties:  (i)
M≈|H|-1 on the signal subspace with |H|=HTH, and(ii)
M≈I or M≈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 vectoreS:=(1,1,…,1)T 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 useeN1:=(1,-1,1,-1,…)T, eN2:=(1,0,-1,0,1,…)T, and eN3:=(0,1,0,-1,0,…)T as typical vectors related to fast oscillations. In case of using probing masks, ŝk=(-1,1,-1) and ŝk=(1,0,-1) are related to eN1 and eN2,  eN3, respectively.For higher dimensional problems, probing vectors typically result from a Kronecker product of 1D probing vectors. The conditions in MSPAI are given byAM≈I in order to derive a good preconditioner and fast convergence on the signal subspace and ρeNTM≈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 matrixH1:=tridiag(1/2,1,1/2) related to the symbol (40)h1(x)=1+eix+e-ix2=1+cos⁡⁡(x).
As approximate inverse preconditioner, we choose the trigonometric polynomial of symbol m4(x)=a+2bcos(x) related to a Toeplitz matrix M=tridiag(b,a,b). The entries a and b should be determined such that the preconditioner M acts both nearly as the inverse on the signal subspace and as zero on the noise subspace. For x∈{0,π/2,π}, this leads to the minimization conditions(41)min⁡a,b⁡{|m4(0)h1(0)-1|,|m4(π2)h1(π2)-1|,ρ|m4(π)|}=min⁡a,b⁡{|2(a+2b)-1|,|a-1|,ρ|a-2b|}.
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 (42)(2a+4b-1)=(a-1)=-ρ(a-2b).
Thus, we are able to derive the optimal regularizing preconditioner(43)M3ρ:=tridiag(ρ-0.52+5ρ,2+2ρ2+5ρ,ρ-0.52+5ρ).
In the limiting case, for only preconditioning on the signal space (ρ=0), we get M30=tridiag(-1/4,1,-1/4), and for only regularizing on the noise subspace (ρ=∞) we get M3∞=tridiag(1/5,2/5,1/5)=(2/5)H1=(2/5)H1T. 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 k is described by the three degrees of freedom m̂k:=(mk-1,k,mk,k,mk+1,k)T. 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 ŝk yields (44)R3L1:min⁡m̂k⁡|(-1,1,-1)m̂k|.
It is possible to transfer this individual regularization condition to the global condition eTM≈0 with probing vector eN1=(1,-1,1,-1,…)T, representing high-frequency noisy components. Thus, we force M to act as approximately zero on the noise subspace. We refer to the high-frequency global probing condition as
(45)R3G1:min⁡P(M)=P⁡‖eN1TM‖2.As a first numerical example, we examine the given 1D blur operatorH12k with k>0 for a given vector x representing the original data. We use the MSPAI regularization property to reconstruct two different signals perturbed by random noise. By using powers of H, 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 MMT after the computation of M; that is, the preconditioner is (MMT)k. 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 x and the reconstruction x̃ 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 operatorH18 and original data x1 of size n=103 affected by random noise of order 0.01%. Here, x1 denotes the smooth signal
(46)x1:=(sin⁡⁡(4πjn))j=1,…,n.
The measurement of the error ∥x1-x̃1∥2 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 M3ρ 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∥x1-x̃1∥2 for the 1D operator H18 and original data x1 of size n=103 affected by random noise of order 0.01%.
Regularization methodOptimal valueReached at iterationCG(H18)2.1870·10-45CG(H18TH18)1.8998·10-462PCG(H18) with M3ρρ=0.42.3699·10-44ρ=11.9746·10-412ρ=51.9236·10-434ρ=1021.9226·10-448Figure 3
Reconstruction error against the iteration for the 1D operatorH18 and original data x1 of size n=103 affected by random noise of order 0.01%. CG is compared to PCG using the optimal preconditioner M3ρ (43) for different weights ρ.In our second 1D example, we observe the behavior for the blur operatorH14 and a signal which has strongly increasing entries near the boundaries and an additional nonanalytic point in the middle
(47)x2:=(1j-1n+1-j+1n/2+1/4-j+sin⁡⁡(4jπn))j=1,…,n.
The problem has size n=103 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 n/2, but deteriorates near the boundary. Therefore, we introduce some correction at the boundary by changing the components (eN1)1=τ and (eN1)n=±τ, 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 (eN1)n/2=0 and (eN1)n/2±1=1/2(eN1)n/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 x̃2 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∥x2-x̃2∥2 for the 1D blur operator H14 and the original data x2. The problem has size n=103 and is affected by random noise of order 0.1%. ρ=102.
Regularization methodOptimal valueReached at iter.Tikhonov(γ=0.3)0.1286—CG(H14)0.12408PCG(H14)ρR3G1no boundary correction0.1210130τ=0.250.043358τ=0.25, middle correction0.026352Figure 4
Reconstruction error against the iteration for the 1D blur operatorH14 and the original data x2 of size n=103 affected by random noise of order 0.1%. CG is compared to CG for the normal equations and to PCG using MSPAI satisfying R3G1 with ρ=102.Likewise, it is possible to use subspace corrections within mask probing. We change the individual probing conditionR3L1 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 M3ρ 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 x 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 ofx1 based on the blur operator B2 of size n=103 affected by random noise of order 1% and 0.1%, respectively. B2 denotes the tridiagonal system with varying coefficients whose kth row is given by (48)(B2)k,::=(0,…,0,wk-1,wk-1+wk,wk,0,…,0),forw:=(2j2n2+1)j=0,…,n,
with k∈{1,…,n}. 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 R3G1 to CG(B2), we are able to reconstruct x1 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.Reconstruction error against the iteration for the 1D blur operatorB2 and the original data x1 of size n=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 B2 and B2TB2 using MSPAI with R3G1.
(a)
Noise level 0.1%(b)
Noise level 1%


### 3.2. A 2D Model Problem

In 2D, we consider the matrix(49)H2:=H1⊗I+I⊗H12=blocktridiag(0,14,0|14,1,14|0,14,0),
with its corresponding symbol h2(x,y)=1+cos(x)/2+cos(y)/2 and preconditioner m5(x,y)=a+2b(cos(x)+cos(y)). With cos(x)+cos(y)=:h(x,y)∈[-1,2]=:I3, we have to minimize(50)min⁡a,b⁡max⁡x,y∈I3⁡|1-(a+2bh(x,y))(1-h(x,y)2)|,
relative to the signal subspace and mina,bρ|a-4b| for h(x,y)=-2 relative to the noise subspace. This yields the optimal block tridiagonal preconditioner(51)M5ρ=blocktridiag(0,2ρ-32,0|2ρ-32,9+8ρ,2ρ-32|0,2ρ-32,0).The individual probing masks for a reduced inner column of the preconditionerm̂k:=(0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T 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 (52)R5L1:min⁡m̂k⁡|((-1,1,-1)⊗(-1,1,-1))m̂k|,R5G1:min⁡P(M)=P⁡‖(eN1⊗eN1)TM‖2,
respectively.We observe the behavior for the 2D problem with blur operatorH24 and original data x1 of size n=502 affected by random noise of order 0.1%. Again, for M5ρ, with increasing values of ρ, we approximate the convergence of CG using normal equations (refer to Figure 6). Applying MSPAI with R5G1 to PCG it is almost possible to reach the value of CG(H24TH24) but in less iterations. Higher values of ρ lead to smooth and broad convergence curves similar to M5ρ. CG has its optimal value after 9 iterations with error 2.219·10-3, MSPAI using R5G1, ρ=1 after 12 iterations with value 1.496·10-3, and CG using normal equations reaches the error 1.282·10-3 after 44 iterations.Reconstruction error against the iteration for the 2D problem with blur operatorH24 and original data x1 of size n=502 affected by random noise of order 0.1%. CG is compared to CG for the normal equations and to M5ρ (a) and MSPAI satisfying ρR5G1 (b), respectively.
(a)
PCG using the optimal preconditionerM5ρ(51)(b)
PCG using MSPAI satisfyingρR5G1


### 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 thederiv2 example which is a discretization of a first kind Fredholm integral equation. We choose deriv2(n,case) with case=2. Our problem has size n=2·103 and is affected with random noise of order 0.001% and 0.01%, respectively. Note that the system matrix K is dense. We apply both CG and PCG to KTK and force the MSPAI to act as approximately zero on the noise subspace by using the high-frequency probing conditions R3L1 and R3L2:minm̂k|(1,0,-1)m̂k| simultaneously. We weigh the subspace with ρ=103and apply the symmetric preconditioner M+MT to the normal equations. To avoid deterioration at the boundary, we adjust M by resetting the values M1,1=M2,2, M2,1=M3,2, Mn,n=Mn-1,n-1, and Mn-1,n=Mn-2,n-1.Following Figure7, we are able to achieve better reconstruction g̃ of the original data g and in fewer iterations when applying MSPAI. For other noise levels and for deriv(n,1), we observed similar behavior.Reconstruction error against the iteration for thederiv2 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 R3L1, R3L2, and ρ=103 for the normal equations. MSPAI is symmetrized via M+MT.
(a)
Noise level 0.001%(b)
Noise level 0.01%Let us focus on theblur [24] test problem as a second example which is deblurring images degraded by atmospheric turbulence blur. The matrix G is an n2-by-n2 symmetric, doubly block Toeplitz matrix that models blurring of an n-by-n image by a Gaussian point spread function. The parameter σ controls the width of G and thus the amount of smoothing and ill-posedness. G=A1⊗A1 is symmetric block banded and possibly positive definite depending on n and σ. We choose the problem to be of size n=150 with bandwidth 4 and σ=1, that is, we invoke blur(150,4,1) and G is of size 1502×1502. The original data vector is denoted by p.We compare the unpreconditioned CG to PCG both forG and for the normal equations GTG. In case of preconditioning G, MSPAI is symmetrized via M+MT and for GTG the preconditioner MTM is applied. For MSPAI, we impose the blocktridiagonal pattern m̂k:=(0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T. In view of the structure of G, we build the high-frequency subspace by Kronecker products of oscillatory probing vectors in the regularizing global conditions R5G1 and the new ones R5G2:min𝒫(M)=𝒫∥(eN1⊗eN2)TM∥2 and R5G3:min𝒫(M)=𝒫∥(eN2⊗eN1)TM∥2, all weighed with ρ=1.Following Figure8 and Table 5, we obtain better reconstruction p̃ when applying MSPAI in contrast to the unpreconditioned CG for G or GTG. We observed similar behavior for other values of band and σ.Table 5
Optimal reconstruction error∥p-p̃∥2 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.
Regularization methodNoise level 0.01%Noise level 0.1%∥p-p̃∥2at it.∥p-p̃∥2at it.CG(G)25.637341.8281CG(GTG)24.1231937.8072PCG(G), MSPAI with R5G1,G2,G323.620639.0071PCG(GTG), MSPAI with R5G1,G2,G323.7785138.0233Regularization methodNoise level 1%Noise level 10%∥p-p̃∥2at it.∥p-p̃∥2at it.CG(G)591.0113.820·1031CG(GTG)82.511.253·1031PCG(G), MSPAI with R5G1,G2,G3189.814.104·1031PCG(GTG), MSPAI with R5G1,G2,G369.710.962·1031Reconstruction error against the iteration for theblur 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 R5G1, R5G2, R5G3, and PCG using MSPAI for normal equations with ρ=1. We use M+MT as preconditioner within PCG(G) and MTM for PCG(GTG).
(a)
Noise level 0.01%(b)
Noise level 1%As a last example we reduce the 2Dblur example of Hansen to the 1D case. For the blur operator we take the 1D analogon of G and reduce the 2D right hand side p to appropriate size. The consideration of example H14 with original data x2 (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 M=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 jth row (0,…,0,rj-1,1,rj,0,…,0). Here, rj≈1/2 near continuous components xj, but rj≈0 near discontinuities. In case that we know the original data x we define Mr with jth row (53)(Mr)j,::=(0,…,0,rj-1,1,rj,0,…,0),forr:=(12⋅11+ρ|xj-xj+1|)j=1,…,n-1.
Otherwise we define the preconditioner Mr̃ with jth row (54)(Mr̃)j,::=(0,…,0,r̃j-1,1,r̃j,0,…,0),viar̃:=(12⋅11+(ρ|bj-bj+1|)k)j=1,…,n-1.
The parameters ρ and k have to be chosen in such a way that discontinuities are revealed as good as possible.Following Figure9, the reconstruction is strongly improved if the preconditioner is adjusted relative to the discontinuities of x. Also, using only the observed data improves the reconstruction. Therefore, we could also consider an iterative process, where a first approximation x1 on x is used to define the tridiagonal preconditioner M1 delivering an improved approximation x2 which again gives a new preconditioner M2, and so on.Reconstruction error against iterations for the 1Dblur problem of size n=103 and band=4 affected by random noise of order 0.1%. G has σ=1 (a) and σ=10 (b), respectively. CG is compared to PCG using Mr and Mr̃, respectively, for G and GTG.
(a)
G with σ=1(b)
G with σ=10


## 3.1. A 1D Model Problem

Let us consider the matrixH1:=tridiag(1/2,1,1/2) related to the symbol (40)h1(x)=1+eix+e-ix2=1+cos⁡⁡(x).
As approximate inverse preconditioner, we choose the trigonometric polynomial of symbol m4(x)=a+2bcos(x) related to a Toeplitz matrix M=tridiag(b,a,b). The entries a and b should be determined such that the preconditioner M acts both nearly as the inverse on the signal subspace and as zero on the noise subspace. For x∈{0,π/2,π}, this leads to the minimization conditions(41)min⁡a,b⁡{|m4(0)h1(0)-1|,|m4(π2)h1(π2)-1|,ρ|m4(π)|}=min⁡a,b⁡{|2(a+2b)-1|,|a-1|,ρ|a-2b|}.
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 (42)(2a+4b-1)=(a-1)=-ρ(a-2b).
Thus, we are able to derive the optimal regularizing preconditioner(43)M3ρ:=tridiag(ρ-0.52+5ρ,2+2ρ2+5ρ,ρ-0.52+5ρ).
In the limiting case, for only preconditioning on the signal space (ρ=0), we get M30=tridiag(-1/4,1,-1/4), and for only regularizing on the noise subspace (ρ=∞) we get M3∞=tridiag(1/5,2/5,1/5)=(2/5)H1=(2/5)H1T. 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 k is described by the three degrees of freedom m̂k:=(mk-1,k,mk,k,mk+1,k)T. 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 ŝk yields (44)R3L1:min⁡m̂k⁡|(-1,1,-1)m̂k|.
It is possible to transfer this individual regularization condition to the global condition eTM≈0 with probing vector eN1=(1,-1,1,-1,…)T, representing high-frequency noisy components. Thus, we force M to act as approximately zero on the noise subspace. We refer to the high-frequency global probing condition as
(45)R3G1:min⁡P(M)=P⁡‖eN1TM‖2.As a first numerical example, we examine the given 1D blur operatorH12k with k>0 for a given vector x representing the original data. We use the MSPAI regularization property to reconstruct two different signals perturbed by random noise. By using powers of H, 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 MMT after the computation of M; that is, the preconditioner is (MMT)k. 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 x and the reconstruction x̃ 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 operatorH18 and original data x1 of size n=103 affected by random noise of order 0.01%. Here, x1 denotes the smooth signal
(46)x1:=(sin⁡⁡(4πjn))j=1,…,n.
The measurement of the error ∥x1-x̃1∥2 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 M3ρ 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∥x1-x̃1∥2 for the 1D operator H18 and original data x1 of size n=103 affected by random noise of order 0.01%.
Regularization methodOptimal valueReached at iterationCG(H18)2.1870·10-45CG(H18TH18)1.8998·10-462PCG(H18) with M3ρρ=0.42.3699·10-44ρ=11.9746·10-412ρ=51.9236·10-434ρ=1021.9226·10-448Figure 3
Reconstruction error against the iteration for the 1D operatorH18 and original data x1 of size n=103 affected by random noise of order 0.01%. CG is compared to PCG using the optimal preconditioner M3ρ (43) for different weights ρ.In our second 1D example, we observe the behavior for the blur operatorH14 and a signal which has strongly increasing entries near the boundaries and an additional nonanalytic point in the middle
(47)x2:=(1j-1n+1-j+1n/2+1/4-j+sin⁡⁡(4jπn))j=1,…,n.
The problem has size n=103 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 n/2, but deteriorates near the boundary. Therefore, we introduce some correction at the boundary by changing the components (eN1)1=τ and (eN1)n=±τ, 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 (eN1)n/2=0 and (eN1)n/2±1=1/2(eN1)n/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 x̃2 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∥x2-x̃2∥2 for the 1D blur operator H14 and the original data x2. The problem has size n=103 and is affected by random noise of order 0.1%. ρ=102.
Regularization methodOptimal valueReached at iter.Tikhonov(γ=0.3)0.1286—CG(H14)0.12408PCG(H14)ρR3G1no boundary correction0.1210130τ=0.250.043358τ=0.25, middle correction0.026352Figure 4
Reconstruction error against the iteration for the 1D blur operatorH14 and the original data x2 of size n=103 affected by random noise of order 0.1%. CG is compared to CG for the normal equations and to PCG using MSPAI satisfying R3G1 with ρ=102.Likewise, it is possible to use subspace corrections within mask probing. We change the individual probing conditionR3L1 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 M3ρ 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 x 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 ofx1 based on the blur operator B2 of size n=103 affected by random noise of order 1% and 0.1%, respectively. B2 denotes the tridiagonal system with varying coefficients whose kth row is given by (48)(B2)k,::=(0,…,0,wk-1,wk-1+wk,wk,0,…,0),forw:=(2j2n2+1)j=0,…,n,
with k∈{1,…,n}. 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 R3G1 to CG(B2), we are able to reconstruct x1 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.Reconstruction error against the iteration for the 1D blur operatorB2 and the original data x1 of size n=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 B2 and B2TB2 using MSPAI with R3G1.
(a)
Noise level 0.1%(b)
Noise level 1%


## 3.2. A 2D Model Problem

In 2D, we consider the matrix(49)H2:=H1⊗I+I⊗H12=blocktridiag(0,14,0|14,1,14|0,14,0),
with its corresponding symbol h2(x,y)=1+cos(x)/2+cos(y)/2 and preconditioner m5(x,y)=a+2b(cos(x)+cos(y)). With cos(x)+cos(y)=:h(x,y)∈[-1,2]=:I3, we have to minimize(50)min⁡a,b⁡max⁡x,y∈I3⁡|1-(a+2bh(x,y))(1-h(x,y)2)|,
relative to the signal subspace and mina,bρ|a-4b| for h(x,y)=-2 relative to the noise subspace. This yields the optimal block tridiagonal preconditioner(51)M5ρ=blocktridiag(0,2ρ-32,0|2ρ-32,9+8ρ,2ρ-32|0,2ρ-32,0).The individual probing masks for a reduced inner column of the preconditionerm̂k:=(0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T 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 (52)R5L1:min⁡m̂k⁡|((-1,1,-1)⊗(-1,1,-1))m̂k|,R5G1:min⁡P(M)=P⁡‖(eN1⊗eN1)TM‖2,
respectively.We observe the behavior for the 2D problem with blur operatorH24 and original data x1 of size n=502 affected by random noise of order 0.1%. Again, for M5ρ, with increasing values of ρ, we approximate the convergence of CG using normal equations (refer to Figure 6). Applying MSPAI with R5G1 to PCG it is almost possible to reach the value of CG(H24TH24) but in less iterations. Higher values of ρ lead to smooth and broad convergence curves similar to M5ρ. CG has its optimal value after 9 iterations with error 2.219·10-3, MSPAI using R5G1, ρ=1 after 12 iterations with value 1.496·10-3, and CG using normal equations reaches the error 1.282·10-3 after 44 iterations.Reconstruction error against the iteration for the 2D problem with blur operatorH24 and original data x1 of size n=502 affected by random noise of order 0.1%. CG is compared to CG for the normal equations and to M5ρ (a) and MSPAI satisfying ρR5G1 (b), respectively.
(a)
PCG using the optimal preconditionerM5ρ(51)(b)
PCG using MSPAI satisfyingρR5G1


## 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 thederiv2 example which is a discretization of a first kind Fredholm integral equation. We choose deriv2(n,case) with case=2. Our problem has size n=2·103 and is affected with random noise of order 0.001% and 0.01%, respectively. Note that the system matrix K is dense. We apply both CG and PCG to KTK and force the MSPAI to act as approximately zero on the noise subspace by using the high-frequency probing conditions R3L1 and R3L2:minm̂k|(1,0,-1)m̂k| simultaneously. We weigh the subspace with ρ=103and apply the symmetric preconditioner M+MT to the normal equations. To avoid deterioration at the boundary, we adjust M by resetting the values M1,1=M2,2, M2,1=M3,2, Mn,n=Mn-1,n-1, and Mn-1,n=Mn-2,n-1.Following Figure7, we are able to achieve better reconstruction g̃ of the original data g and in fewer iterations when applying MSPAI. For other noise levels and for deriv(n,1), we observed similar behavior.Reconstruction error against the iteration for thederiv2 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 R3L1, R3L2, and ρ=103 for the normal equations. MSPAI is symmetrized via M+MT.
(a)
Noise level 0.001%(b)
Noise level 0.01%Let us focus on theblur [24] test problem as a second example which is deblurring images degraded by atmospheric turbulence blur. The matrix G is an n2-by-n2 symmetric, doubly block Toeplitz matrix that models blurring of an n-by-n image by a Gaussian point spread function. The parameter σ controls the width of G and thus the amount of smoothing and ill-posedness. G=A1⊗A1 is symmetric block banded and possibly positive definite depending on n and σ. We choose the problem to be of size n=150 with bandwidth 4 and σ=1, that is, we invoke blur(150,4,1) and G is of size 1502×1502. The original data vector is denoted by p.We compare the unpreconditioned CG to PCG both forG and for the normal equations GTG. In case of preconditioning G, MSPAI is symmetrized via M+MT and for GTG the preconditioner MTM is applied. For MSPAI, we impose the blocktridiagonal pattern m̂k:=(0,mk-n,k,0|mk-1,k,mk,k,mk+1,k|0,mk+n,k,0)T. In view of the structure of G, we build the high-frequency subspace by Kronecker products of oscillatory probing vectors in the regularizing global conditions R5G1 and the new ones R5G2:min𝒫(M)=𝒫∥(eN1⊗eN2)TM∥2 and R5G3:min𝒫(M)=𝒫∥(eN2⊗eN1)TM∥2, all weighed with ρ=1.Following Figure8 and Table 5, we obtain better reconstruction p̃ when applying MSPAI in contrast to the unpreconditioned CG for G or GTG. We observed similar behavior for other values of band and σ.Table 5
Optimal reconstruction error∥p-p̃∥2 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.
Regularization methodNoise level 0.01%Noise level 0.1%∥p-p̃∥2at it.∥p-p̃∥2at it.CG(G)25.637341.8281CG(GTG)24.1231937.8072PCG(G), MSPAI with R5G1,G2,G323.620639.0071PCG(GTG), MSPAI with R5G1,G2,G323.7785138.0233Regularization methodNoise level 1%Noise level 10%∥p-p̃∥2at it.∥p-p̃∥2at it.CG(G)591.0113.820·1031CG(GTG)82.511.253·1031PCG(G), MSPAI with R5G1,G2,G3189.814.104·1031PCG(GTG), MSPAI with R5G1,G2,G369.710.962·1031Reconstruction error against the iteration for theblur 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 R5G1, R5G2, R5G3, and PCG using MSPAI for normal equations with ρ=1. We use M+MT as preconditioner within PCG(G) and MTM for PCG(GTG).
(a)
Noise level 0.01%(b)
Noise level 1%As a last example we reduce the 2Dblur example of Hansen to the 1D case. For the blur operator we take the 1D analogon of G and reduce the 2D right hand side p to appropriate size. The consideration of example H14 with original data x2 (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 M=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 jth row (0,…,0,rj-1,1,rj,0,…,0). Here, rj≈1/2 near continuous components xj, but rj≈0 near discontinuities. In case that we know the original data x we define Mr with jth row (53)(Mr)j,::=(0,…,0,rj-1,1,rj,0,…,0),forr:=(12⋅11+ρ|xj-xj+1|)j=1,…,n-1.
Otherwise we define the preconditioner Mr̃ with jth row (54)(Mr̃)j,::=(0,…,0,r̃j-1,1,r̃j,0,…,0),viar̃:=(12⋅11+(ρ|bj-bj+1|)k)j=1,…,n-1.
The parameters ρ and k have to be chosen in such a way that discontinuities are revealed as good as possible.Following Figure9, the reconstruction is strongly improved if the preconditioner is adjusted relative to the discontinuities of x. Also, using only the observed data improves the reconstruction. Therefore, we could also consider an iterative process, where a first approximation x1 on x is used to define the tridiagonal preconditioner M1 delivering an improved approximation x2 which again gives a new preconditioner M2, and so on.Reconstruction error against iterations for the 1Dblur problem of size n=103 and band=4 affected by random noise of order 0.1%. G has σ=1 (a) and σ=10 (b), respectively. CG is compared to PCG using Mr and Mr̃, respectively, for G and GTG.
(a)
G with σ=1(b)
G with σ=10


## 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 columnmk of the preconditioner M and can be written in the form minm̂k|ŝkm̂k-fk| with reduced vectors ŝk, m̂k, and scalar fk∈ℝ. Mask probing has the advantage that for each sparse column mk we can use a different sparse probing vector sk. Furthermore, we have introduced probing conditions based on probing vectors that are global for the whole matrix M in the form min𝒫(M)=𝒫∥eTM-fT∥2 and min𝒫(M)=𝒫∥eTAM-fT∥2, respectively, with e,f∈ℝn. The probing vectors are related to the low-frequency or smooth subspace, represented by eS=(1,1,…,1)T, or to the high-frequency oscillatory subspace, represented, for example, by eN1=(1,-1,1,-1,…)T, eN2=(1,0,-1,0,1,…)T, or eN3=(0,1,0,-1,0,…)T.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 onA, only satisfying the approximation eTAM≈eT, 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 vectorx. So, the deterioration of the reconstruction at discontinuities of x 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.



---

*Source: 930218-2010-12-09.xml*