Arnoldi Projection Fractional Tikhonov for Large Scale Ill-Posed Problems

,Xu

1.College of Science,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,P.R.China;

2.Institute of Automation,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,P.R.China

Abstract:It is well known that Tikhonov regularization in standard form may determine approximate solutions that are too smooth for ill-posed problems,so fractional Tikhonov methods have been introduced to remedy this shortcoming.And Tikhonov regularization for large-scale linear ill-posed problems is commonly implemented by determining a partial Arnoldi decomposition of the given matrix.In this paper,we propose a new method to compute an approximate solution of large scale linear discrete ill-posed problems which applies projection fractional Tikhonov regularization in Krylov subspace via Arnoldi process.The projection fractional Tikhonov regularization combines the fractional matrices and orthogonal projection operators.A suitable value of the regularization parameter is determined by the discrepancy principle.Numerical examples with application to image restoration are carried out to examine that the performance of the method.

Key words:ill-posed problems;fractional matrix;Tikhonov regularization;orthogonal projection operator;image restoration

0 Introduction

This paper is concerned with the solution of least-square problem

with a large square matrix A of ill-determined rank.In particular,such a matrix is severely illconditioned and may be singular by which its singular values decrease to zero gradually and without obvious interval.The vector b represents the available data that is usually with a discrete error or measurement error e∈Rn,i.e.

In view of the ill-condition of A and the error e in b,the straightforward solution generally yields a meaningless approximation,so it is essential that the computation is stabilized by regularization.Tikhonov regularization is one of the most popular regularization methods for properties and application.Based on Tikhonov regularization,we consider a penalized least-squares problem

where the scalarμ＞0 is referred to the regularization parameter and the matrix L∈Rl×nis the regularization operator[2-3].The method of this paper requires L to be a square matrix.Calvetti et al.[4]and Hansen et al.[5]described a variety of square regularization operators.For the purpose of obtaining an accurate approximate solution of,the least-squares problem(1)is replaced by the minimization problem(4).The number of rows in L,l≤n,but regularization matrices with l＞n were also applied.Let R(K)and N(K)denote the range and null space of the matrix K,respectively.The matrices A and L to be chosen are assumed to satisfy

Then the Tikhonov minimization problem(4)has the unique solution

for anyμ＞0 and the superscript“Τ”denotes transposition of the matrix[6].

This paper solves the minimization problem(4)by simplifying it to standard form as well as uses a fractional power of the matrixas weighting matrix to measure the residual error in standard with a semi-norm.And then,using a few steps of the Arnoldi process,this paper reduces the problem(3)to a problem of smaller size,which is solved by using the projection fractional Tikhonov,and the regularization parameters a andμare determined.At last,the illustrative numerical examples are also reported,and concluding remark can be found.

1 Projection Fractional Tikhonov

In this section,we discuss the method which combines the fractional matrices and orthogonal projection operators.Projection fractional Tikhonov regularization provides that the penalized least-squares problem(4)can be simplified to standard form and uses a fractional power as weighting matrix to measure the residual error in standard with a semi-norm.

1.1 Form simplification

The penalized least-squares problem(4)can be simplified to standard form with the orthogonal projection

which is well suited for using in Tikhonov regularization.In Eq.(6),L is used as regularization operator.It is convenient to consider the relation of the choice of the matrix L and the matrix P,and actually the choice of P determines the choice of L.Moreover,the choice of matrix P can be carried out in many different ways,some of which may yield regularization operators,and they can give more accurate approximations ofthan the general finite difference-based regularization operators[7].

Give the A-weighted pseudo-inverse of L as

where L?∈ Rn×ldenotes the Moore-Penrose pseudoinverse of the regularization operator L,and I is the identity matrix.

Suppose that Eq.(6)holds and introduce the QR-factorization shown as

where R∈Rl×lis upper triangular and Q∈Rn×lhas orthonormal columns.Using the properties of the Moore-Penrose pseudo-inverse and orthogonal projection,we have the following identities for L

So yield that

Substituting Eqs.(8),(10)into Eq.(7),we get

which simplifies to

Transforming the matrix and vectors of Tikhonov minimization problem(4)by the following substitutions

where

When L is an orthogonal projection operator,Eqs.(13),(14),can be expressed in a simple manner as

An attractive property of this transformation is that thedefined by Eq.(7)is of simple form which makes the orthogonal projection(6)easy to use.For anyμ＞0,letλ=1/μ,and then the minimization problem of Eq.(18)is

Given anyλ＞0,x(λ)has a certain value and is satisfied as

Then

is defined.Consequently,F(λ)is continuous in[0,∞),and some properties of F(λ)are given in the following.

Proposition 1 F(λ)is infinitely differentiable,and has the following properties:

(2)For anyλ＞0,the first and second order derivatives of F(λ)are as follows

Proof:

(1)Computing the inner product of the formula(20)with x(λ)yields

which implies that

According to this estimate and Eq.(22),we obtain that

Thus the conclusion(1)can be drawn from the definition of F(λ).

(2)Implicit differentiation of Eq.(19)with respect toλcombining with Eq.(20)yields

thus the conclusion(2)is proved.

Proof:We consider

Computing the inner product of the formula(23)with x′(λ)yields

In view of Eq.(20),we obtain that

Then we prove that the equal-sign in the above equation does not hold.Assume that＞0 satisfies,then we have.Due to Eq.(23),is obtained,then note that the form(20)yields

Proposition 3 F(λ)satisfies the differential relationship

Proof:Implicit differentiation of Eq.(20)with respect toλyields

Computing the inner product of the above equation with x(λ)yields

and combining with Eq.(20)yields

i.e.

Therefore,Proposition 3 has been proved.

1.2 Fractional Tikhonov

In this section,we use a fractional power of the matrixas weighting matrix to measure the residual error in standard form(18)with a semi-norm[8].We will replace the penalized leastsquares problem(18)by a minimization problem of the form

where the matrix H is symmetric positive semidefinite and

for any M.It is quite natural that the value ofμ counts for a great deal that determines how sensitive the solution of Eq.(24)is to the error e in.The minimization problem(24)has a unique solutionfor anyμ,such as the penalized leastsquares problem(18).

Assuming that

for a＞0.When a＜1,we define H as the Moore-Penrose pseudo-inverse of.The choice of a is the key to determine,which makes the approximate solution more accurate.We refer to the minimization problem(24)as the fractional Tikhonov method(the weighted Tikhonov method)[9].When a=1,we can obtain the standard Tikhonov regularization.

The normal equation associated with the penalized least-squares problem(24)is given by

Then introduce the singular value decomposition(SVD)of,shown as

where

and

are orthogonal matrices and

whose diagonal elements are arranged in the following order

where the index r is the rank of.

Substituting the singular value decomposition(Eq.(28))into Eq.(27)yields

Then the solution of Eq.(27)can be written as

which is equivalent to

where

The solution xμof Eq.(5)can be recovered from the solution of Eq.(33)according to

In addition,the filter function for some a＞0 is given by Eq.(34),it has the following asymptotics

and

Then we consider the filter function of standard Tikhonov regularization shown as

It is easy to show that the filter function(34)is less smoothing than(σ)for 0＜a＜1,and the singular values are damped less by the filter function(34)than by(σ),which means that the approximate solution(35)has higher quality than that with the exact solution.

2 Arnoldi-Projection Fractional Tikhonov

The regularization method is based on the singular value decomposition of the coefficient matrix.However,the singular value decomposition requires a very large amount of computation for the large-scale matrix.Therefore,we choose to project the large-scale problem to the low-dimensional Krylov subspace.Lewis and Reichel proposed Arnoldi Tikhonov regularization method[11]in 2009,and introduced the method in detail.Moreover,Global Arnoldi Tikhonov and Augmented Arnoldi Tikhonov Regularization Methods were successively proposed[12-13].

We propose to reduce the problem(3)to a problem of smaller size by application of the Arnoldi process applied to A with initial vectorThis yields the decomposition

where Vk=[v1,v2,…,vk]∈Rn×kis the first k columns of Vk+1,and Vk+1∈Rn×(k+1)has orthonormal columns,which span the Krylov subspace

We assume that k is chosen sufficiently small so thatis an upper Hessenberg matrix with nonvanishing subdiagonal entries.Thenis of rank k.We seek to determine an approximate solution xμ,kof Eq.(4)in the Krylov subspace(39).

Substituting

into Eq.(4)and using Eq.(38)yields the reduced minimization problem

whose solution is denoted by yμ,k.And the reduced minimization problem(40)solved using the projection fractional Tikhonov regularization methods is described in Section 1,then

is an approximate solution of Eq.(4).

3 Parameters Selection

This section discusses the determination of the regularization parameter.We first consider the effects of parameters a andμon.It follows from the solution that

and

Conjugating

we have

and we assume that an estimate of the norm of the error

Then we can apply the discrepancy principle to determine a suitable value of the regularization parameterμ.Let a＞0 be fixed and define that

whereη＞1 is a user-supplied constant independent ofε.We determineμ＞0,so that the solution xμof Eq.(4)satisfies

Then the vector xμis asked to satisfy the discrepancy principle[15].Solution of Eq.(48)about μis equivalent to the positive zero of the function

where r is the rank of A.Thus

and

We consider the initial approximate solution μ0:=0 for Newton method withμ=μ0-φ′a(μ)/φ″a(μ)to compute the positive zero of the functionφa(μ).The iterations with Newton′s method are terminated as soon as a value ofμ,such that

has been determined.The factor 1/100 in Eq.(52)is used in our implementation,but other positive factors strictly smaller than 1 can be also used.

4 Numerical Examples

We use three text examples to illustrate the performance of the Arnoldi projection fractional Tikhonov(APFT)regularization and compare them to Arnoldi fractional Tikhonov(AFT)and Arnoldi Tikhonov(AT)for large scale linear dis-crete ill-posed problems.The orthogonal projection with

has the same null space as the regularization operator

which will be applied in the following examples.All computations were carried out in MATLAB with about 16 significant decimal digits.

Example 1 Considering the Fredholm integral equation of the first kind shown as

the MATLAB code Shaw produces a discretization A∈R1000×1000and the right-hand side∈R1000by a Galerkin method with orthonormal box functions[16].The noise-levelλis defined by.Then,we will give a comparison of the approximate solution by the APFT regularization method and exact solution when taking the different value of the error vector e.

Fig.1 illustrates that the approximate solution obtained by the APFT method can approximate the exact solution well,which means that APFT regularization method is effective.

Example 2 The Fredholm integral equation of the first kind is

and the MATLAB code discretes Barrt,Shaw,Phillips,Gravity,Foxgod and Deriv2 by a Galerkin method with orthonormal box functions about the matrix order n=1 000.The noise-levelλis defined by

The regularization parameterμis determined by the discrepancy principle.The tables report relative errorsfor several noiselevel and show that the method we proposed improves the accuracy of the computed solutions.

Fig.1 Recovery results of Phillips with diverse noise-level

Tables 1 and 2 show the qualities of AT,AFT and APFT for various examples(n=1 000).The following results show that APFT usually renders solutions of high quality.In other words,we can see that APFT is superior to AFT and AT.

Table 1 Qualities of these methods with the error-level(λ=1%)

Table 2 Qualities of these methods with the error-level(λ=10%)

Example 3 We show the performance of the method about the restoration of a discrete image which has been contaminated by blur and noise.Our task is to deblur the two-dimensional images degraded by additive noise and spatially invariant blur.The restoration problems were proposed by the US Air Force Phillips Laboratory.The twodimensional image restoration problem can be modeled by a linear system of equations Ax=b.The matrix A is a discrete blurring operator referred to as a discrete point spread function.Then the components of the vectors b and^x are the lexicographically-ordered pixel values of distorted images and the exact,respectively.We efficiently compute matrix-vector products without explicitly forming A by using the fast discrete Fourier transform and the discrete point spread function.

Fig.2 displays the noise-and blur-free images,the contaminated image,as well as restored images of Lena which determined by the AFT and APFT methods.Meanwhile,the images above illustrate that APFT gives better reconstructions than AFT.

Fig.2 Original,blurred,and restored Lena images

Fig.3 displays the noise-and blur-free images,the contaminated image,as well as restored images of“MATH”which are determined by the AT and APFT methods.The approximate solutions abtained by the APFT method are nearly optimal for this example.Actually,the computed solutions are close to the orthogonal projection of the exact solution into the range-restricted subspace.However,the AT produces an approxi-mate solution of lower quality than the APFT method.

Fig.3 Original,blurred,and restored MATH images

5 Conclusions

In this paper,we propose the APFT regularization method for solving the large scale linear discrete ill-posed problems.Our method is easy to realize and numerical examples show that the proposed method is effective by which we can give a more accurate approximation than AT and AFT methods.

Acknowledgements

This work was supported by the National Natural Science Foundations of China(Nos.11571171 and 61473148).