A wavelet multiscale iterative regularization ... - Semantic Scholar
Neurocomputing 104 (2013) 138–145
Contents lists available at SciVerse ScienceDirect
Neurocomputing journal homepage: www.elsevier.com/locate/neucom
A wavelet multiscale iterative regularization method for the parameter estimation problems of partial differential equations Hongsun Fu a,n, Bo Han b,n, Hongbo Liu c,n a
Department of Mathematics, Dalian Maritime University, Dalian 116026, China Department of Mathematics, Harbin Institute of Technology, Harbin 150001, China c School of Information, Dalian Maritime University, Dalian 116026, China b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 12 March 2012 Received in revised form 19 August 2012 Accepted 8 October 2012 Communicated by A. Abraham Available online 12 November 2012
A wavelet multiscale iterative regularization method is proposed for the parameter estimation problems of partial differential equations. The wavelet analysis is introduced and a wavelet multiscale method is constructed based on the idea of hierarchical approximation. The inverse problem is decomposed into a sequence of inverse problems which rely on the scale variables and are solved successively according to the size of scale from the longest to the shortest. By combining multiscale approximations, the problem of local minimization is overcomed, and the computational cost is reduced. At each scale, based on the wavelet approximation, the problem of inverting the parameter is transformed into the problem of estimating the finite wavelet coefficients in the scale space. A novel iterative regularization method is constructed. The efficiency of the method is illustrated by solving the coefficient inverse problems of one- and two-dimensional elliptical partial differential equations. & 2012 Elsevier B.V. All rights reserved.
Keywords: Wavelet multiscale Iterative regularization method Partial differential equation Parameter estimation
1. Introduction Inverse problems of partial differential equations (PDEs) originate from various practical problems, and belong to the category of cross subjects. Its meaningful to carry out studies on the theory of inverse problems and their applications [1–3]. Generally speaking, to solve inverse problems is very difficult, the main reason lies in: Ill-posedness: The solution does not depend continuously on the observed data. A minor disturbance of the observed data may cause large change on the solution of the inverse problems. Nonlinearity: Nonlinear dependence of the observation with respect to the parameter to be estimated causes the presence of numerous local minimum, a good initial estimate is crucial for any numerical methods. Large computational cost: An inversion process needs tremendous forward computations, the computational cost is usually very large. Therefore, the key problem of an inverse problem is how to quickly find a stable solution in a wide range.
n
Corresponding authors. E-mail addresses: [email protected] (H. Fu), [email protected] (B. Han), [email protected] (H. Liu). 0925-2312/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.neucom.2012.10.007
In order to overcome the three main difficulties, we investigate an inversion strategy, called wavelet multiscale iterative regularization method, for the general parameter estimation problem of partial differential equation in time–space domain. This strategy synthesizes the property of easy computation of the iterative method with the advantages of minor dependence on the initial model, fast convergence, and stability of results of multiscale inversion. We first decompose the parameter to be inverted by a wavelet multiscale decomposition algorithm, then the original problems are transformed into a sequence of inverse problems depending on the scale variables. At each scale, based on the wavelet approximation, the problem of estimating the parameter in a differential equation is reduced to the problem of estimating the finite coefficients in the approximation spaces. Therefore, when the scale is coarsened for one step, the scope of the inverse problem will be reduced to one-half of it. Consequently the number of unknown coefficients is greatly reduced, and the speed of computation is quickened. In the whole computation from the coarsest scale to the finest one the terminative value of a former scale is taken as an initial estimate for the next finer scale. In this way, the risk of trapping in a local minimum is reduced, and the computational cost may be saved as well. Besides, in order to overcome the ill-posedness, Tikhonov regularization and the perturbation method are combined into a new and stable iterative method, which may be used as an optimistic algorithm to search out the global optimal solution at each scale. The rest of this paper is organized as follows. Related works about the inverse problems of partial differential equations are
H. Fu et al. / Neurocomputing 104 (2013) 138–145
139
reviewed in Section 2. The relevant properties of wavelet are summarized in Section 3. We present the wavelet multiscale inversion method and the iterative method for the general parameter estimation problem of differential equations in Section 4. In order to test our method’s effectiveness, numerical simulations are carried out for inverse problems of recovering coefficients in elliptic differential equations in Section 5. Finally, in Section 6 some conclusions are given.
by combining the wavelet multiscale analysis and the Gauss–Newton iteration. Lu et al.[21] introduced and investigated a numerical method for the treatment of inverse problems based on an adaptive wavelet-Galerkin method. Qian [22] considered a Cauchy problem for an elliptic equation with variable coefficients, and provided a regularization method—wavelet dual least squares method to solve the problem. All these works showed the effectiveness of wavelet on inverse problems.
2. Related works
3. Wavelet
Multiscale inversion is a fairly developed inversion strategy to accelerate convergence, enhance stability of inversion, overcome disturbance of local minima, and by which one can search out the global minimum [4]. The main essence of multiscale inversion is to decompose the original problem into different scales or frequency bands, and to carry out inversion beginning at the longest scale. Since at a longer scale, the objective function shows stronger convexity and has less minima. Therefore, one has a good chance of finding the global minimum. The global minimum of the longest scale is in the neighborhood of the global minimum for the problem at the second longest scale, thus it can be used as the initial guess for the second optimization problem. Successively fining upward in this way succeeds in finding the global minimum of the original fine scale problem [5]. Wavelet theory [6,7] was originally applied as a powerful tool for signal and image processing. Wavelet analysis allows to obtain an efficient sparse representation of some function thanks to a particularly good localization both in space and in scale of each element of a given wavelet basis. In fact, in the wavelet expansion of a function many coefficients are negligible and by discarding these coefficients one can obtain a sparse but accurate approximate representation. In the last years, the good features of wavelet have generated a huge interest in different areas of applied mathematics, physics and engineering. More recently wavelets have been applied to the numerical solution of partial differential [8–11]. Based on the Galerkin principle, a differential equation is solved in the wavelet coefficient space. At present, wavelet methods have already been used for the treatment of inverse problems. For instance, Dicken et al. [12] proposed a method based on wavelet function combined optimal convergence rates with algorithmic efficiency. The proof in that paper utilized the approximation properties of wavelets and results from the general theory of regularization methods. Dahlke et al. [13] discussed some ideas on how adaptive wavelet schemes can be applied to the treatment of certain inverse problems, and showed the main essence of the method through an application to the inverse heat equation. Cohen et al. [14] introduced numerical methods for the treatment of linear inverse problems, based on an adaptive wavelet Galerkin discretization. Wavelet analysis is a multiresolution method which has the properties of adaptability and localization. This allows the idea of wavelet multiresolution to be used for inverse problems. Liu [15] numerically studied wavelet multiresolution method for distributed parameter estimation, and a simple elliptic problem was taken as a suitable model. Several authors have also used wavelet to analyze linear ill-posed problems [16,17]. Recently, wavelet has been applied to solving nonlinear ill-posed problems. For instance, Ghanem et al. [18] presented a wavelet-based approach for identifying the mechanical parameter of zero-memory non-linear discrete structural systems. The presented approach allows both the parameter estimation of a prior known dynamical models and the identification of classes of suitable non-linear models based on input–output data. Chiao et al. [19] presented a wavelet multiresolution procedure for identifying the parameter of non-linear geophysics. Fu et al. [20] consider the problem of estimating the velocity in a two-dimensional acoustic wave equation, and designed a widely convergent method
In this section we summarize the main properties of wavelet as far as they are needed in the following section. For a more detailed introduction to the theory of wavelet we refer to [6,12,23]. 3.1. Wavelet bases Wavelet bases have been documented in numerous textbooks and survey papers (see [11] for a general treatment). With a little effort they can be adapted to fairly general domain O Rd (see [24] for a discussion of the characterizations of function spaces on O by wavelet coefficients). A wavelet basis consists of two type of functions: scaling function fl and wavelet function cl . The index concatenates the usual scale and space parameters j and k. Thus for standard wavelet bases on R, it is defined in Definition 1. Definition 1 (Orthogonal wavelet). A function c A L2 ðRÞ is called an orthogonal wavelet if the set fcj,k ðxÞ ¼ 2j=2 cð2j xkÞg
ð3:1Þ 2
forms an orthonormal basis for L ðRÞ. A wavelet basis on R is labeled by two indices, the shift-index k and the scale-index j. If we collect all basis functions on the same scale in the subspace W j ¼ spanfcj,k 9k A Zg, then we obtain the orthogonal multi-scale decomposition L2 ðRÞ ¼ W j : jAZ
This type of multi-resolution analysis was first introduced by Mallat (see, e.g., [25]). The subspaces Wj are called the scale of the decomposition, and the projection X /f , cj,k SL2 cj,k Q jf ¼ is said to represent the details of f on the scale j. Definition 2 (Multi-resolution analysis). A multi-resolution analysis is a sequence of closed subspaces Vj in L2 ðRÞ, called scaling spaces, satisfying (1) (2) (3) (4) (5) (6)
V j V j þ 1 , for all j A Z, S V is dense in L2 ðRÞ, Tj A Z j j A Z V j ¼ f0g, f ðxÞ A V j if and only if f ð2xÞ A V j þ 1 , 8j A Z, f ðxÞ A V 0 if and only if f ðxkÞ A V 0 , 8k A Z, there exists f A V 0 such that ff0,k : kA Zg
is an orthonormal basis in V0. The function is called the scaling function of the multi-resolution analysis, simply called ‘‘MRA’’. Obviously, the subspaces Vj are spanned by V j ¼ spanffj,k ðxÞ ¼ 2j=2 fð2j xkÞ9k A Zg:
ð3:2Þ
140
H. Fu et al. / Neurocomputing 104 (2013) 138–145
Let us denote Vj ¼ Wl lrj
the subspace of L2 ðRÞ which contains all details of size j or coarser one. However, the notation cl takes into account the possible adaptations of wavelet to multivariate bounded domains in which case the functions cl and fl usually change form near the boundary. At the scale level j which corresponds to resolution 2j , the scaling functions ðfl Þl A Lj span a space Vj within a hierarchy V 0 V 1 L2 ðOÞ of nested approximation spaces. The space Vj should be thought of a finite element space Vh with uniform mesh size h 2j , and local basis ðfl Þl A Lj . The wavelets ðcl Þl A Lj span a complement Wj of Vj into V j þ 1 . With these notations, the wavelet decomposition takes the form X XX f¼ cl fl þ dl cl , ð3:3Þ
on some fine scale J. Then employing the scaling Eqs. (3.4) and (3.5), we obtain recursively X cjk ¼ hn2k cjkþ 1 , ð3:8Þ mAZ j
dk ¼
X
g n2k cjkþ 1 :
ð3:9Þ
mAZ
It remains to compute the values of cJk on a fine scale J. In a realistic situation we are only given discrete values f k , k A Z, of f than the function itself. Not all scaling functions have the interpolating property like the Haar function. But they necessarily R have a non-vanishing mean value [16], fðxÞ dx a 0, hence J J 2 fð2 xÞ approximates the d-distribution, and we simply put f k ¼ cJk on an appropriate scale J.
j Z l l A ,j
l A Ll
where ðfl Þl A Ll is the scaling function basis spanning the approximation at level l, ðcl Þl A Lj is the wavelet basis spanning the details at level j, and cl ¼ cl ðf Þ and dl ¼ dl ðf Þ are wavelet coefficients of f, respectively.
4. Inversion strategies 4.1. Model and wavelet multiscale inversion strategy Consider the following partial differential equation defined in the time–space domain:
3.2. Multiscale decomposition and reconstruction
LðcðxÞ,tÞuðx,tÞ ¼ 0,
For simplicity, we illustrate multiscale decomposition and reconstruction only for one dimension. The case of multidimension is analogous. Let fV j ðxÞ; j A Zg be a one-dimensional multiresolution analysis (MRA), its scale function is fðxÞ, the wavelet function is cðxÞ. Since f A V 0 V 1 and c A W 0 V 1 , it follows that both the wavelet c and scaling function f satisfy the scaling equations pffiffiffi X fðxÞ ¼ 2 hk fð2xkÞ, ð3:4Þ kAZ
pffiffiffi X
cðxÞ ¼ 2
g k fð2xkÞ,
ð3:5Þ
kAZ
T
fðxÞ ¼ fð2xÞ þ fð2x1Þ, cðxÞ ¼ fð2xÞfð2x1Þ:
Euðx,0Þ ¼ gðxÞ, x A O:
Buðx,tÞ ¼ f ðx,tÞ,
x A @O, 0 r t rT,
x A G @O, 0 rt r T,
cj ðxÞ ¼ P j ðcðxÞÞ:
j
ð3:6Þ
is particularly easy and fast. Assume that we know the values of the scalar products cJk ¼ /f , fJ,k SL2
ð3:7Þ
ð4:4Þ
with A the auxiliary condition operator and G a part of @O, (4.1)–(4.4) constitutes the general parameter estimation problem of partial differential equations. Let c A L2 ðRÞ is a wavelet belonging to the Daubechies family of compactly supported wavelets, f A L2 ðRÞ is the corresponding scale function. For j A Z, let
g k ¼ ð1Þk h2N1k ,
dk ¼ /f , cj,k SL2 ,
ð4:3Þ
where E, B are the initial condition operator and boundary condition operator, respectively. If c(x), g(x), f ðx,tÞ are known, (4.1)–(4.3) constitutes a forward problem to determine uðx,tÞ. If c(x) is unknown, given the auxiliary condition
V j ¼ spanffj,i ¼ 2j=2 fð2j xiÞ9i A Zg,
For actual computations one is only interested in wavelets and scaling functions which obey such finite scaling equations. In this case the computation of the wavelet decomposition, i.e., the computation of the scalar products
ð4:2Þ
The boundary condition is
The Haar wavelet is the only obvious compactly supported orthogonal wavelet. It was the accomplishment of Daubechies [6] to construct a family ffN g of orthogonal wavelets with suppðfN Þ ¼ ½0,2N1 and linearly increasing regularity. Here c1 equals the Haar wavelet. The corresponding filter coefficients gk in (3.4) are given by k ¼ 0,1, . . . ,2N1:
ð4:1Þ
p
where x ¼ ðx1 ,x2 , . . . ,xp Þ , O R , @O is the boundary of O, uðx,tÞ is a sufficiently smooth function defined in O ð0,TÞ. L is a differential operator. The initial condition is
Auðx,tÞ ¼ jðx,tÞ,
with real scaling coefficients hk, gk. For example, the Haar function fðxÞ ¼ w½0,1 ðxÞ and the Haar wavelet cðxÞ ¼ w½1,1=2 ðxÞw½1=2,1 ðxÞ satisfy the above conditions with
x A O, 0 r t r T,
W j ¼ spanfcj,i ¼ 2j=2 fð2j xiÞ9i A Zg: Suppose Pj is the projector from L2 ðRÞ to Vj, Pj : L2 ðRÞ-V j , and denote the projection of cðxÞ A L2 ðRÞ by ð4:5Þ
j
Obviously c ðxÞ A V j , and lim cj ðxÞ ¼ cðxÞ:
j- þ 1
ð4:6Þ
Since, in practical computations, the solving domain is often bounded, and the differential equation, the initial condition, the boundary condition, and the auxiliary condition have to be discretized, the original inverse problem will be approximated by a finite dimensional model, and c(x) only has finite resolution. Therefore, in the inversion process of wavelet multiscale method, the value of scale j is taken in a range 0 rj rJ. The approximation
H. Fu et al. / Neurocomputing 104 (2013) 138–145
where dcjl ðxÞ is a small perturbation, which can be determined by the minimization of the following nonlinear functional:
of the function c(x) can now be written as X
cJ ðxÞ ¼
J X X
c0i f0,i ðxÞ þ
i A L0
j
di cj,i ðxÞ,
ð4:7Þ
j ¼ 0 i A Lj
where Lj is an index set. Therefore, the problem to determine cJ ðxÞ is reduced to solving a finite dimensional wavelet coefficients vector. And according to the principle of multiscale inversion, the undetermined parameter c(x) will be decomposed onto different subspaces Vj, such that the problem of determining c(x) is transformed into the problem of determining the projection cj ðxÞ in the subspaces Vj, j ¼ 0,1, . . . ,J. We start from the coarsest level of resolution and progressively move to the finest level. Let us rewrite approximation (4.7) for any scale j ð0 r j r JÞ as X j cj ðxÞ ¼ cj1 ðxÞ þ di cj,i ðxÞ, 0
c ðxÞ ¼
i A Lj
X
c0i f0,i ðxÞ:
ð4:8Þ
Ja ½dcj ðxÞ ¼ JA½cjl ðxÞ þ dcjl jðx,tÞJ2G½0,T þ aJdcjl ðxÞJ2 :
At each scale, cj ðxÞ will be expressed as an expansion of the wavelet basis, and the problem of determining cj ðxÞ is reduced to determine the finite coefficients of representation. The whole inversion will start at the longest scale. At the longest scale j¼0, use the iterative regularization method to make correction on the initial guess until the global minimum is obtained, and the approximate solution will be used as the initial guess for the inversion at scale j ¼1. Successively fining upward in this way until we finally find a good approximate solution cJn ðxÞ at the finest scale j¼J. cJn ðxÞ can be regarded as the optimal solution of the original problem.
to denote the solution of the forward We use problem (4.1)–(4.3) when cjl ðxÞ is given. The corresponding solution of (4.1)–(4.3) is denoted by uðcjl ðxÞ þ dcjl ðxÞ; x,tÞ after being added a minor perturbation X
dcjl ðxÞ ¼
j X X
T ddj,l i cj,i ðxÞ9dK l FðxÞ
ð4:14Þ
j ¼ 0 i A Lj
to cjl ðxÞ. In this way, the problem to determine dcjl ðxÞ is reduced to solve dK l . Whereas dK l can be determined by the minimization of the following functional: F½dK l ¼ JAuðcjl ðxÞ þ dcjl ðxÞ; x,tÞjðx,tÞJ2G½0,T þ aJdK l J2 :
ð4:15Þ
is a minor perturbation, from the multivariate Since Taylor expansion formula, we have T
uðcjl ðxÞ þ dcjl ðxÞ; x,tÞ ¼ uðcjl ðxÞ; x,tÞ þ XK l uðcjl ðxÞ; x,tÞ dK l þ OðJdcjl ðxÞJ2 Þ,
j X X
c0i f0,i ðxÞ þ
i A L0
T
F½dK l ¼ JAuðcjl ðxÞ; x,tÞjðx,tÞ þ XK l Auðcjl ðxÞ; x,tÞ dK l J2G½0,T þ aJdK l J2 :
If taking 2j discrete points ðxm ,t m Þ ðm ¼ 0,1, . . . ,2J 1Þ, then from numerical integration theory, we have
j
A½c ðxÞ ¼ A2 ½A1 ½cð ðxÞ ¼ jðx,tÞ:
ð4:10Þ
j
2
J a ½c ðxÞ ¼ JA½c ðxÞjðx,tÞJG½0,T þ aJc ðxÞJ ,
ð4:18Þ
J
Generally, it is ill-posed in the sense of Hadamard and should be regularized. Tikhonov regularization [26] of such ill-posed problem is achieved by replacing the operator equation (4.8) by the minimization problem ð4:11Þ
where a is the regularization parameter. We will use the perturbation method and linearization technique to establish an iterative method. Let cjl þ 1 ðxÞ ¼ cjl ðxÞ þ dcjl ðxÞ,
rm ½Auðcjl ðxÞ; xm ,tm Þjðxm ,tm Þ
þ XK l Auðcjl ðxÞ; xm ,t m Þ dK l 2 þ aJdK l J2 ,
ð4:9Þ 0
2
J 2X 1
T
j
di cj,i ðxÞ9K T FðxÞ
in the subspace V j ð0 r j rJÞ, where K T ¼ fðc0i Þi A L0 ; ðdi Þi A L0 ; 1 j ðdi Þi A L1 ; ; ðdi Þi A Lj g is the wavelet coefficient in the approximation subspace Vj, FðxÞ ¼ fðf0,i ðxÞÞi A L0 , ðc0,i ðxÞÞi A L0 ,ðc1,i ðxÞÞi A L1 , . . . , ðcj,i ðxÞÞi A Lj gT is the wavelet basis function on the scale j. Therefore, the problem to determine cj ðxÞ is reduced to solving a finite dimensional vector K. On the other hand, since the solution uðx,tÞ of (4.1)–(4.3) is nonlinearly dependent on cj ðxÞ, if we define a nonlinear operator A1 which maps cj ðxÞ to uðcj ðxÞ; x,tÞ, then A1 ½cj ðxÞ ¼ uðcj ðxÞ; x,tÞ. From (4.4) we know that there exists a nonlinear operator A2 mapping from the solution uðcj ðxÞ; x,tÞ to the data jðx,tÞ such that A2 ½uðcj ðxÞ; x,tÞ ¼ jðx,tÞ. Therefore, the inverse problem is equivalent to solve the following nonlinear operator equation:
j
ð4:17Þ
m¼0
j ¼ 0 i A Lj
j
ð4:16Þ
consequently
F½dK l ¼
Approximation (4.7) can be expanded as
j
dc0,l f0,i ðxÞ þ i
i A L0
4.2. An iterative regularization method
cj ðxÞ ¼
ð4:13Þ
uðcjl ðxÞ; x,tÞ
dcjl ðxÞ
i A L0
X
141
ð4:12Þ
where frm g2m 1 ¼ 1 are some adequate weights. Let 2 3 r0 Auðcjl ðxÞ; x0 ,t 0 Þ 6 7 6 7 r1 Auðcjl ðxÞ; x1 ,t1 Þ 6 7 U¼6 7, 6 7 ^ 4 5 r2J 1 Auðcjl ðxÞ; x2J 1 ,t2J 1 Þ 2 3 r0 jðx0 ,t0 Þ 6 7 r1 jðx1 ,t1 Þ 6 7 7, Un ¼ 6 6 7 ^ 4 5 r2J 1 jðx2J 1 ,t2J 1 Þ H ¼ ðhmi Þ2J dimðLj Þ , hmi ¼
@ r Auðcjl ðxÞ; xm ,tm Þ: @ki m
Then (4.16) can be simplified to F½dK l ¼ dK Tl HT HdK l þ 2dK Tl HT ðUU n Þ þ ðUU n ÞT ðUU n Þ þ aJdK l J2 :
ð4:19Þ
Obviously the minimum point of F½dK l satisfies ðHT H þ aIÞdK l ¼ HT ðUU n Þ:
ð4:20Þ
It is just Euler’s equation of (4.17). Solving dK l from (4.18) and substituting it into (4.12), we can finally obtain the perturbation dcjl ðxÞ, and cjl þ 1 ðxÞ further. Repeat the above process until a satisfied precision is obtained.
142
H. Fu et al. / Neurocomputing 104 (2013) 138–145
if f A H1 ðOÞ for almost all x A O, we denote
4.3. The process of the wavelet multiscale iterative regularization method
X ¼ fcðxÞ A L1 ðOÞ : 0 o g rcðxÞ r g , suppðcðxÞc0 ðxÞÞ Og
In summary, we present the main process of the wavelet multiscale iterative regularization method for the parameter estimation problem of differential equation:
and define F : X -Y ¼ H10 ðOÞ,
In each case the forward operator was discretized by using finite elements on a uniform grid (triangular, in the case of two dimensions).
(1) We first make the multiscale decomposition of c(x) by using Daubechies compact supported orthogonal wavelet to obtain cj ðxÞ ðj ¼ J,J1, . . . ,0Þ. The size of J is determined by the resolution ratio of practical data. Given an initial estimate c0 ðxÞ of the true parameter cn ðxÞ. Thus the initial guess c0 ðxÞ P can be approximated as c0 ðxÞ ¼ i A L0 c0i f0,i ðxÞ. We solve the 0 numerical solution uðcl ; xm ,t m Þ of the problem (4.1) and (4.3) by the finite-difference method, and obtain the value Auðc0l ; xm ,t m Þ ðm ¼ 0,1, . . . ,2J 1Þ in the auxiliary condition. (2) Compute the elements hmi of the matrix H by numerical differentiation
5.1. One-dimensional examples In our first example, we take 8 n c ðxÞ ¼ 10099:9ð12xÞ2 , > > > > > uðxÞ ¼ xð1xÞ, > > > > < f ðxÞ ¼ 200599:4ð12xÞ2 , > g ¼ 0:1, g ¼ 200, > > > > > c > 0 ðxÞ ¼ 100, > > : O ¼ ½0,1:
@ Auðc0l ; xm ,t m Þ @ki 1 ¼ ½Auðc0l ðxÞ þ tf0,i ðxÞ; xm ,t m ÞAuðc0l ðxÞ; xm ,t m Þ:
hmi ¼
(3) By the iterative method given in the above subsection, we can obtain the global minimum point c0n at scale 0. (4) Conduct the process (1)–(3) at scale 1 by taking c0n ðxÞ as an initial value for c1n ðxÞ. (5) Repeat the process (1)–(4), until we finally obtain the global minimum point cJn ðxÞ. 5. Numerical examples
Remark. Because the parameter c(x) that we are looking for is required only to be in O ¼ ½0,1, we shall use the simplest wavelet and scaling function, namely the Haar wavelet and its scaling function.
To show the feasibility of our approach and illustrate its numerical behavior, we consider the elliptical differential equation
r ðcðxÞruÞ ¼ f
in O Rd ,
u¼0
on @O,
The second model is depicted in Fig. 2. In the experiments, we try to identify a piecewise continuous c(x) with several jump discontinuities by using the proposed method. We take (1 x A ½0, 15 [ ½25 , 35 [ ½45 ,1, 2, cðxÞ ¼ ð5:3Þ 1, x A ½15 , 25 [ ½35 , 45:
ð5:1Þ
The wavelet basis and parameters a, h are the same as the first example.
110
110
100
100
90
90
80
80
70
70
c (x)
c (x)
where f A L2 ðOÞ. The observations of the solution u are used to recover the coefficient c(x). Our task is to determine an approximation of parameter c(x) from a measurement ud in the model problem (5.1) for d¼1,2. Since the model problem is only well-posed for cðxÞ Z g 40, and
60 50
60 50
exact solution itertion solution
40
exact solution itertion solution
40
30 20
ð5:2Þ
The solution of the forward problem (5.1) was based on linear finite elements on a uniform grid with mesh size h ¼ 1=26 . In this case, we use wavelet multiscale decomposition method to decompose the parameter c(x) into six scales, thus we obtain a sequence of inverse problems which depend on the scale variables j ¼ 0,1, . . . ,6. At each scale j, we take regularization parameter a ¼ 2j 104 , the iteration will terminate when the iteration number kj r20 or the relative error is equal or less than 1=ð10ðj þ 1ÞÞ 100%. To test the noisy suppression ability of our methods, we invert the parameter c(x) when there is 10% random noise added to the data u. Fig. 1 illustrates the numerical results.
t
(
cðxÞ/u ðsolution of ð5:1ÞÞ:
30 0
0.2
0.4
x
0.6
0.8
1
20
0
0.2
0.4
x
0.6
0.8
1
Fig. 1. Inversion results for the first example in one-dimensional. (a) There is no noise. (b) There is 10% noise added to data.
H. Fu et al. / Neurocomputing 104 (2013) 138–145
3
3 exact solution itertion solution
2.5
2 c (x)→
c (x)→
exact solution itertion solution
2.5
2 1.5
1.5
1
1
0.5
0.5
0
143
0
0.2
0.4
x→
0.6
0.8
1
0
0
0.2
0.4
0.6 x→
0.8
1
Fig. 2. Inversion results for the second example in one-dimensional. (a) There is no noise. (b) There is 10% noise added to data.
Table 1 The numerical results of the wavelet multiscale iterative regularization method with noisy data (J ¼ 6, C ¼ 1). k
d%
JcJl cJn J JcJn J
d%
k
JcJl cJn J JcJn J
1:00 101
58
1:023 100
1:00 100
60
0:834 100
0:50 101
49
0:815 100
0:50 100
55
0:658 100
1:00 100
42
0:364 100
1:00 100
46
0:243 100
0:50 100
39
0:522 101
0:50 101
41
0:276 101
1:00 101
36
0:204 101
1:00 101
39
0:185 101
1
0:50 10
33
0:418 10
2
1
0:50 10
37
0:214 102
0:25 101
32
0:221 102
0:25 101
32
0:102 102
Jcjl cjn J=Jcjn J
In the process of computation, we take ðj ¼ 0, 1,2, . . . ,6Þ as relative error levels. Fig. 2 illustrates the inversion results. Experiment results with noisy data are displayed in Table 1. Here we use uniformly distributed random noise. For each noise level, we carry out four test runs and compute the average relative error JcJl cJn J=JcJn JðJ ¼ 6Þ. The left three columns in Table 1 show the results for the first example with the initial value c0 ðxÞ 100. In the right three columns we display the results for the second example with the initial value c0 ðxÞ 1. Here, k denotes the iteration total number by using the wavelet multiscale iterative regularization method.
5.2. Two-dimensional examples To illustrate the applicability of the proposed method to higher-dimensional parameter identification problems, we also consider the model problem in two-dimension. Here, we take O ¼ ½0,12 , f ðx1 ,x2 Þ 1, and " # " # ðx1 0:5Þ2 ðx2 0:5Þ2 cðx1 ,x2 Þ ¼ exp exp , 0:05 0:05 see Fig. 3(a). Due to the specific geometry of O, we can use periodized Daubechies wavelets of period 1. In particular, the space Vj ¼ V j V j Vj þ 1 ,
jAZ
form a separable two-dimensional MRA with scale function
Fðx1 ,x2 Þ ¼ fðx1 Þfðx2 Þ. The orthogonal complement of Vj in Vj þ 1 is Wj ¼ Whj Wvj Wdj ,
where W0 is defined by the horizontal, vertical and diagonal wavelet functions
Ch ðx1 ,x2 Þ ¼ fðx1 Þcðx2 Þ, Cv ðx1 ,x2 Þ ¼ cðx1 Þfðx2 Þ, Cd ðx1 ,x2 Þ ¼ cðx1 Þcðx2 Þ: The [0,1]-periodic scaling function and wavelets are defined by the formulae XX j Fk1 ,k2 ðx1 l1 ,x2 l2 Þ l1 A Z l2 A Z
and XX
Cj,k1l,k2 ðx1 l1 ,x2 l2 Þ, l ¼ h,v,d,
l1 A Z l2 A Z
respectively. To obtain un and ud , we solve (5.1) by the finite element method on a mesh with N ¼ 64 64 elements. In numerical simulations, we use two-dimensional fast wavelet transform to decompose the parameter c(x) into six scales. We first consider the case when there is no noise, i.e. ud ¼ un . Fig. 3(b) and (c) shows un and the identified coefficients cJn ðx1 ,x2 ÞÞ from the noise-free observation. It shows that our algorithms have reconstructed the parameter c(x) very accurately. Second, we consider noisy data. The noise level is 12.14%. The reconstruction results, shown in Fig. 3(d), approximate well to the true solution. Numerical results indicate (1) when there is no noise in the data, for the two different examples in one dimension, the inversion results have high precisions. Besides, at scale 0, only one coefficient needs to be inverted, and only two iterations are performed to obtain an approximate value with 10% relative error. Although iterative number will increase with the increment of scale, since the final approximation at a former scale will be used as a good initial value for the latter, even at the finest scale—scale 6, only 13 and 15 iterations are needed to obtain the global optimal solution for this two examples, respectively. It indicates the advantages of the wavelet multiscale inversion method on the aspects of computational precision, efficiency and convergent range. (2) When the noise is added, the number of iteration will increase, and the precision of inversion results are lowered a bit. But because wavelet multiscale decomposition has a good ability of noise suppression, the effectiveness of noise will gradually drop. Therefore, we can still obtain almost perfect reconstruction results even 12.14% noise is added in two dimension.
144
H. Fu et al. / Neurocomputing 104 (2013) 138–145
Fig. 3. Inversion results for the 2D example. (a) Exact solution. (b) Exact data u. (c) There is no noise. (d) There is noise added to data.
6. Conclusions For the general parameter estimation problem of partial differential equations, a wavelet multiscale iterative regularization method is proposed in this paper. Numerical simulations for the inversion problem of one- and two-dimensional elliptical equation showed the method’s effectiveness in the aspects of stability, global convergence, noise suppression and fast computation. Wavelet multiscale method widens the scope of research of multiscale inversion. To conduct wavelet multiscale decomposition for the undetermined parameter will lower the dependence of the inverse problems on initial estimate, and overcome the effect of local minimum in the objective functional. The computational cost is greatly reduced by transforming the undetermined parameter into the finite wavelet coefficients. Moreover, because there are ready-made fast algorithms of decomposition and reconstruction, the operation is simple, fast, and stable. It is helpful to solve more complicated inverse problems. However, many theoretical and computational problems are still open.
Acknowledgments This work is partly supported by the National Natural Science Foundation of China (Grant nos. 41074088 and 61173035), the Doctoral Science Technology Foundation of Liaoning Province (Grant no. 20091007), the China Postdoctoral Science Foundation (Grant no. 20110491533), the Fundamental Research Funds for the Central Universities (Grant no. 2012TD027), the Program for New Century Excellent Talents in University (Grant no. NCET-110861).
References [1] V. Isakov, Inverse Problems for Partial Differential Equations, vol. 127, Springer-Verlag, 2006.
[2] B. Choi, J. Lee, Comparison of generalization ability on solving differential equations using backpropagation and reformulated radial basis function networks, Neurocomputing 73 (1–3) (2009) 115–118. [3] D. Luengo, S. Monzon, A. Artes-Rodriguez, Novel fast random search clustering algorithm for mixing matrix identification in mimo linear blind inverse problems with sparse inputs, Neurocomputing 87 (2012) 62–78. [4] D. Benaouda, F. Murtagh, J. Starck, O. Renaud, Wavelet-based nonlinear multiscale decomposition model for electricity load forecasting, Neurocomputing 70 (1–3) (2006) 139–154. [5] M. Mu, Q. Ruan, S. Guo, Shift and gray scale invariant features for palmprint identification using complex directional wavelet and local binary pattern, Neurocomputing 74 (17) (2011) 3351–3360. [6] I. Daubechies, Ten Lectures on Wavelet, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1992. [7] Y. Meyer, D. Salinger, C.U. Press, Wavelets And Operators, vol. 2, Cambridge University Press, 1992. [8] G. Beylkin, On the representation of operators in bases of compactly supported wavelets, SIAM J. Numer. Anal. 6 (6) (1992) 1716–1740. [9] W. Dahmen, Wavelet methods for PDEs—some recent developments, J. Comput. Appl. Math. 128 (1) (2001) 133–185. ¨ [10] J. Frohlich, K. Schneider, An adaptive wavelet—vaguelette algorithm for the solution of pdes, J. Comput. Phys. 130 (2) (1997) 174–190. [11] V. Comincioli, G. Naldi, T. Scapolla, A wavelet-based method for numerical solution of nonlinear evolution equations, Appl. Numer. Math. 33 (1–4) (2000) 291–297. [12] V. Dicken, P. Maaß, Wavelet-Galerkin methods for ill-posed problems, J. Inver. Ill-Posed Prob. 4 (3) (1996) 203–222. [13] S. Dahlke, P. Maaß, An outline of adaptive wavelet Galerkin methods for Tikhonov regularization of inverse parabolic problems, in: Recent Development in Theories and Numerics. Proceedings of the International Conference on Inverse Problems, 2003, pp. 56–66. [14] A. Cohen, M. Hoffmann, M. Reiss, Adaptive wavelet Galerkin methods for linear inverse problems, SIAM J. Numer. Anal. 42 (4) (2005) 1479–1501. [15] J. Liu, A multiresolution method for distributed parameter estimations, SIAM J. Sci. Comput. 14 (2) (1991) 389–405. [16] D. Donoho, Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition, Appl. Comput. Harmonic Anal. 2 (2) (1995) 101–126. [17] T. Reginska, L. Elde´n, Stability and convergence of a wavelet-Galerkin method for the sideways heat equation, J. Inver. Ill-Posed Prob. 8 (3) (2000) 31–49. [18] R. Ghanem, F. Romeo, A wavelet-based approach for model and parameter identification of non-linear systems, Int. J. Non-lin. Mech. 36 (5) (2001) 835–859. [19] L. Chiao, W. Liang, Multiresolution parameterization for geophysical inverse problems, Geophysics 68 (1) (2003) 199–209. [20] H. Fu, B. Han, A wavelet multiscale method for the inverse problems of a twodimensional wave equation, Inver. Prob. Sci. Eng. 12 (6) (2004) 643–656. [21] X. Lu, Y. Sun, G. Bai, Adaptive wavelet-Galerkin methods for limited angle tomography, Image Vis. Comput. 28 (4) (2010) 696–703.
H. Fu et al. / Neurocomputing 104 (2013) 138–145
[22] A. Qian, A new wavelet method for solving an ill-posed problem, Appl. Math. Comput. 203 (2) (2008) 635–640. [23] A. Louis, P. Maaß, A. Rieder, Wavelets: Eine Einfuhrung in Theorie und Anwendungen, Teubner Verlag, Stuttgart, 1994. [24] A. Cohen, Numerical analysis of wavelet methods, Stud. Math. Appl. 32 (2003) 1–335. [25] S. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2 ðRÞ, Trans. Am. Math. Soc. 315 (1) (1989) 69–87. [26] A. Tikhonov, V. Arsenin, F. John, Solutions of Ill-Posed Problems, Winston Washington, DC, 1977.
Hongsun Fu received MS and PhD at Harbin Institute of Technology, China. She is an Associate Professor of the Department of Mathematics, Dalian Maritime University. She serves as a referee for the inverse problems in Science & Engineering.
Bo Han received MS and PhD at the Harbin Institute of Technology, China. He is a Professor of the Department of Mathematics, Harbin Institute of Technology. He is an Executive Director of Society for Industrial and Applied Mathematics in China. He is also the President of Society for Industrial and Applied Mathematics in Heilongjiang Province, China. He serves as a Referee for the Mathematical Reviews of Chinese Mathematical Society. Also he serves as a Referee for International Journals: Multiscale Modeling and Simulation (SIAM); International Journal on Numerical Analysis and Modeling; Selected Publications of Chinese Universities: Mathematics; Journal on Information and Computational Science; Communications in Computational Physics; Journal of Computational Mathematics.
145 Hongbo Liu received his three level educations (BSc, MSc, PhD) at the Dalian University of Technology, China. Currently he is a Professor at School of Information Engineering and Science, Dalian Maritime University, with an affiliate appointment in the Computer Science and Biomedical Engineering Division at the Dalian University of Technology. His research interests are in system modeling and optimization involving soft computing, probabilistic modeling, cognitive computing, machine learning, data mining, etc. He participates and organizes actively International Conference and Workshop and International Journals/Publications.
Contents lists available at SciVerse ScienceDirect
Neurocomputing journal homepage: www.elsevier.com/locate/neucom
A wavelet multiscale iterative regularization method for the parameter estimation problems of partial differential equations Hongsun Fu a,n, Bo Han b,n, Hongbo Liu c,n a
Department of Mathematics, Dalian Maritime University, Dalian 116026, China Department of Mathematics, Harbin Institute of Technology, Harbin 150001, China c School of Information, Dalian Maritime University, Dalian 116026, China b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 12 March 2012 Received in revised form 19 August 2012 Accepted 8 October 2012 Communicated by A. Abraham Available online 12 November 2012
A wavelet multiscale iterative regularization method is proposed for the parameter estimation problems of partial differential equations. The wavelet analysis is introduced and a wavelet multiscale method is constructed based on the idea of hierarchical approximation. The inverse problem is decomposed into a sequence of inverse problems which rely on the scale variables and are solved successively according to the size of scale from the longest to the shortest. By combining multiscale approximations, the problem of local minimization is overcomed, and the computational cost is reduced. At each scale, based on the wavelet approximation, the problem of inverting the parameter is transformed into the problem of estimating the finite wavelet coefficients in the scale space. A novel iterative regularization method is constructed. The efficiency of the method is illustrated by solving the coefficient inverse problems of one- and two-dimensional elliptical partial differential equations. & 2012 Elsevier B.V. All rights reserved.
Keywords: Wavelet multiscale Iterative regularization method Partial differential equation Parameter estimation
1. Introduction Inverse problems of partial differential equations (PDEs) originate from various practical problems, and belong to the category of cross subjects. Its meaningful to carry out studies on the theory of inverse problems and their applications [1–3]. Generally speaking, to solve inverse problems is very difficult, the main reason lies in: Ill-posedness: The solution does not depend continuously on the observed data. A minor disturbance of the observed data may cause large change on the solution of the inverse problems. Nonlinearity: Nonlinear dependence of the observation with respect to the parameter to be estimated causes the presence of numerous local minimum, a good initial estimate is crucial for any numerical methods. Large computational cost: An inversion process needs tremendous forward computations, the computational cost is usually very large. Therefore, the key problem of an inverse problem is how to quickly find a stable solution in a wide range.
n
Corresponding authors. E-mail addresses: [email protected] (H. Fu), [email protected] (B. Han), [email protected] (H. Liu). 0925-2312/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.neucom.2012.10.007
In order to overcome the three main difficulties, we investigate an inversion strategy, called wavelet multiscale iterative regularization method, for the general parameter estimation problem of partial differential equation in time–space domain. This strategy synthesizes the property of easy computation of the iterative method with the advantages of minor dependence on the initial model, fast convergence, and stability of results of multiscale inversion. We first decompose the parameter to be inverted by a wavelet multiscale decomposition algorithm, then the original problems are transformed into a sequence of inverse problems depending on the scale variables. At each scale, based on the wavelet approximation, the problem of estimating the parameter in a differential equation is reduced to the problem of estimating the finite coefficients in the approximation spaces. Therefore, when the scale is coarsened for one step, the scope of the inverse problem will be reduced to one-half of it. Consequently the number of unknown coefficients is greatly reduced, and the speed of computation is quickened. In the whole computation from the coarsest scale to the finest one the terminative value of a former scale is taken as an initial estimate for the next finer scale. In this way, the risk of trapping in a local minimum is reduced, and the computational cost may be saved as well. Besides, in order to overcome the ill-posedness, Tikhonov regularization and the perturbation method are combined into a new and stable iterative method, which may be used as an optimistic algorithm to search out the global optimal solution at each scale. The rest of this paper is organized as follows. Related works about the inverse problems of partial differential equations are
H. Fu et al. / Neurocomputing 104 (2013) 138–145
139
reviewed in Section 2. The relevant properties of wavelet are summarized in Section 3. We present the wavelet multiscale inversion method and the iterative method for the general parameter estimation problem of differential equations in Section 4. In order to test our method’s effectiveness, numerical simulations are carried out for inverse problems of recovering coefficients in elliptic differential equations in Section 5. Finally, in Section 6 some conclusions are given.
by combining the wavelet multiscale analysis and the Gauss–Newton iteration. Lu et al.[21] introduced and investigated a numerical method for the treatment of inverse problems based on an adaptive wavelet-Galerkin method. Qian [22] considered a Cauchy problem for an elliptic equation with variable coefficients, and provided a regularization method—wavelet dual least squares method to solve the problem. All these works showed the effectiveness of wavelet on inverse problems.
2. Related works
3. Wavelet
Multiscale inversion is a fairly developed inversion strategy to accelerate convergence, enhance stability of inversion, overcome disturbance of local minima, and by which one can search out the global minimum [4]. The main essence of multiscale inversion is to decompose the original problem into different scales or frequency bands, and to carry out inversion beginning at the longest scale. Since at a longer scale, the objective function shows stronger convexity and has less minima. Therefore, one has a good chance of finding the global minimum. The global minimum of the longest scale is in the neighborhood of the global minimum for the problem at the second longest scale, thus it can be used as the initial guess for the second optimization problem. Successively fining upward in this way succeeds in finding the global minimum of the original fine scale problem [5]. Wavelet theory [6,7] was originally applied as a powerful tool for signal and image processing. Wavelet analysis allows to obtain an efficient sparse representation of some function thanks to a particularly good localization both in space and in scale of each element of a given wavelet basis. In fact, in the wavelet expansion of a function many coefficients are negligible and by discarding these coefficients one can obtain a sparse but accurate approximate representation. In the last years, the good features of wavelet have generated a huge interest in different areas of applied mathematics, physics and engineering. More recently wavelets have been applied to the numerical solution of partial differential [8–11]. Based on the Galerkin principle, a differential equation is solved in the wavelet coefficient space. At present, wavelet methods have already been used for the treatment of inverse problems. For instance, Dicken et al. [12] proposed a method based on wavelet function combined optimal convergence rates with algorithmic efficiency. The proof in that paper utilized the approximation properties of wavelets and results from the general theory of regularization methods. Dahlke et al. [13] discussed some ideas on how adaptive wavelet schemes can be applied to the treatment of certain inverse problems, and showed the main essence of the method through an application to the inverse heat equation. Cohen et al. [14] introduced numerical methods for the treatment of linear inverse problems, based on an adaptive wavelet Galerkin discretization. Wavelet analysis is a multiresolution method which has the properties of adaptability and localization. This allows the idea of wavelet multiresolution to be used for inverse problems. Liu [15] numerically studied wavelet multiresolution method for distributed parameter estimation, and a simple elliptic problem was taken as a suitable model. Several authors have also used wavelet to analyze linear ill-posed problems [16,17]. Recently, wavelet has been applied to solving nonlinear ill-posed problems. For instance, Ghanem et al. [18] presented a wavelet-based approach for identifying the mechanical parameter of zero-memory non-linear discrete structural systems. The presented approach allows both the parameter estimation of a prior known dynamical models and the identification of classes of suitable non-linear models based on input–output data. Chiao et al. [19] presented a wavelet multiresolution procedure for identifying the parameter of non-linear geophysics. Fu et al. [20] consider the problem of estimating the velocity in a two-dimensional acoustic wave equation, and designed a widely convergent method
In this section we summarize the main properties of wavelet as far as they are needed in the following section. For a more detailed introduction to the theory of wavelet we refer to [6,12,23]. 3.1. Wavelet bases Wavelet bases have been documented in numerous textbooks and survey papers (see [11] for a general treatment). With a little effort they can be adapted to fairly general domain O Rd (see [24] for a discussion of the characterizations of function spaces on O by wavelet coefficients). A wavelet basis consists of two type of functions: scaling function fl and wavelet function cl . The index concatenates the usual scale and space parameters j and k. Thus for standard wavelet bases on R, it is defined in Definition 1. Definition 1 (Orthogonal wavelet). A function c A L2 ðRÞ is called an orthogonal wavelet if the set fcj,k ðxÞ ¼ 2j=2 cð2j xkÞg
ð3:1Þ 2
forms an orthonormal basis for L ðRÞ. A wavelet basis on R is labeled by two indices, the shift-index k and the scale-index j. If we collect all basis functions on the same scale in the subspace W j ¼ spanfcj,k 9k A Zg, then we obtain the orthogonal multi-scale decomposition L2 ðRÞ ¼ W j : jAZ
This type of multi-resolution analysis was first introduced by Mallat (see, e.g., [25]). The subspaces Wj are called the scale of the decomposition, and the projection X /f , cj,k SL2 cj,k Q jf ¼ is said to represent the details of f on the scale j. Definition 2 (Multi-resolution analysis). A multi-resolution analysis is a sequence of closed subspaces Vj in L2 ðRÞ, called scaling spaces, satisfying (1) (2) (3) (4) (5) (6)
V j V j þ 1 , for all j A Z, S V is dense in L2 ðRÞ, Tj A Z j j A Z V j ¼ f0g, f ðxÞ A V j if and only if f ð2xÞ A V j þ 1 , 8j A Z, f ðxÞ A V 0 if and only if f ðxkÞ A V 0 , 8k A Z, there exists f A V 0 such that ff0,k : kA Zg
is an orthonormal basis in V0. The function is called the scaling function of the multi-resolution analysis, simply called ‘‘MRA’’. Obviously, the subspaces Vj are spanned by V j ¼ spanffj,k ðxÞ ¼ 2j=2 fð2j xkÞ9k A Zg:
ð3:2Þ
140
H. Fu et al. / Neurocomputing 104 (2013) 138–145
Let us denote Vj ¼ Wl lrj
the subspace of L2 ðRÞ which contains all details of size j or coarser one. However, the notation cl takes into account the possible adaptations of wavelet to multivariate bounded domains in which case the functions cl and fl usually change form near the boundary. At the scale level j which corresponds to resolution 2j , the scaling functions ðfl Þl A Lj span a space Vj within a hierarchy V 0 V 1 L2 ðOÞ of nested approximation spaces. The space Vj should be thought of a finite element space Vh with uniform mesh size h 2j , and local basis ðfl Þl A Lj . The wavelets ðcl Þl A Lj span a complement Wj of Vj into V j þ 1 . With these notations, the wavelet decomposition takes the form X XX f¼ cl fl þ dl cl , ð3:3Þ
on some fine scale J. Then employing the scaling Eqs. (3.4) and (3.5), we obtain recursively X cjk ¼ hn2k cjkþ 1 , ð3:8Þ mAZ j
dk ¼
X
g n2k cjkþ 1 :
ð3:9Þ
mAZ
It remains to compute the values of cJk on a fine scale J. In a realistic situation we are only given discrete values f k , k A Z, of f than the function itself. Not all scaling functions have the interpolating property like the Haar function. But they necessarily R have a non-vanishing mean value [16], fðxÞ dx a 0, hence J J 2 fð2 xÞ approximates the d-distribution, and we simply put f k ¼ cJk on an appropriate scale J.
j Z l l A ,j
l A Ll
where ðfl Þl A Ll is the scaling function basis spanning the approximation at level l, ðcl Þl A Lj is the wavelet basis spanning the details at level j, and cl ¼ cl ðf Þ and dl ¼ dl ðf Þ are wavelet coefficients of f, respectively.
4. Inversion strategies 4.1. Model and wavelet multiscale inversion strategy Consider the following partial differential equation defined in the time–space domain:
3.2. Multiscale decomposition and reconstruction
LðcðxÞ,tÞuðx,tÞ ¼ 0,
For simplicity, we illustrate multiscale decomposition and reconstruction only for one dimension. The case of multidimension is analogous. Let fV j ðxÞ; j A Zg be a one-dimensional multiresolution analysis (MRA), its scale function is fðxÞ, the wavelet function is cðxÞ. Since f A V 0 V 1 and c A W 0 V 1 , it follows that both the wavelet c and scaling function f satisfy the scaling equations pffiffiffi X fðxÞ ¼ 2 hk fð2xkÞ, ð3:4Þ kAZ
pffiffiffi X
cðxÞ ¼ 2
g k fð2xkÞ,
ð3:5Þ
kAZ
T
fðxÞ ¼ fð2xÞ þ fð2x1Þ, cðxÞ ¼ fð2xÞfð2x1Þ:
Euðx,0Þ ¼ gðxÞ, x A O:
Buðx,tÞ ¼ f ðx,tÞ,
x A @O, 0 r t rT,
x A G @O, 0 rt r T,
cj ðxÞ ¼ P j ðcðxÞÞ:
j
ð3:6Þ
is particularly easy and fast. Assume that we know the values of the scalar products cJk ¼ /f , fJ,k SL2
ð3:7Þ
ð4:4Þ
with A the auxiliary condition operator and G a part of @O, (4.1)–(4.4) constitutes the general parameter estimation problem of partial differential equations. Let c A L2 ðRÞ is a wavelet belonging to the Daubechies family of compactly supported wavelets, f A L2 ðRÞ is the corresponding scale function. For j A Z, let
g k ¼ ð1Þk h2N1k ,
dk ¼ /f , cj,k SL2 ,
ð4:3Þ
where E, B are the initial condition operator and boundary condition operator, respectively. If c(x), g(x), f ðx,tÞ are known, (4.1)–(4.3) constitutes a forward problem to determine uðx,tÞ. If c(x) is unknown, given the auxiliary condition
V j ¼ spanffj,i ¼ 2j=2 fð2j xiÞ9i A Zg,
For actual computations one is only interested in wavelets and scaling functions which obey such finite scaling equations. In this case the computation of the wavelet decomposition, i.e., the computation of the scalar products
ð4:2Þ
The boundary condition is
The Haar wavelet is the only obvious compactly supported orthogonal wavelet. It was the accomplishment of Daubechies [6] to construct a family ffN g of orthogonal wavelets with suppðfN Þ ¼ ½0,2N1 and linearly increasing regularity. Here c1 equals the Haar wavelet. The corresponding filter coefficients gk in (3.4) are given by k ¼ 0,1, . . . ,2N1:
ð4:1Þ
p
where x ¼ ðx1 ,x2 , . . . ,xp Þ , O R , @O is the boundary of O, uðx,tÞ is a sufficiently smooth function defined in O ð0,TÞ. L is a differential operator. The initial condition is
Auðx,tÞ ¼ jðx,tÞ,
with real scaling coefficients hk, gk. For example, the Haar function fðxÞ ¼ w½0,1 ðxÞ and the Haar wavelet cðxÞ ¼ w½1,1=2 ðxÞw½1=2,1 ðxÞ satisfy the above conditions with
x A O, 0 r t r T,
W j ¼ spanfcj,i ¼ 2j=2 fð2j xiÞ9i A Zg: Suppose Pj is the projector from L2 ðRÞ to Vj, Pj : L2 ðRÞ-V j , and denote the projection of cðxÞ A L2 ðRÞ by ð4:5Þ
j
Obviously c ðxÞ A V j , and lim cj ðxÞ ¼ cðxÞ:
j- þ 1
ð4:6Þ
Since, in practical computations, the solving domain is often bounded, and the differential equation, the initial condition, the boundary condition, and the auxiliary condition have to be discretized, the original inverse problem will be approximated by a finite dimensional model, and c(x) only has finite resolution. Therefore, in the inversion process of wavelet multiscale method, the value of scale j is taken in a range 0 rj rJ. The approximation
H. Fu et al. / Neurocomputing 104 (2013) 138–145
where dcjl ðxÞ is a small perturbation, which can be determined by the minimization of the following nonlinear functional:
of the function c(x) can now be written as X
cJ ðxÞ ¼
J X X
c0i f0,i ðxÞ þ
i A L0
j
di cj,i ðxÞ,
ð4:7Þ
j ¼ 0 i A Lj
where Lj is an index set. Therefore, the problem to determine cJ ðxÞ is reduced to solving a finite dimensional wavelet coefficients vector. And according to the principle of multiscale inversion, the undetermined parameter c(x) will be decomposed onto different subspaces Vj, such that the problem of determining c(x) is transformed into the problem of determining the projection cj ðxÞ in the subspaces Vj, j ¼ 0,1, . . . ,J. We start from the coarsest level of resolution and progressively move to the finest level. Let us rewrite approximation (4.7) for any scale j ð0 r j r JÞ as X j cj ðxÞ ¼ cj1 ðxÞ þ di cj,i ðxÞ, 0
c ðxÞ ¼
i A Lj
X
c0i f0,i ðxÞ:
ð4:8Þ
Ja ½dcj ðxÞ ¼ JA½cjl ðxÞ þ dcjl jðx,tÞJ2G½0,T þ aJdcjl ðxÞJ2 :
At each scale, cj ðxÞ will be expressed as an expansion of the wavelet basis, and the problem of determining cj ðxÞ is reduced to determine the finite coefficients of representation. The whole inversion will start at the longest scale. At the longest scale j¼0, use the iterative regularization method to make correction on the initial guess until the global minimum is obtained, and the approximate solution will be used as the initial guess for the inversion at scale j ¼1. Successively fining upward in this way until we finally find a good approximate solution cJn ðxÞ at the finest scale j¼J. cJn ðxÞ can be regarded as the optimal solution of the original problem.
to denote the solution of the forward We use problem (4.1)–(4.3) when cjl ðxÞ is given. The corresponding solution of (4.1)–(4.3) is denoted by uðcjl ðxÞ þ dcjl ðxÞ; x,tÞ after being added a minor perturbation X
dcjl ðxÞ ¼
j X X
T ddj,l i cj,i ðxÞ9dK l FðxÞ
ð4:14Þ
j ¼ 0 i A Lj
to cjl ðxÞ. In this way, the problem to determine dcjl ðxÞ is reduced to solve dK l . Whereas dK l can be determined by the minimization of the following functional: F½dK l ¼ JAuðcjl ðxÞ þ dcjl ðxÞ; x,tÞjðx,tÞJ2G½0,T þ aJdK l J2 :
ð4:15Þ
is a minor perturbation, from the multivariate Since Taylor expansion formula, we have T
uðcjl ðxÞ þ dcjl ðxÞ; x,tÞ ¼ uðcjl ðxÞ; x,tÞ þ XK l uðcjl ðxÞ; x,tÞ dK l þ OðJdcjl ðxÞJ2 Þ,
j X X
c0i f0,i ðxÞ þ
i A L0
T
F½dK l ¼ JAuðcjl ðxÞ; x,tÞjðx,tÞ þ XK l Auðcjl ðxÞ; x,tÞ dK l J2G½0,T þ aJdK l J2 :
If taking 2j discrete points ðxm ,t m Þ ðm ¼ 0,1, . . . ,2J 1Þ, then from numerical integration theory, we have
j
A½c ðxÞ ¼ A2 ½A1 ½cð ðxÞ ¼ jðx,tÞ:
ð4:10Þ
j
2
J a ½c ðxÞ ¼ JA½c ðxÞjðx,tÞJG½0,T þ aJc ðxÞJ ,
ð4:18Þ
J
Generally, it is ill-posed in the sense of Hadamard and should be regularized. Tikhonov regularization [26] of such ill-posed problem is achieved by replacing the operator equation (4.8) by the minimization problem ð4:11Þ
where a is the regularization parameter. We will use the perturbation method and linearization technique to establish an iterative method. Let cjl þ 1 ðxÞ ¼ cjl ðxÞ þ dcjl ðxÞ,
rm ½Auðcjl ðxÞ; xm ,tm Þjðxm ,tm Þ
þ XK l Auðcjl ðxÞ; xm ,t m Þ dK l 2 þ aJdK l J2 ,
ð4:9Þ 0
2
J 2X 1
T
j
di cj,i ðxÞ9K T FðxÞ
in the subspace V j ð0 r j rJÞ, where K T ¼ fðc0i Þi A L0 ; ðdi Þi A L0 ; 1 j ðdi Þi A L1 ; ; ðdi Þi A Lj g is the wavelet coefficient in the approximation subspace Vj, FðxÞ ¼ fðf0,i ðxÞÞi A L0 , ðc0,i ðxÞÞi A L0 ,ðc1,i ðxÞÞi A L1 , . . . , ðcj,i ðxÞÞi A Lj gT is the wavelet basis function on the scale j. Therefore, the problem to determine cj ðxÞ is reduced to solving a finite dimensional vector K. On the other hand, since the solution uðx,tÞ of (4.1)–(4.3) is nonlinearly dependent on cj ðxÞ, if we define a nonlinear operator A1 which maps cj ðxÞ to uðcj ðxÞ; x,tÞ, then A1 ½cj ðxÞ ¼ uðcj ðxÞ; x,tÞ. From (4.4) we know that there exists a nonlinear operator A2 mapping from the solution uðcj ðxÞ; x,tÞ to the data jðx,tÞ such that A2 ½uðcj ðxÞ; x,tÞ ¼ jðx,tÞ. Therefore, the inverse problem is equivalent to solve the following nonlinear operator equation:
j
ð4:17Þ
m¼0
j ¼ 0 i A Lj
j
ð4:16Þ
consequently
F½dK l ¼
Approximation (4.7) can be expanded as
j
dc0,l f0,i ðxÞ þ i
i A L0
4.2. An iterative regularization method
cj ðxÞ ¼
ð4:13Þ
uðcjl ðxÞ; x,tÞ
dcjl ðxÞ
i A L0
X
141
ð4:12Þ
where frm g2m 1 ¼ 1 are some adequate weights. Let 2 3 r0 Auðcjl ðxÞ; x0 ,t 0 Þ 6 7 6 7 r1 Auðcjl ðxÞ; x1 ,t1 Þ 6 7 U¼6 7, 6 7 ^ 4 5 r2J 1 Auðcjl ðxÞ; x2J 1 ,t2J 1 Þ 2 3 r0 jðx0 ,t0 Þ 6 7 r1 jðx1 ,t1 Þ 6 7 7, Un ¼ 6 6 7 ^ 4 5 r2J 1 jðx2J 1 ,t2J 1 Þ H ¼ ðhmi Þ2J dimðLj Þ , hmi ¼
@ r Auðcjl ðxÞ; xm ,tm Þ: @ki m
Then (4.16) can be simplified to F½dK l ¼ dK Tl HT HdK l þ 2dK Tl HT ðUU n Þ þ ðUU n ÞT ðUU n Þ þ aJdK l J2 :
ð4:19Þ
Obviously the minimum point of F½dK l satisfies ðHT H þ aIÞdK l ¼ HT ðUU n Þ:
ð4:20Þ
It is just Euler’s equation of (4.17). Solving dK l from (4.18) and substituting it into (4.12), we can finally obtain the perturbation dcjl ðxÞ, and cjl þ 1 ðxÞ further. Repeat the above process until a satisfied precision is obtained.
142
H. Fu et al. / Neurocomputing 104 (2013) 138–145
if f A H1 ðOÞ for almost all x A O, we denote
4.3. The process of the wavelet multiscale iterative regularization method
X ¼ fcðxÞ A L1 ðOÞ : 0 o g rcðxÞ r g , suppðcðxÞc0 ðxÞÞ Og
In summary, we present the main process of the wavelet multiscale iterative regularization method for the parameter estimation problem of differential equation:
and define F : X -Y ¼ H10 ðOÞ,
In each case the forward operator was discretized by using finite elements on a uniform grid (triangular, in the case of two dimensions).
(1) We first make the multiscale decomposition of c(x) by using Daubechies compact supported orthogonal wavelet to obtain cj ðxÞ ðj ¼ J,J1, . . . ,0Þ. The size of J is determined by the resolution ratio of practical data. Given an initial estimate c0 ðxÞ of the true parameter cn ðxÞ. Thus the initial guess c0 ðxÞ P can be approximated as c0 ðxÞ ¼ i A L0 c0i f0,i ðxÞ. We solve the 0 numerical solution uðcl ; xm ,t m Þ of the problem (4.1) and (4.3) by the finite-difference method, and obtain the value Auðc0l ; xm ,t m Þ ðm ¼ 0,1, . . . ,2J 1Þ in the auxiliary condition. (2) Compute the elements hmi of the matrix H by numerical differentiation
5.1. One-dimensional examples In our first example, we take 8 n c ðxÞ ¼ 10099:9ð12xÞ2 , > > > > > uðxÞ ¼ xð1xÞ, > > > > < f ðxÞ ¼ 200599:4ð12xÞ2 , > g ¼ 0:1, g ¼ 200, > > > > > c > 0 ðxÞ ¼ 100, > > : O ¼ ½0,1:
@ Auðc0l ; xm ,t m Þ @ki 1 ¼ ½Auðc0l ðxÞ þ tf0,i ðxÞ; xm ,t m ÞAuðc0l ðxÞ; xm ,t m Þ:
hmi ¼
(3) By the iterative method given in the above subsection, we can obtain the global minimum point c0n at scale 0. (4) Conduct the process (1)–(3) at scale 1 by taking c0n ðxÞ as an initial value for c1n ðxÞ. (5) Repeat the process (1)–(4), until we finally obtain the global minimum point cJn ðxÞ. 5. Numerical examples
Remark. Because the parameter c(x) that we are looking for is required only to be in O ¼ ½0,1, we shall use the simplest wavelet and scaling function, namely the Haar wavelet and its scaling function.
To show the feasibility of our approach and illustrate its numerical behavior, we consider the elliptical differential equation
r ðcðxÞruÞ ¼ f
in O Rd ,
u¼0
on @O,
The second model is depicted in Fig. 2. In the experiments, we try to identify a piecewise continuous c(x) with several jump discontinuities by using the proposed method. We take (1 x A ½0, 15 [ ½25 , 35 [ ½45 ,1, 2, cðxÞ ¼ ð5:3Þ 1, x A ½15 , 25 [ ½35 , 45:
ð5:1Þ
The wavelet basis and parameters a, h are the same as the first example.
110
110
100
100
90
90
80
80
70
70
c (x)
c (x)
where f A L2 ðOÞ. The observations of the solution u are used to recover the coefficient c(x). Our task is to determine an approximation of parameter c(x) from a measurement ud in the model problem (5.1) for d¼1,2. Since the model problem is only well-posed for cðxÞ Z g 40, and
60 50
60 50
exact solution itertion solution
40
exact solution itertion solution
40
30 20
ð5:2Þ
The solution of the forward problem (5.1) was based on linear finite elements on a uniform grid with mesh size h ¼ 1=26 . In this case, we use wavelet multiscale decomposition method to decompose the parameter c(x) into six scales, thus we obtain a sequence of inverse problems which depend on the scale variables j ¼ 0,1, . . . ,6. At each scale j, we take regularization parameter a ¼ 2j 104 , the iteration will terminate when the iteration number kj r20 or the relative error is equal or less than 1=ð10ðj þ 1ÞÞ 100%. To test the noisy suppression ability of our methods, we invert the parameter c(x) when there is 10% random noise added to the data u. Fig. 1 illustrates the numerical results.
t
(
cðxÞ/u ðsolution of ð5:1ÞÞ:
30 0
0.2
0.4
x
0.6
0.8
1
20
0
0.2
0.4
x
0.6
0.8
1
Fig. 1. Inversion results for the first example in one-dimensional. (a) There is no noise. (b) There is 10% noise added to data.
H. Fu et al. / Neurocomputing 104 (2013) 138–145
3
3 exact solution itertion solution
2.5
2 c (x)→
c (x)→
exact solution itertion solution
2.5
2 1.5
1.5
1
1
0.5
0.5
0
143
0
0.2
0.4
x→
0.6
0.8
1
0
0
0.2
0.4
0.6 x→
0.8
1
Fig. 2. Inversion results for the second example in one-dimensional. (a) There is no noise. (b) There is 10% noise added to data.
Table 1 The numerical results of the wavelet multiscale iterative regularization method with noisy data (J ¼ 6, C ¼ 1). k
d%
JcJl cJn J JcJn J
d%
k
JcJl cJn J JcJn J
1:00 101
58
1:023 100
1:00 100
60
0:834 100
0:50 101
49
0:815 100
0:50 100
55
0:658 100
1:00 100
42
0:364 100
1:00 100
46
0:243 100
0:50 100
39
0:522 101
0:50 101
41
0:276 101
1:00 101
36
0:204 101
1:00 101
39
0:185 101
1
0:50 10
33
0:418 10
2
1
0:50 10
37
0:214 102
0:25 101
32
0:221 102
0:25 101
32
0:102 102
Jcjl cjn J=Jcjn J
In the process of computation, we take ðj ¼ 0, 1,2, . . . ,6Þ as relative error levels. Fig. 2 illustrates the inversion results. Experiment results with noisy data are displayed in Table 1. Here we use uniformly distributed random noise. For each noise level, we carry out four test runs and compute the average relative error JcJl cJn J=JcJn JðJ ¼ 6Þ. The left three columns in Table 1 show the results for the first example with the initial value c0 ðxÞ 100. In the right three columns we display the results for the second example with the initial value c0 ðxÞ 1. Here, k denotes the iteration total number by using the wavelet multiscale iterative regularization method.
5.2. Two-dimensional examples To illustrate the applicability of the proposed method to higher-dimensional parameter identification problems, we also consider the model problem in two-dimension. Here, we take O ¼ ½0,12 , f ðx1 ,x2 Þ 1, and " # " # ðx1 0:5Þ2 ðx2 0:5Þ2 cðx1 ,x2 Þ ¼ exp exp , 0:05 0:05 see Fig. 3(a). Due to the specific geometry of O, we can use periodized Daubechies wavelets of period 1. In particular, the space Vj ¼ V j V j Vj þ 1 ,
jAZ
form a separable two-dimensional MRA with scale function
Fðx1 ,x2 Þ ¼ fðx1 Þfðx2 Þ. The orthogonal complement of Vj in Vj þ 1 is Wj ¼ Whj Wvj Wdj ,
where W0 is defined by the horizontal, vertical and diagonal wavelet functions
Ch ðx1 ,x2 Þ ¼ fðx1 Þcðx2 Þ, Cv ðx1 ,x2 Þ ¼ cðx1 Þfðx2 Þ, Cd ðx1 ,x2 Þ ¼ cðx1 Þcðx2 Þ: The [0,1]-periodic scaling function and wavelets are defined by the formulae XX j Fk1 ,k2 ðx1 l1 ,x2 l2 Þ l1 A Z l2 A Z
and XX
Cj,k1l,k2 ðx1 l1 ,x2 l2 Þ, l ¼ h,v,d,
l1 A Z l2 A Z
respectively. To obtain un and ud , we solve (5.1) by the finite element method on a mesh with N ¼ 64 64 elements. In numerical simulations, we use two-dimensional fast wavelet transform to decompose the parameter c(x) into six scales. We first consider the case when there is no noise, i.e. ud ¼ un . Fig. 3(b) and (c) shows un and the identified coefficients cJn ðx1 ,x2 ÞÞ from the noise-free observation. It shows that our algorithms have reconstructed the parameter c(x) very accurately. Second, we consider noisy data. The noise level is 12.14%. The reconstruction results, shown in Fig. 3(d), approximate well to the true solution. Numerical results indicate (1) when there is no noise in the data, for the two different examples in one dimension, the inversion results have high precisions. Besides, at scale 0, only one coefficient needs to be inverted, and only two iterations are performed to obtain an approximate value with 10% relative error. Although iterative number will increase with the increment of scale, since the final approximation at a former scale will be used as a good initial value for the latter, even at the finest scale—scale 6, only 13 and 15 iterations are needed to obtain the global optimal solution for this two examples, respectively. It indicates the advantages of the wavelet multiscale inversion method on the aspects of computational precision, efficiency and convergent range. (2) When the noise is added, the number of iteration will increase, and the precision of inversion results are lowered a bit. But because wavelet multiscale decomposition has a good ability of noise suppression, the effectiveness of noise will gradually drop. Therefore, we can still obtain almost perfect reconstruction results even 12.14% noise is added in two dimension.
144
H. Fu et al. / Neurocomputing 104 (2013) 138–145
Fig. 3. Inversion results for the 2D example. (a) Exact solution. (b) Exact data u. (c) There is no noise. (d) There is noise added to data.
6. Conclusions For the general parameter estimation problem of partial differential equations, a wavelet multiscale iterative regularization method is proposed in this paper. Numerical simulations for the inversion problem of one- and two-dimensional elliptical equation showed the method’s effectiveness in the aspects of stability, global convergence, noise suppression and fast computation. Wavelet multiscale method widens the scope of research of multiscale inversion. To conduct wavelet multiscale decomposition for the undetermined parameter will lower the dependence of the inverse problems on initial estimate, and overcome the effect of local minimum in the objective functional. The computational cost is greatly reduced by transforming the undetermined parameter into the finite wavelet coefficients. Moreover, because there are ready-made fast algorithms of decomposition and reconstruction, the operation is simple, fast, and stable. It is helpful to solve more complicated inverse problems. However, many theoretical and computational problems are still open.
Acknowledgments This work is partly supported by the National Natural Science Foundation of China (Grant nos. 41074088 and 61173035), the Doctoral Science Technology Foundation of Liaoning Province (Grant no. 20091007), the China Postdoctoral Science Foundation (Grant no. 20110491533), the Fundamental Research Funds for the Central Universities (Grant no. 2012TD027), the Program for New Century Excellent Talents in University (Grant no. NCET-110861).
References [1] V. Isakov, Inverse Problems for Partial Differential Equations, vol. 127, Springer-Verlag, 2006.
[2] B. Choi, J. Lee, Comparison of generalization ability on solving differential equations using backpropagation and reformulated radial basis function networks, Neurocomputing 73 (1–3) (2009) 115–118. [3] D. Luengo, S. Monzon, A. Artes-Rodriguez, Novel fast random search clustering algorithm for mixing matrix identification in mimo linear blind inverse problems with sparse inputs, Neurocomputing 87 (2012) 62–78. [4] D. Benaouda, F. Murtagh, J. Starck, O. Renaud, Wavelet-based nonlinear multiscale decomposition model for electricity load forecasting, Neurocomputing 70 (1–3) (2006) 139–154. [5] M. Mu, Q. Ruan, S. Guo, Shift and gray scale invariant features for palmprint identification using complex directional wavelet and local binary pattern, Neurocomputing 74 (17) (2011) 3351–3360. [6] I. Daubechies, Ten Lectures on Wavelet, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1992. [7] Y. Meyer, D. Salinger, C.U. Press, Wavelets And Operators, vol. 2, Cambridge University Press, 1992. [8] G. Beylkin, On the representation of operators in bases of compactly supported wavelets, SIAM J. Numer. Anal. 6 (6) (1992) 1716–1740. [9] W. Dahmen, Wavelet methods for PDEs—some recent developments, J. Comput. Appl. Math. 128 (1) (2001) 133–185. ¨ [10] J. Frohlich, K. Schneider, An adaptive wavelet—vaguelette algorithm for the solution of pdes, J. Comput. Phys. 130 (2) (1997) 174–190. [11] V. Comincioli, G. Naldi, T. Scapolla, A wavelet-based method for numerical solution of nonlinear evolution equations, Appl. Numer. Math. 33 (1–4) (2000) 291–297. [12] V. Dicken, P. Maaß, Wavelet-Galerkin methods for ill-posed problems, J. Inver. Ill-Posed Prob. 4 (3) (1996) 203–222. [13] S. Dahlke, P. Maaß, An outline of adaptive wavelet Galerkin methods for Tikhonov regularization of inverse parabolic problems, in: Recent Development in Theories and Numerics. Proceedings of the International Conference on Inverse Problems, 2003, pp. 56–66. [14] A. Cohen, M. Hoffmann, M. Reiss, Adaptive wavelet Galerkin methods for linear inverse problems, SIAM J. Numer. Anal. 42 (4) (2005) 1479–1501. [15] J. Liu, A multiresolution method for distributed parameter estimations, SIAM J. Sci. Comput. 14 (2) (1991) 389–405. [16] D. Donoho, Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition, Appl. Comput. Harmonic Anal. 2 (2) (1995) 101–126. [17] T. Reginska, L. Elde´n, Stability and convergence of a wavelet-Galerkin method for the sideways heat equation, J. Inver. Ill-Posed Prob. 8 (3) (2000) 31–49. [18] R. Ghanem, F. Romeo, A wavelet-based approach for model and parameter identification of non-linear systems, Int. J. Non-lin. Mech. 36 (5) (2001) 835–859. [19] L. Chiao, W. Liang, Multiresolution parameterization for geophysical inverse problems, Geophysics 68 (1) (2003) 199–209. [20] H. Fu, B. Han, A wavelet multiscale method for the inverse problems of a twodimensional wave equation, Inver. Prob. Sci. Eng. 12 (6) (2004) 643–656. [21] X. Lu, Y. Sun, G. Bai, Adaptive wavelet-Galerkin methods for limited angle tomography, Image Vis. Comput. 28 (4) (2010) 696–703.
H. Fu et al. / Neurocomputing 104 (2013) 138–145
[22] A. Qian, A new wavelet method for solving an ill-posed problem, Appl. Math. Comput. 203 (2) (2008) 635–640. [23] A. Louis, P. Maaß, A. Rieder, Wavelets: Eine Einfuhrung in Theorie und Anwendungen, Teubner Verlag, Stuttgart, 1994. [24] A. Cohen, Numerical analysis of wavelet methods, Stud. Math. Appl. 32 (2003) 1–335. [25] S. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2 ðRÞ, Trans. Am. Math. Soc. 315 (1) (1989) 69–87. [26] A. Tikhonov, V. Arsenin, F. John, Solutions of Ill-Posed Problems, Winston Washington, DC, 1977.
Hongsun Fu received MS and PhD at Harbin Institute of Technology, China. She is an Associate Professor of the Department of Mathematics, Dalian Maritime University. She serves as a referee for the inverse problems in Science & Engineering.
Bo Han received MS and PhD at the Harbin Institute of Technology, China. He is a Professor of the Department of Mathematics, Harbin Institute of Technology. He is an Executive Director of Society for Industrial and Applied Mathematics in China. He is also the President of Society for Industrial and Applied Mathematics in Heilongjiang Province, China. He serves as a Referee for the Mathematical Reviews of Chinese Mathematical Society. Also he serves as a Referee for International Journals: Multiscale Modeling and Simulation (SIAM); International Journal on Numerical Analysis and Modeling; Selected Publications of Chinese Universities: Mathematics; Journal on Information and Computational Science; Communications in Computational Physics; Journal of Computational Mathematics.
145 Hongbo Liu received his three level educations (BSc, MSc, PhD) at the Dalian University of Technology, China. Currently he is a Professor at School of Information Engineering and Science, Dalian Maritime University, with an affiliate appointment in the Computer Science and Biomedical Engineering Division at the Dalian University of Technology. His research interests are in system modeling and optimization involving soft computing, probabilistic modeling, cognitive computing, machine learning, data mining, etc. He participates and organizes actively International Conference and Workshop and International Journals/Publications.
