Aristotelian prior boundary conditions 1 Introduction

19 Pages·2012·1.08 MB·English

(cid:0)(cid:2)(cid:1)
International Journal of Mathematics M
CS
and Computer Science, 1(2006), 63-81
Aristotelian prior boundary conditions
Daniela Calvetti1, Jari P. Kaipio2, Erkki Someralo3
1 Case Western Reserve University, Department of Mathematics
10900 Euclid Ave., Cleveland, OH 44106, USA.
e-mail: [email protected]
2 University of Kuopio, Department of Applied Physics
PO Box 1627, FIN-70211 Kuopio, Finland.
e-mail: jari.kaipio@uku.fl
3 Helsinki University of Technology, Institute of Mathematics
PO Box 1100, FIN-02015 HUT, Finland.
e-mail: erkki.somersalo@hut.fl
(Received June 23, 2005, Accepted October 21, 2005)
Abstract
The selection of boundary conditions when restoring a flnite 1D or 2D signal
has received a lot of interest in recent times. Dirichlet, periodic, re(cid:176)ective
and antire(cid:176)ective boundary conditions may be appropriate selection for some
problems, while wholly unsuitable for others. In this paper we show that, in
an Aristotelian approach to knowledge, when it is not known a priori which
boundary conditions should be chosen, by admitting our lack of information
it is possible to let the data itself determine them. We consider 1D and 2D
smoothness prior conditions. The application of this Aristotelian approach to
boundary conditions to the solution of linear discrete ill-posed problems with
Tikhonovregularizationortruncatediterativemethodsisdiscussed. Computed
examples of deblurring and limited angle tomography showing the efiectiveness
of Aristotelian boundary conditions are presented.
1 Introduction
The selection of suitable boundary conditions when restoring a flnite signal
from a distorted specimen may play an important role in the success of the
outcome. The issue of which boundary conditions are most appropriate has
been addressed repeatedly in the image deblurring literature, where Dirichlet,
Neumann, peridic, re(cid:176)ective and antire(cid:176)ective boundary conditions have been
proposed; see, e.g., [3, 6] and references therein. It is clear from the examples
reported in the literature that it is always possible to flnd problems for which
Research of the flrst author was supported in part by grant NIH GM-66309-01.
ISSN 1814-0424 (cid:176)c 2006, http://ijmcs.futureintech.net
Aristotelian prior boundary conditions 64
one type of boundary conditions work much better than the others. For ex-
ample, if the signal that we want to recover does not vanish on the boundary,
imposing Dirichlet boundary conditions may cause the computed solution to
show oscillatory behavior near the boundary. Similar misbehavior of the com-
puted solution can be observed with the other boundary conditions when used
inappropriately.
In this paper we consider the case in which we have no deflnite a priori knowl-
edge as to which type of boundary conditions might be most appropriate for
our problem. Instead of guessing which boundary conditions might be most
appropriate, or selecting the boundary conditions according to how much such
selectionwillsimplifythecomputations,weacknowledgeour\tabularasa"state
and let our experience of the signal, i.e., the data itself, determine what is the
most likely behavior of the solution at the boundary. We call the boundary
conditions determined in this fashion Aristotelian boundary conditions because
we, like the Greek philosopher, assume that knowledge is a process guided by
experiences,builtlayerafterlayerstartingfromacleanslate. 1 Inthispaperwe
propose an approach to boundary condition selection that allows a great deal
of (cid:176)exibility in terms of expressing the various degrees of knowledge, or lack
thereof, about the behavior of the solution at the boundary, and an algorithm
which is computationally e–cient also for problems of large dimensions.
We remark that while the basic ideas behind Aristotelian boundary conditions
areessentiallythesameinone,twoorthreedimensionalspace,thetechnicalde-
tailsbecomemoretediousinhigherdimensions. Thereforewewillflrstdescribe
the algorithm in one dimension, then explain how it should be modifled for ap-
plicationtohigherdimensionalproblems. Furthermore,sincealotoftherecent
discussions on the choice of boundary conditions has focused on the problem of
signal and image deblurring, we will also assume flrst that we want to recover
a signal x from a blurred, and possibly noisy copy b. After discretization we
have the linear system
b=Ax+e; x Rn; b; e Rm; (1)
2 2
where the matrix A Rm£n is of ill-determined rank, that is, it has singular
2
values of difierent orders of magnitude close to zero. In order to keep the
ampliflederrorcomponentsfromdominatingthecomputedsolution, someform
of regularization is typically needed when solving (1). Here we consider only
two regularization methods. The flrst one, Tikhonov regularization, replaces
the solution of (1) with the solution of the minimization problem
min b Ax 2+fi Lx 2 ; (2)
fk ¡ k k k g
where denotestheEuclideannormandLisaregularizationoperatorwhich
k¢k
deflnes a seminorm. The minimizer exists and is unique if ker(A) ker(L) =
\
1Aristotle’spositionwastheoppositeofthatofhisteacherPlato,whoregardedknowledge
asatraitinnateinthehumanmind.
Aristotelian prior boundary conditions 65
0 . In particular, if L Rn£n is invertible, the problem (2) is well-posed.
f g 2
The second method considered here is truncated iteration. Regularization by
truncated iteration consists of applying a few steps of an iterative method to
the linear system (1), or to its normal equations, stopping the iteration prior
to convergence to avoid the harmful efiects of the amplifled error components
on the computed solution. The use of smoothing preconditioners in connection
with truncated iteration may improve the quality of the computed solution. In
particular, assuming for simplicity that the matrix A is square and that the
matrix L in (2) is invertible, the solution by truncated iteration of the right-
preconditioned linear system
AL¡1w=b; Lx=w; (3)
may be computationally more e–cient than the solution of (1), see [2].
In this paper we will show how to express our a priori knowledge, or lack
thereof, about the sought solution and its behavior at the boundary by means
ofaninvertiblepositivedeflnitematrixLtobeusedasaregularizationoperator
for Tikhonov or as a right preconditioner in truncated iteration. The matrix L
may, for example, express the fact that the solution is smooth along a portion
of the boundary, without forcing any particular type of boundary conditions.
We will call the matrix L Aristotelian boundary preconditioner when used in
connection with iterative methods for linear systems.
Thepaperisorganizedasfollows. InSection2weshowhowstatisticalinversion
provides the tools both for expressing various levels of uncertainty about the
behavior of the solution at the boundary and for having the data itself deter-
mine what this behavior should be. In Section 3 we describe how to construct
the matrix L in the one-dimensional model case. Section 4 describes its ex-
tension to higher dimensions. We also consider smoothness priors related to
structural information, and we present computed examples where Aristotelian
prior boundary conditions are used in the context of limited angle tomography.
2 Inverse problems, regularization and statistics
The construction of the smoothness prior with Aristotelian boundary condi-
tions is based on the Bayesian interpretation of regularization. Now we brie(cid:176)y
recall the basic concepts of the statistical theory of inverse problems. For a
comprehensive discussion of this topic, we refer to [4].
Assume that x Rn represents a quantity we are interested in. Since we are
2
not sure about its value, we model it as a random variable. The randomness
expresses our degree of uncertainty, which in turn is encoded in its probability
distribution. We adopt the convention that a random variable is denoted by
a capital letter, X, and its realizations by lower case letters, x. Further, we
Aristotelian prior boundary conditions 66
denote the probability density2 of X by …(x).
Consider the linear inverse problem (1) of estimating x from the noisy data b.
Assume that prior to the measurement of b, we have some information about
the possible distribution of the unknown x. This information is encoded in
the prior density of the corresponding random variable X, denoted by … (x).
pr
Furthermore,assumethattheerroreisarealizationofanoiseprocessEwhose
probability density is … (e). Then, assuming for the moment that X=x is
noise
flxed, the equation (1) implies that b is a realization of a random variable B
whose probability density must be
…(b x)=… (b Ax): (4)
noise
j ¡
The conditional density (4) is known as the likelihood density. The joint prob-
ability density of the variables B and X is then
…(x;b)=…(b x)… (x)=… (b Ax)… (x):
pr noise pr
j ¡
Bayes formula states that the posterior density of x is proportional to the joint
probability density, that is,
…(x b) … (b Ax)… (x):
noise pr
j / ¡
AparticularbutverycommoncaseiswhenthevariablesXandEareindepen-
dent and Gaussian. Assume, for simplicity, that
E (0;(cid:190)2I); X (0;¡);
»N »N
where I is the m m identity matrix and ¡ Rn£n is a symmetric positive
£ 2
deflnitematrix. ThenBayesformulaimpliesthattheposteriordensityis, upto
a norming constant,
1 1
…(x b) exp b Ax 2+xT¡¡1x : (5)
j / ¡2 (cid:190)2k ¡ k
(cid:181) ‰ (cid:190)¶
Let ¡¡1 = (cid:176)LTL. Such decomposition exists because ¡ is symmetric positive
deflnite and therefore admits a Cholesky factorization. Then the maximum a
posteriori estimate of x, x , is the maximizer of the posterior density (5) or,
MAP
equivalently,
x =argmin b Ax 2+(cid:190)2(cid:176) Lx 2 :
MAP
k ¡ k k k
Hence the regularization matrix L'for Tikhonov regulariza“tion corresponds to
the covariance matrix of the prior probability density. Thus this matrix re(cid:176)ects
our prior information about the unknown. The regularization parameter fi in
(2) is related to the statistical parameters via fi = (cid:190)2(cid:176), giving insight of its
statistical interpretation in terms of the noise level and scaling of the prior
covariance. The determination of fi is usually based on discrepancy principle,
2To avoid measure theoretic considerations, all probability densities are assumed to be
absolutelycontinuouswithrespecttotheLebesguemeasure.
Aristotelian prior boundary conditions 67
while the statistical parameters can be estimated by hierarchical models. For
further discussion of these topics, see [4].
For the purposes of this article, the key observation is that the statistical ap-
proachgivesusmeansofinterpretingtheregularization. Indeed, ifL Rn£n is
2
an invertible matrix appearing inside the Tikhonov functional (2) or as a pre-
conditioner in (3), the matrix B = LTL deflnes, up to a scaling constant, the
inverse of a covariance. Exploration of the corresponding prior density e.g. by
random draws, gives us a way to check if this matrix corresponds to what we
assume to be known a priori about the unknown x.
3 One dimensional smoothness prior
We begin this section by assuming that the desired solution is a smooth one-
dimensional signal f(t) with support on the unit interval [0;1] and we denote
by x=[x ;x ;:::;x ]T its values at the points of a uniformly spaced grid.
1 2 n
Assuming that we believe a priori that f C([0;1]) C2(]0;1[), we seek to
2 \
construct a prior density of x that re(cid:176)ects this information. Observe that no
assumption whatsoever of the boundary conditions of f is made.
Thenaturalstartingpointforasmoothnessprioristoconsideraflnitedifierence
approximation of f00(t) at the interior points t ;:::;t . We have
2 n¡1
f00(t ) 1 2 1 x
2 1
¡
f00(t3) 1 1 2 1 x2 1
2 ... 3… h2 2 ¡... ... ... 32 ... 3= h2Lx; (6)
6 7 6 76 7
6 f00(t ) 7 6 1 2 1 76 x 7
6 n¡1 7 6 ¡ 76 n 7
4 5 4 54 5
whereh=1=(n 1). Thisapproximationgivesrisetoacandidatepriordensity,
¡
which we refer to as smoothness preprior, deflned by
1 1
… (x) exp fi Lx 2 =exp fixTBx ;
pre
/ ¡2 k k ¡2
(cid:181) ¶ (cid:181) ¶
where
B =LTL Rn£n; L R(n¡2)£n:
2 2
We observe that since
rank(B)=n 2;
¡
… is not a proper Gaussian probability density. The rank deflciency of B
pre
alone does not prevent from using L as regularizing matrix, although care must
be taken if the null spaces of L and A coincide. For certain techniques such as
right priorconditioning ([2]), the invertibility of B is essential. Consequently,
we want to augment L in such a manner that the resulting prior is a proper
density.
Aristotelian prior boundary conditions 68
0.04
0.05
0.03 0.04
0.02 0.03
0.01 (s D)j 0.02 (s A)j
0.01
0 0
−0.01
−0.01
−0.02
−0.02 −0.03
−0.03 −0.04
−0.05
−0.04
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Figure 1: Six random draws from the prior … (left) and … (right) with the
D A
parametervaluefi=1. ThedrawsaregeneratedbysolvingthesystemLx=w,
L = L or L = L , where w is a white noise Gaussian random draw. The
D A
standard deviations of the elements are plotted as a bold curve.
AstraightforwardwayofaugmentingLsoastoobtainaninvertiblen nmatrix
£
is to assume boundary conditions at t=0 and t=1, e.g., Dirichlet, Neumann,
Robin or re(cid:176)ecting boundary conditions. Another way is to assume that f
admits a C2-extension outside the interval [0;1] and the extension vanishes in
the extended grid. Such an assumption leads naturally to the extension of L as
2 1
¡
1 2 1 1 2 1
¡ 2 ¡ 3
1 2 1 1 2 1
L=2 ¡... ... ... 3!66 ¡... ... ... 77=LD
6 7 6 7
6 1 2 1 7 6 1 2 1 7
6 ¡ 7 6 ¡ 7
4 5 6 1 2 7
6 ¡ 7
4 5
The corresponding probability density is
1
… (x) exp fi L x 2 ;
D D
/ ¡2 k k
(cid:181) ¶
which is a proper density since L , and consequently LTL , is invertible.
D D D
To understand which prior assumptions are implied by this extension of L, we
make few random draws from the resulting density. The random draws are
facilitated by the observation that if W is Gaussian white noise, i.e., W
»
(0;I),thenX=(pfiL )¡1Wisdistributedaccordingto… . InFigure1,six
D D
N
independent random draws from this density are shown. The boundary values
arenearzero,indicatingthatbyaugmentingLinthisfashion,wehaveimposed
thatthevaluesofthesolutionattheboundaryaresmall. Tobetterunderstand
the relation between the matrix L and the constraint on the solution at the
boundary,considerthevariancesofindividualcomponentsofX … . Wehave
D
»
(cid:190) 2 =E X2 = ¡ ;
D j j D jj
¡ ¢ ' “ ¡ ¢
Aristotelian prior boundary conditions 69
where the covariance matrix ¡ is
D
¡ = fiLTL ¡1:
D D D
The standard deviations (cid:190) are p£lotted in⁄Figure 1. It is clear from the plot
D j
of (cid:190) that the variance at the endpoints is much smaller than in the interior of
D ¡ ¢
the interval. We are now going to relax the constraint at the boundary without
afiecting the smoothness properties inside the interval.
The extension L L was based on some implicit assumptions about the
D
!
boundary values of f. We now analyze the dependence of the interior values on
theboundaryvalues. Wepermutetheentriesofthevectorxsothattheinterior
points come flrst, followed by the boundary points and we use the notation
x n 2
x=[x2;x3;:::;xn¡1;xn;x1]T = x1 ¡2 :
2
• ‚
Accordingly, we partition the matrix L as
L= L1 L2 ; L1 R(n¡2)£(n¡2); L2 R(n¡2)£2:
2 2
Observe that the b£lock L⁄ is invertible. The partition of L induces a partition
1
of B of the form
LTL LTL B B
B =LTL= 1 1 1 2 = 11 12 :
LTL LTL B B
• 2 1 2 2 ‚ • 21 22 ‚
Since the matrix L is invertible, so is B . If the boundary values of f and
1 11
thus the vector x were known, we could calculate the conditional probability
2
density of x conditioned on this information using the preprior. From now on,
1
we assume that fi=1 in the preprior. We have
xTBx = xTB x +xTB x +xTB x +xTB x
1 11 1 1 12 2 2 21 1 2 22 2
= x +B¡1B x TB x +B¡1B x (7)
1 11 12 2 11 1 11 12 2
¡ ¢ +x¡T B B B¡1¢B x :
2 22¡ 21 11 12 2
Notice that since the last term in (7) does not¡depend on x , it on¢ly afiects the
1
norming constant in the exponential deflning the preprior. The quantity
B =B B B¡1B
11 22¡ 21 11 12
is the Schur complement of B . In this setting, the matrix B vanishes. In
b 11 11
fact, from the relationship between the partition of B and that of L, we have
that b
B B B¡1B =LTL LTL LTL ¡1LTL =0:
22¡ 21 11 12 2 2¡ 2 1 1 1 1 2
Consequently, the conditional distribution of X£1 cond⁄itioned on X2 = x2 is a
proper Gaussian density,
… (x x ) B¡1B x ; B¡1 :
pre 1 j 2 »N ¡ 11 12 2 11
¡ ¢
Aristotelian prior boundary conditions 70
We recall that this density corresponds to the assumption that the Dirichlet
boundary values are known, while in reality, they should be treated as ran-
dom variables. Assume that the boundary values are Gaussian, with marginal
distribution
x (0;C);
2
»N
where C R2£2 is a symmetric positive deflnite matrix whose inverse has a
2
factorization of the form
C¡1 =KTK:
The marginal probability density of X is then
2
1 1
… (x ) exp xTC¡1x =exp Kx 2 :
0 2 / ¡2 2 2 ¡2k 2k
(cid:181) ¶ (cid:181) ¶
By using the identity
…(x ;x )=…(x x )…(x ); (8)
1 2 1 2 2
j
we can now deflne a proper prior probability density … as
pr
… (x) = … (x x )… (x )
pr pre 1 2 0 2
j
exp 1 x +B¡1B x TB x +B¡1B x + Kx 2
/ ¡2 1 11 12 2 11 1 11 12 2 k 2k
(cid:181) ¶
£¡ ¢ ¡ ¢ ⁄
1 B B x
= exp xT xT 11 12 1 :
(cid:181)¡2 1 2 • B21 B21B1¡11B12+C¡1 ‚• x2 ‚¶
£ ⁄
Furthermore, the matrix appearing above allows a factorization in terms of L ,
1
L and K of the form
2
B B LTL LTL
11 12 = 1 1 1 2
B B B¡1B +C¡1 LTL LTL +KTK
• 21 21 11 12 ‚ • 2 1 2 2 ‚
T
L L L L
= 1 2 1 2 =LTL ;
0 K 0 K A A
• ‚ • ‚
which, in turn, implies that
1
… (x) exp L x 2 :
pr A
/ ¡2k k
(cid:181) ¶
Note that
det(L )=det(L )det(K)=0:
A 1
6
This is equivalent to assuming that we have an observation model of the form
X =AX +E; E 0;B¡1 ;
1 2 »N 1;1
where ¡ ¢
A= B¡1B
¡ 11 12
Aristotelian prior boundary conditions 71
is the discretization of an operator that continues the boundary values of a
difierential equation inside the domain.
If the boundary values are independent Gaussian random variables with equal
variance, then
1 0 1 1 0
C =(cid:190)2 ; K = :
0 1 (cid:190) 0 1
• ‚ • ‚
SubstitutingtheexpressionforLintothematrixL andrepermutingtheentries
A
of the vector x so that the pixels appear in the natural ordering we have that
1=(cid:190) 0
1 2 1
2 ¡ 3
1 2 1
LA =66 ¡... ... ... 77:
6 7
6 1 2 1 7
6 ¡ 7
6 0 1=(cid:190) 7
6 7
4 5
The selection of the value of the parameter (cid:190) can be done so that the variance
of the entries of the vector x is as uniform as possible. This can be achieved by
letting
(cid:190) = max (cid:190) ;
1•j•n D j
where (cid:190) is computed from the density …¡. I¢n fact,
D D
(cid:190) 2 =var X = LTL ¡1 ;
D j j D D jj
so denoting by e the¡unit¢vector w¡ith¢the¡£jth com⁄po¢nent equal to one,
j
(cid:190) 2 =eT LTL ¡1e = L¡Te 2:
D j j D D j k D jk
¡ ¢ £ ⁄
We choose j = [n=2], since the maximum variance occurs at the center of the
interval.
Figure 1 shows six random draws from the density … , as well as the standard
pr
deviation of the components X , i.e., the square root of the diagonal of the
j
matrix LTL ¡1. Therandomdrawsaregeneratedusingthesamerealizations
A A
of the white noise as in Figure 1.
£ ⁄
We now compare the efiects of the matrices L and L with an example where
D A
L and L are used as regularization operators for Tikhonov regularization.
A D
To show how the matrix L performs as a Tikhonov penalty, we consider the
A
following deblurring problem. Let
1
b = K(s ;t)f(t)dt+e ; 1 j m;
j j j
• •
Z0
where the kernel K is deflned as
K(s;t)=c(s)s(1 s)e¡(t¡s)2=2w2; w =0:03;
¡
Aristotelian prior boundary conditions 72
1
true
Zero extension
0.8
Aristotelian
0.6
0.4
0.2
0
−0.2
−0.4
0 0.2 0.4 0.6 0.8 1
Figure 2: Tikhonov regularized solutions of the deblurring problem using the
regularization matrices L and L , respectively.
D A
and the scaling function c chosen so that max K(s;t) = 1. The sampling
t2[0;1]
points are
s =(j 1)=(m 1); m=40:
j
¡ ¡
The problem is discretized by approximating
1 1 n
K(s ;t)f(t)dt K(s ;t )x ;
j j ‘ ‘
… n
Z0 ‘=1
X
where the points t and s coincide, hence n = m = 40. The noise is assumed
j j
to be Gaussian white noise,
E (0;(cid:176)2I); (cid:176) = 1% of max. of the noiseless signal.
»N
To avoid the \inverse crime" (see [4]) of computing the data using the same
discretization grid in which the inverse problem is solved, we generate the data
using a uniformly spaced grid of 100 discretization points. The inverse problem
is solved using the Tikhonov regularization scheme, i.e., by solving
x =argmin b Ax 2+fi Lx 2 ;
fi
k ¡ k k k
where either L=L or L=L . T¡he parameter fi>0 in¢ both cases is selected
D A
using the Morozov discrepancy principle, i.e., from the condition
b Ax 2 =E E 2 =m(cid:176)2:
fi
k ¡ k k k g
The results are shown in Figure 2. T'he improvements near the boundaries
obtained when using L are clearly visible.
A