- Email: [email protected]
Tikhonov regularisation in standard form for polynomial basis conversion
Tikhonov regularisation in standard form for polynomial basis conversion Joab R. Winkler Department of Computer Science, The University of Shefield, Shefield,
UK
Polynomial basis conversion is an important industrial requirement for high-integn’ty data exchange between computer-aided design systems. These conversions may be numen’cally ill-conditioned, and this paper considers the use of Tikhonov regulatisation in standard form for their stabilisation. It is shown that the method yields a good approximation to the theoretically exact solution if it is known a priori that the desired solution satisfies the discrete Picard condition. An example of the use of Tikhonov regularisation for the stabilisation of the conversion between the Bernstein and B-spline polynomial bases is given. This basis transformation can also be performed by the Oslo algorithm, but it is difficult to determine whether this algotithm ir unstable for given knot vectors. However when this basis transformation is performed by Tikhonov regularisation, it is readily determined whether the linear algebraic equation is unstable, and if it is, regulan’sation follows immediately if the discrete Picard condition is satisfied. 0 1997 by Elsevier Science inc.
Keywords: computer-aided
design, polynomial
basis conversion, Tikhonov regularisation,
1. Introduction Polynomial basis conversion is an important industrial requirement because it allows the exchange of data between computer-aided design systems that use different mathematical representations of curves and surfaces. Even when the conversion is theoretically exact, numerical difficulties may occur when the linear algebraic equation that defines the transformation is ill-conditioned, in which case the solution is computationally unreliable. Conventional methods such as Gaussian elimination and Cholesky decomposition cannot be used, and special techniques, called regularisation methods, must be employed. One popular regularisation method is truncated singular value decomposition, and its effectiveness in regularising a specific class of ill-conditioned linear algebraic equation has been considered.’ Another regularisation method, known as Tikhonov regularisation in standard form (henceforth called Tikhonov regularisation), is considered in this paper, and it is shown that this method is closely related to truncated singular value decomposition. An example of the use of Tikhonov regularisation for the transformation of a Bezier curve segment into its B-spline’ form is given.
Address reprint requests to Dr. Joab R. Winkler at The University of Sheffield, Regent Court, Department of Computer Science, 211 Portobello Street, Sheffield S14DP, UK. Received
27 January
1997; revised
5 June
1997; accepted
Appl. Math. Modelling 1997, 21:651-662, October 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
3 July 1997.
B-splines
It is shown in Section 2 that the condition number of the coefficient matrix of the linear algebraic equation that defines the transformation of a Bezier curve segment into its B-spline form increases rapidly as the degree of the polynomial increases. Although the transformation between piecewise polynomial curves is of greater industrial relevance, it is advantageous to consider first the transformation of individual curve segments because attention can be restricted to the application of regularisation to the simplest class of polynomial basis conversion problem. This is an important point because Tikhonov regularisation and truncated singular value decomposition are closely related to the numerical solution of Fredholm integral equations of the first kind, and these equations, and their discretised equivalents, are known apriori to be ill-conditioned. However polynomial basis transformation problems are not associated with integral equations, and thus the linear algebraic equation that defines the basis transformation may be ill-conditioned. Tikhonov regularisation and the discrete Picard condition are introduced in Section 3, and it is shown in Section 4 that if this condition is satisfied, the equation is ill-conditioned but that the converse is not true. This leads to the generalised singular value decomposition, and this is considered in Section 5. Regularisation must, of necessity, introduce an error, and this error is quantified in Section 6. It is shown that if the discrete Picard condition is satisfied, this error is guaranteed to be small, that is, the regularised solution is a good approximation to the theoretically exact solution.
0307-904x/97/$17.00 PII s0307-904x(97)00081-4
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
An important issue in regularisation is the degree of smoothing that is imposed, because if too little smoothing is applied, the computed solution may be sensitive to small perturbations in the input data, but if too much regularisation is applied, the computed solution is too smooth, that is, some of the desired underlying solution has been removed from the computed solution. A graphical method for the determination of the degree of smoothing to be imposed is described in Section 7, and Section 8 contains an example of the conversion of a Bezier curve segment of degree 12 to its B-spline form.
2. Polynomial
basis transformation
The transformation between a Bezier curve segment its B-spline form, both of degree n, is’ JH=GF
and
(1)
where J, G E !Rmxm, H, FE !Rmmx2,m = n + 1, and rank J = rank G = m. Explicit forms for the elements of G and J are in Ref. 1. The matrices H and F store the coordinates of the vertices of the control polygons of the B-spline curve and Bezier curve, respectively. The conversion between the Bernstein and B-spline polynomial bases can be performed by the generalized Oslo algorithm, but it is noted in Ref. 3, page 152, that “. . . it is possible to determine that the generalized Oslo algorithm is stable in whether the algorithm is some cases.. . . To determine unstable for given knot vectors.. . seems like a more difficult task.” It is shown in Ref. 1 that the condition number of J increases rapidly with the degree of the polynomial, and it follows that there exist situations in which the basis conversion problem is ill-conditioned. It will be shown that when the transformation is performed numerically, these ill-conditioned situations are readily apparent, and that in some, but not all, situations, Tikhonov regularisation leads to a computationally reliable solution that is a good approximation to the theoretically exact answer. It is shown in Ref. 1 that it is necessary to rewrite equation (1) as Jh, = (GF)l
Jh, = (GF);?
where H = [h, h,] and GF = [(GF), these equations separately. It follows linear algebraic equation
(2a,b)
GF),l, and that
3. Tikhonov
regularisation
Since equation (3) is ill-conditioned, it is necessary to regularise the equation, that is, make it less sensitive to noise. Most numerical methods for the regularisation of a linear algebraic equation seek to overcome the problems of ill-conditioning by replacing it with a “nearby” wellconditioned problem whose solution approximates the desired solution and is more satisfactory than is the simple solution. The latter goal is achieved by the addition of more information about the sought solution, often that the computed solution should be smooth. These methods are called regularisation procedures, and they always include a parameter that controls the degree of smoothing that is imposed on the computed solution. For example, the regularisation parameter in truncated singular value decomposition is a discrete index that controls the number of nonzero singular values that are retained in the expression for the solution x. The solution that is obtained from a regularisation procedure must possess the following two properties: (i) It must be a good approximation to the theoretically exact solution; and (ii) it must be less sensitive to noise than is the simple (direct) solution of the equation that is being regularised. It will be shown that both properties are satisfied if it is known a prior-f that the theoretically exact solution satisfies the discrete Picard condition. The regularised solution of Tikhonov regularisation4-6 is defined by the unique solution of min {IlAx - bll’ + ~u~llx1l~l X
(4)
where (Y is the regularisation parameter and controls the degree with which the constraint on ((x1(is imposed on x = x( CY). 11(1 denotes the 2-norm. The solution x(a) of this minimisation problem is given by the solution of the equation 64% + (Y21)X(a) = ATb
(5)
This equation is most conveniently considered by using the singular value decomposition7 of A, given by USVT where U, S, V E sm”““, U and V are orthogonal matrices, and S is the diagonal matrix whose diagonal elements, sii = si, i = l..m, are the singular values of A, arranged in nonincreasing order. It follows that
solve
the general
x(cy)=v(s2+(y21)-1sc= 5
L-v.
ci
i=l S12+(Y2 Si ’
Ax=b
(3)
can be considered. A E smxm, x, b E %“‘, and A is nonsingular and ill-conditioned. Equation (2a) defines the transformation between the x-coordinates of the vertices of the control polygons of the Bezier and B-spline forms of a curve, and equation (2b) defines the transformation between the y-coordinates of the vertices of the control polygons of the curves.
652
Appl.
(6)
Math.
Modelling,
1997,
Vol. 21, October
where c = {ci} = UTb, and v, is the ith column of V. U and V are the matrices of the left and right singular vectors, respectively. The theoretically exact solution of equation (3) is
x(0) = vs-‘c
m ci = c ;vi i=, I
(7)
Tikhonov regularisation Equation
X(a)
(6) can be written
in the form
is solved. The solution
Y,(,,%‘~
= ~ i=l
is given by
2
1
&W=*=
2’
1
1+
i = l..m
(9)
” i s, 1
It is clear that
=
J. R. Winkler
is
ci + SC; 5 ~ v L s;
i
i=l
i
Since the ratios Ici + 6cil/si and Ic,l/si usually span several orders of magnitude, it is more convenient to consider the Picard ratios logic; + 6ci(/si and logIciJ/s, flog = log,,). Since it is assumed that the desired solution satisfies the discrete Picard condition, lcil + 0 as i + m. However the components lc, + 6ciI/si typically decrease to a minimum level that is determined by the noise and then increase, as shown in Figure 2. In particular, if p is the index at the minimum point of the curve in the figure, then
1
if (Y<< si
lcil B lSc,l
l/2
if a=si
lcil=lSc,I
if (YB si
lcil e 16cil
i 0
conversion:
of this equation
x(0) + 6x(O) =
(8)
I
where the filter factor g,(a)
g,(a)
for polynomial
for
i
for
i=p i >p
for
Figure 2 shows that The variation of gi( CY>with (Y/S; is shown in Figure 1. As cr increases, the filter factor g;(c~y) decreases, and the terms in equation (8) that are associated with the small singular values are removed from the expression for X((Y). It follows that if X((Y) is to be a good approximation to x(O), it is necessary that x(O) be dominated by the large singular values, which implies that the ratio lcil/si must decrease towards zero as the index i increases. Since the singular values are arranged in nonincreasing order, it is necessary that the components lcil decrease towards zero more rapidly than do the singular values. This is the discrete Picard condition.*
x(0) + 6x(O) =
+vm
(11)
m
and thus the simple solution of equation (10) is dominated by the noise and all knowledge of the exact solution is lost. However the regularisation parameter CY can be chosen such that g,(ff)
=
1 0 i
i=l..p-1 i=p..m
in which case Definition: The discrete Picard condition. The unperturbed right-hand side b of equation (3) satisfies the discrete Picard condition if for all numerically nonzero singular values si, the coefficients lcil decay towards zero, on average, more rapidly than do the singular values. In practical situations the right-hand side vector b is corrupted by noise such that the perturbed equation A(x+6x)=b+6b
(10)
This solution is a good approximation to x(O) if the decay of the ratio Icil/sj is sufficiently rapid. It is seen that the satisfaction of the discrete Picard condition guarantees that the regularised solution x(cu) + OX is a good approximation to the theoretically exact solution x(01.
Picard ratio Ii
filter factor
,
i
P
Figure
1.
The variation
of the filter
factor
g,(a)
with
a/s,
Figure 2.
Appl.
Math.
The variation
Modelling,
of the Picard
1997,
ratio logjc,+
Gc,Vs, with
Vol. 21, October
653
i.
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
The difference between truncated singular value decomposition and Tikhonov regularisation is the manner in which the singular values are removed; in the former method, they are removed in discrete steps by reducing the value of an integer, but in the latter method they are reduced towards zero gradually by increasing the value of the continuous parameter CY.The close relationship between these two regularisation methods will be shown by example in Section 8. Tikhonov regularisation and truncated singular value decomposition are based on the following heuristic:5
Heuristic. The number of oscillations in the left and right singular vectors ui and vi tends to increase with increasing i. The number of oscillations in a vector z is a measure of its spectral components; as the number of oscillations increases, the components of z exhibit a greater number of sign changes. The heuristic6 is clearly evident in Figure 3 in Hansen (1992j6, taken from an example of image restoration. It is clear that when the heuristic is satisfied, the regularised solution is smooth because it is a linear combination of nonoscillatory vectors. This oscillatory feature is satisfied when equation (3) is derived from the discretisation of a Fredholm integral equation of the first kind using orthonormal basis functions.6 However this property may not be satisfied for polynomial basis conversion problems. For example, Figure 3 shows the left singular matrix U for the transformation matrix P between the power and Bernstein polynomial bases for univariate polynomials of degree twelve, where Py = r and y and r store, respectively, the coefficients of the polynomial in the Bernstein and power bases.’ Similarly Figure 4 shows the right singular matrix V of P. Figures 5 and 6 show, respectively, the matrices U and V for the transformation matrix J in equation (1). It is clear from the figures that the heuristic is not satisfied for the transformation matrices P and J. It is interesting to note that only the left singular matrix of J has a strong pattern. Even though the
column Figure 4. matrix P.
index The right singular
matrix V for the transformation
heuristic may not be satisfied for polynomial basis conversion problems, it is shown in Section 8 that Tikhonov regularisation is effective in regularising a specific class of ill-conditioned linear algebraic equation.
4. The stability
of linear algebraic
equations
It was shown in the previous section that if the discrete Picard condition is satisfied, the solution of equation (3) is susceptible to the effects of random noise, from which it follows that the equation is ill-conditioned. This qualitative statement is quantified in this section, and it is shown that Tikhonov regularisation can only solve a particular class of ill-conditioned linear algebraic equation. The ratio relative error in x( (Y> relative error in b
U
column Figure 3. matrix P.
654
index
column
”
The left singular
Appl.
Math.
matrix
Modelling,
U for the transformation
1997,
Vol. 21, October
Figure 5. matrix J.
index The left singular
matrix U for the transformation
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
It is important to note that this measure cannot be used if the discrete Picard condition is satisfied. Consider the effective condition number of equation (10)
K,&,
IIC+ 6cll
b + 6 b) = L s,
Ils-‘(c
+ Sc)II
Since it is assumed that the discrete Picard condition is satisfied and that llbll x=-IlSbll, it follows that IIS-‘(c + 6c)ll = llS-'Sell = ISc,l/s, and Ilc + 6cll= llcll = Icll. Thus ~,,(-kb column Figure 6. ma’trix J.
index The right singular
matrix V for the transformation
is called the sensitivity function, and a formula for it is developed in Theorem 1. This ratio is important for the determination of the error magnification that may occur.
This shows that the effective condition number is dominated by the effects of noise, and it therefore provides totally misleading answers. However the measure in equation (15) is used for the derivation of the following important result. Further study of the discrete Picard condition is facilitated by postulating the following simple relationship between the magnitude of the components ci and singular values si
Theorem 1. Let A E ‘Smxm,x, b E S”, and rank A = m. The sensitivity function 5(x( a), b) for equation (5) is
Jcil =
sf,i = 1..m,
830
The parameter 0 determines ratio IciI/si with i:
S(x(a),b) =
ICll + ab) = ~sc,l
llW + a21)%cII IISCII
Ml ll(S2 + c2I)%ll
(12)
0 < 0 < 1 j IciI/si increases
(16) the dependence
of the
monotonically, i= l...m
where
SC = U% b.
O=l*(ci(/si=l,i=l...m
Proofi It follows from equation
1 < 0 * Icil/si decreases (6) that
8x((Y) =V(S2 + cr21)-1s6c The relative
error in X((Y) is given by ll(S2 + cY21)-k6cJI
IlSx((r)ll Ax(a)
= IlX(Q)II
since V is an orthogonal given by Ab = ll6bll -=llbll
=
ll(S2 + (y21)-1Scll
matrix. The relative
(13)
error in b is
IlScll (14)
llcll
Equations (13) and (14) lead directly to equation (12). This completes the proof. The effective condition number of equation (3) is obtained by setting CY= 0 in equation (12) and taking the maximum value over all perturbations SC.” This yields +&‘&b)
I
llcll
= sc~a& Sk(O), b) = y- m IIS-‘cll
(15)
monotonically,
i = 1. . . m
This model distinguishes clearly between the right-hand side vectors that do and do not satisfy the discrete Picard condition. It should be noted that although the model in equation (16) is simpler than is likely to occur in a real problem, it is adequate for the sequel. It is shown in Ref. 11 that when this model of c is substituted into equation (15), equation (3) is ill-conditioned if the discrete Picard condition is satisfied (0 > l), but that the converse is not true. This emphasises the importance of’establishing that the discrete Picard condition is satisfied and shows that Tikhonov regularisation can only solve a particular class of ill-conditioned linear algebraic equation. The satisfaction of the discrete Picard condition can be verified by plotting the terms Ici + &,I and si against i. If the coefficients Jci + &,I decay towards zero more rapidly than do the singular values si, it can be assumed that the discrete Picard condition is satisfied. However if the coefficients level off at a value that is sufficiently small, then this value is an estimate of the noise level in c. The imposition of mild statistical assumptions about the noise 8 b sheds further insight into the discrete Picard condition, and moreover, it reveals the degradation in the solution x(a) that may occur if this condition is not satisfied. The square of the sensitivity function of equa-
Appl.
Math.
Modelling,
1997,
Vol. 21, October
655
Tikhonov
regularisation
for polynomial
conversion:
J. R. Winkler
tion (5) is
into the expression
for $(cr) it follows that
s%(a),b) =
KS2 + a~I)-‘S6cl12
llCl12
l16cl12
i=l
(Sf+(Y2)2
ll(S2+ a21)- ‘Scl12 _
j=l
6cV
E
(17)
(Sf+a2)3
I I
5
i-1 (S,‘+ (Y2)3 j=l Theorem 2. The derivative S 2(x( (Y>,b) satisfies
with
dS(x(a ), b)
S(x(a),b)
da2
= -
respect
to
a2
of
which simplifies
ll~cl12
a’)’
to
Since this expression written as m-l
+ (Y~I)-~SSC
tic(Y)=
(19a, b)
m
c
i = j, it can
is zero when
6c~c;s,‘s;(
s’
-
be
si”)
c
i=l
j=i+l
(S,2+a2)3(S,~+CX2)3
+ a21)-3Sc
T3 = 6cTS(S2
+ a21)-3S6c
m-l
m
j=l
c i=j+l
+c
(19c, d) T4 = crS(S2
+
T42
where
T, = crS(S2
(Sf
11412U,T2 - TJ,) (18)
T, = 6cTS(S2
+:
6c~c;s;s;( (Sf
+
si’ -
c12)3(S;
+
ST) CT’)’
+ LU~I)-~SC
Proofi
It follows from equation (17) that it is necessary evaluate the derivative of the expression
to
On making the substitutions i + j, j + i in the second term on the right-hand side, this expression becomes
ll(S2 + (Y21)-1s6c112
ll(S2+ a2K1Scl12 and then multiply the result by (l~~~~/ll6cJl~. This completes the proof. Let 6 b be a random vector with a diagonal covariance matrix proportional to the identity matrix
If the statistical assumption about 6c in equation (21) is substituted into this expression, it follows that the expected value of $4 a 1 is m-l E{+((Y)}=E~
c i=l
E(6b 6bT) = ~~1
m
s~sj’(s~-q(c;-c~)
c j=i+l
(Sj
+
a2)3(S,!
+
Cf2)3
(20) (22)
where I E !RH”““,c2 is the variance of each component 6 b and E is the expectation operator. It follows that E{ 6c 6cT} = ~~1
of
(21)
Theorem 3. If the right-hand side vector b in equation (3) satisfies the discrete Picard condition and the perturbation 6 b satisfies equation (20), then Tikhonov regularisation reduces, on average, the effect of noise on the solution of equation (3). Proofi
It is required
that
+(a)=T,T2-T3T4<0
because this implies that the derivative of the sensitivity function is negative. Also, it follows from equation (21) that 116~11~ = me2. On substituting equations (19a)-Cd)
656
Appl.
Math.
Modelling,
1997,
Vol. 21, October
If the discrete Picard condition is satisfied, the average value of 1+4(a) is negative and Tikhonov regularisation reduces, on average, the effect of noise on the solution of equation (3). However if the discrete Picard condition is not satisfied, the term on the right-hand side of equation (22) may be positive, in which case Tikhonov regularisation degrades the solution of equation (3) because the effect of the noise 6 b is amplified. This completes the proof. Theorem 3 shows that if lcjl < lcil for j > i, the effect of Tikhonov regularisation is a reduction in the sensitivity of x(a) to noise. However the satisfaction of this condition does not guarantee that the discrete Picard condition is satisfied, and thus the error between X((Y) and the theoretically exact solution may be large, such that x(a) is of little use. It was noted in Section 1 that this paper is concerned in standard form, and it has with Tikhonov regularisation
Tikhonov regularisation been shown that this form of regularisation can only stabilise a particular class of ill-conditioned linear algebraic equation. It is shown in the next section that the generalisation of this form of regularisation enables other classes of ill-conditioned linear algebraic equation to be stabilised.
for polynomial
Theorem 4. The generalised singular value decomposition. Let the matrices A and L be defined in equations (31 and (23b), respectively. Then there exist matrices R E 8 mxm,T E %PpxPwith RTR = I,,,,TTT = I, (the identity matrices of orders m and p, respectively), and a nonsingular matrix X E Ytmxm such that A=RNX-’
5. The generalised
singular value decomposition
conversion: J. R. Winkler
L = T[M,O]X-’
(26a, b)
where
It is clear from equation (4) that Tikhonov regularisation places a constraint on Ilxll. It is possible to place constraints on the norm of other functions of x, such as the norm of its first derivative. Thus equation (4) can be written more generally as
N = diag[ r1 . . ap
1.. 11 E %I”““, (27a)
M=diag[~l..~plE%pPXP and
min{JIAx - b112+ a2~jLx~~2)
(23a)
x
where L E %Ppxm,m >p,rank(L)
=p, N(A) n N(L) = (0)
ai*+j$=l,
(27~)
i = l..p
(23b) The generalised IlLxll is a semi-norm and L is typically a discrete approximation to a derivative operator. N(X) denotes the null space of X, and it is shown in Theorem 5 that the condition that the intersection of the null spaces of A and L is the zero vector guarantees that equation (23a) has a unique minimiser. Examples of L include L, E MC” - ‘jxrn and L, E !@‘-2)xm, given by 1 0
-1
0
0 0
l-1
0 0
0 0
1
&,
0 1.00
0 00.
-1 L, =
0 01. 2
0 0 0 0
l-l
i = l..p
If L = I,, then X-’ = M-ITT and A = R(NM-‘)TT, and the generalised singular values y, reduce to the singular values of A (with pi # 0, but in reverse order. It is convenient to partition the matrices R, N, and X as
-1
0
0 0
2
-1
0
2 -1 .
-1
[
0 0 0
0 0 0
2.
0 -1 I
-1
L, and L, are, respectively, proportional to approximations of the first and second derivative operators on a uniform net. It follows that the minimiser of equation (23a) is given by the solution of the equation (ATA +
a2LTL)X(a)
=
ATb
I
NP O
N=
IO
0
0. I
0 0
-1
-1
0 0
01-l . 0
2
&,I
0
-1 0.
-1
values ‘y, are defined by
Pi
R=[R, L,=
singular
(24)
X=[X,
x,1
(28a, b, c)
where R, and X, have p columns, Np E %PxP and o = m -p. Equations (26a), (26b), (27a1, (27b), (27~1, (28a), (28b1, and (28~) enable the matrix A’ to be simplified to A’ = X,FN;R;
+ X,R;
where Nl is the pseudo-inverse the diagonal matrix
of Np, and F E %Pp*P is
It follows that the expression written as
in equation
(25) can be
and thus X(a) = (ATA + (y2LTL)-lATb
(25)
The matrix A’ = (ATA+ (Y*LTL>-’ AT is most conveniently analyzed by the generalised singular value decomposition’ of the matrix pair (A,L). This decomposition reduces to the singular value decomposition of A when L = I.
x((y)’
2
5 yi i=, ++a*
qTb -xi a,
+
E (rTb)xi i=p+l
(29)
where ri is the ith column of R. It is noted that the generalised singular values are arranged in nondecreasing order, which must be compared with the order in which
Appl. Math.
Modelling,
1997,
Vol. 21, October
657
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
the singular values si are arranged. The second term on the right-hand side of equation (29) is independent of the regular&&on parameter ry, and the vectors xi, i =p + l..m lie in the null space of L. This follows immediately from equations (26b) and (28~) because LX=L[X,
X,1 =T[M
Y?
i = l..p
+Cr2’
‘yi
As the regularisation parameter (Y increases, the terms in equation (29) that are associated with the small generalised singular values are filtered out from the expression for x(a). It follows that the expression in equation (29) is a good approximation to the theoretically exact solution x(O) if the dominant components of x(0) are associated with the large generalised singular values yi, that is, the expression in equation (29) is dominated by the components that are associated with large values of the index i, and this implies that it is necessary that
lr’bl
-=-ai
1 IrTbl +O Pi
asi-
Yi
Since pi + 1 as i + 1, it follows approximation to x(O) if
lr,% -ro
rank[ t] where
01
and thus LX, = 0. Equation (29) shows that the solution x(a) is expressed as a linear combination of the vectors xi and that the filter factor g&a) in equation (9) is replaced by the filter factor
2
is y = 0, and this implies that the rank of the coefficient matrix of this homogeneous equation is m. Since
that
x(a)
is a good
= rank[ aAL]
(Y is an arbitrary rank1 iL]
and
nonzero
= raok{MT
rankA = rankATA
constant,
it follows that
aLTl[ :L]}
= rank(ATA + cy2LTL) = m
(31)
( = ) If rank(ATA + (Y2LTL) = m, it follows from equation (31) that the only solution of equation (30) is y = 0, and this implies that N(A) n N(L) = (0). This completes the proof. If the ratio Icil/si does not decrease towards zero as the index i increases but equation (3) is ill-conditioned because 8 = 1 in equation (16), it is necessary to select a matrix L such that the coefficients lrTb[ decay towards zero more rapidly than do the generalised singular values 7,. The geometric interpretation is that a set of linearly independent vectors xi is sought such that the dominant components of x(O) lie in the subspace of St” that is spanned by those vectors xi that are associated with the large generalised singular values ‘yi. An example of the use of the matrix L, for the regularisation of the equation that arises from the discretisation of a Fredholm integral equation of the first kind is in Ref. 8. It is concluded that the form of regularisation in equation (5) is a particular example of a more general class of regularisation methods, and that this larger class must be considered when solving ill-conditioned linear algebraic equations. The selection of the matrix L is a major research issue, and the incorrect choice yields an erroneous answer.
asi-+l 6. The regularisation
Yi
Since yi -+ 0 as i --) 1, it follows that the coefficients lr’bl must decay towards zero more rapidly than do the generalised singular values yi. This is the discrete Picard condition for equation (241, and its satisfaction is crucial for the correct use of Tikhonov regularisation in general form.
Theorem 5. Let the matrices A, L be defined (3) and (23b), respectively. Then
in equations
error
It is shown in Ref. 1 that truncated singular value decomposition introduces a bias because the regularised solution does not reduce to the theoretically exact solution in the absence of noise. It is shown in this section that a similar result is valid for Tikhonov regularisation. The total absolute error in the solution of equation (3) is x(0) -(x((Y)
+ 6x((Y))
=x(O) - (ATA + (r21)-‘AT(b + 6b) N(A) n N(L) = (0)
ts
rank(A’A + (Y’ L’L) = m
Proofi (a) If N(A) n N(L) = {O}, the only solution equation
[
658
I
A L Y=o
Appt.
and this simplifies
Modelling,
1997,
+ s b)
to
of the
(30)
Math.
= X(O) - (ADA+ n2r)-‘~T(Ax(0)
Vol. 21, October
The first term in this expression is the error, in the absence of noise, due to the regularisation procedure, and the second term includes the effects of noise and the
Tikhonov
regularisation procedure. It is seen that as 6 b -+ 0, the absolute error does not approach the zero vector, which implies that Tikhonov regularisation introduces a bias in the solution X((Y). The regularisation error e(a) is defined as the normalised value of the total absolute error in the absence of noise e(a)
=
=
regularisation
for polynomial
IA
B
I C
I
ll(I - (ADA+ (Y~I)-~A~A)x(O)II
IIx(O)II
Figure
(16) is
Pro08 See Ref. 5, Theorem
3.1. If (Y is sufficiently small compared to s1 and 8 > 1, X((Y) is guaranteed to approximate x(O), such that the regularisation error e(cr> is small, and the larger the value of 8, the better the approximation. However if the discrete Picard condition is not satisfied (0 < 0 Q l), x( LY)is a poor approximation to x(O). As the magnitude of the perturbation increases, the required value of (Y increases because a larger number of the small singular values must be removed from the expression in equation (6). It follows that 0 must increase in order to maintain a small value of e(a).
7. The L-curve The success of Tikhonov regularisation in regularising equation (3) is dependent upon the value of (Y,even if it is known that the theoretically exact solution satisfies the discrete Picard condition. If the chosen value of cr is too small, x(cw> may be corrupted by noise, but if it is too large, a significant portion of the desired solution may be removed from X((Y). It is shown in Refs. 12 and 13 that the L-curve is a good method for the calculation of the optimum value of LY. Equation (4) shows that it is desired to minimise IlAx - bll while constraining Ilx(cr>ll, from which it follows that a method for the optimum determination of a should measure both these quantities. The L-curve is a parametric plot of logllx(cy)II against logllAx(cy) - bll. It is shown in Refs. 12 and 13 that if (i) the unperturbed right-hand side vector b of equation (3) satisfies the discrete Picard condition, (ii) the perturbation 6 b is a zero mean random variable that satisfies equation (20), and (iii)
J. R. Winkler
log magnitude
IIx(0) -x(a)II IIx(O) II
Theorem 6. If the simple model for lcil in equation assumed, the regularisation error e( (Y> satisfies
conversion:
7.
A
schematic
D
E
log residual
diagram
of an L-curve
II6 bll < llbll, the sequence of points that define the L-curve assume the shape shown in Figure 7. The value of (Y increases as one progresses along the curve from A to E. In particular, the section ABC of the curve is characterised by solutions that are dominated by the effects of noise, and the portion CDE is characterised by solutions that have been smoothed too much. The optimum value of (Y is its value at the corner C because it provides the best value that minimises both IlAx - bll and Ilx(cy)II. It is shown in Refs. 12 and 13 that the curve in Figure 7 is a combination of two separate L-curves (Figure 8): (9 the L-curve when b = b, is defined by zero mean random noise that satisfies equation (20); and (ii> the L-curve when b = b, has no noise component and satisfies the discrete Picard condition. When (Y is small, the composite L-curve is dominated by the effects of b,, but as it increases, the curve is dominated by contributions from b,. The region between these two extremes defines the corner C, as shown in Figure 7.
8. Example This section contains an example of the conversion of a curve segment of degree 12 (m = 13) from the Bernstein basis to the B-spline basis. Since high-degree polynomials (up to degree 15 or 21) are used in some computer-aided design systems, l4 there is a need to obtain a computationally reliable answer to this ill-conditioned basis transformation problem. The notation of equation (3) that was used in the previous sections is replaced by the notation of
log magnitude _______ : [ (a)
(b)
I Figure
8.
log residual The
L-curves
that satisfies
the discrete
Math.
Modelling,
Appl.
for (a) only Picard
noise
and (b) only
signal
condition.
1997,
Vol. 21, October
659
Tikhonov regularisation
for polynomial
Bezier
conversion: J. R. Winkler
curve
0.28
0.27
\
-81 Figure (2a).
0.26
0.25
t=o 0.24
...
Figure 9. polygon.
...... 0.29 0.3
0.31
0.32
0.33
0.34
The Bezier curve and the convex hull of its control
equations (2a) and (2b). In order to facilitate comparison between Tikhonov regularisation and truncated singular value decomposition, the data of Example 2 in Ref. 1 is used. Figures 9 and 10, respectively, the Bezier and B-spline curves and the convex hulls of their control polygons. Although the curves are necessarily identical, their control polygons differ. The vertices of the control polygons in these figures, and Figures 13 and 16(b), are denoted by black dots. Zero mean Gaussian noise was added to the right-hand side of equations (2a) and (2b), and its variance was adjusted so that
ll(G~h II II 6(GF)l II
11 (GFh 11 = 2.12 x lo3
B-spline 0.5
II6(GF)2 II
11.
(a) loglc,I/s,
and (b) logic,+ 6c,I/s,
Figures 11 and 12 show the variation of the ratios loglcil/si and logic, + 6ciI/s, with i for equations (2a) and (2b), respectively. They show that the simple solutions of equations (2a) and (2b) is unsatisfactory. This is confirmed by Figure 13, which shows the convex hull of the control polygon of the B-spline curve when the random noise is applied to the right-hand side of the equations and no regularisation is used. Comparison of Figures 10 and 13 shows that all knowledge of the exact solution is lost and that the noise dominates the solution. Figure 14 shows the B-spline curve in this circumstance, computed from
p(t)
+
6p(t) = t(GF + 6(GF)) = tJ(H + 6H)
and comparison of Figure 9 and 14 shows that the change in the curve is negligible. Since 6 p(t) = t 6(GF) it follows that
11 Sp(t)ll,
< lltllxll S(GF>h
=
11 S(GF)lI=
= 1.07 x lo3 for 0
curve
F
4
Picard ratio
(b)
2
Figure 10.
The B-spline curve
and the convex
hull of its con-
Appl.
Math.
Modelling,
1997,
Vol. 21, October
,---,*
-8 t Figure 12.
trol polygon.
660
for equation
(a) loglc,I/sj
and (b) logic;+ Sc,I/s,
for equation G’b).
Tikhonov
B-spline
basis
regularisation
for polynomial
conversion:
J. R. Winkler
log magnitude
log magnitude
Figure 13. The convex hull of the control polygon of the B-spline curve when noise is applied to equations (2a. b).
is small. (Note that the infinity norm clidean norm is used.) For 0 < t Q 1 Ilp(t)ll,
= Il[x(t)
ywlll,
and not the Eu-
= 0.34
and this occurs at t = 0.8. It follows that a measure relative change in the curve due to the perturbation
of the is
II 6pCt) IL < 4.38 x 1O-4 IIp(t) IL which is small. Figures 13 and 14 show that it is necessary to distinguish between the integrity of the transformed curve and the integrity of the transformed convex hull. Figure 15(a) and (b), shows the discrete points that define the L-curves for equations (2a) and (2b), respec-
B-spline
(b) Figure 15. (a) The L curve for equation (2a). -, Tikhonov regularisation, * truncated singular value decomposition. (b) The L-curve for equation (2b). -, Tikhonov regularisation, . truncated singular value decomposition.
basis
tively. The results that are obtained when truncated singular value decomposition is used, reproduced from Ref. 1, are shown, and it is seen that there is excellent agreement between the L-curves for the two methods. Figure 16(a) and (b) show the B-spline curve, and the B-spline curve and convex hull of its control polygon, for (Y, = (Y*= 1.2 x 10e4, where (Y, and (Y* denote the optimum regularisation parameters, as determined from Figure 15(a) and (b), respectively. It is seen that Tikhonov regularisation is very effective in restoring the quality of the convex hull of the B-spline curve segment. Very good answers were also obtained for (Ye= CY*= 1.2 x lo-“. The magnitudes of the residual and solution at these two values of (Y are
0.301
log,,II Jh,( a) - (GFh 11= -3.841 Ly= 1.2 x lo-“:
0.28.
log,,llh,(a)ll
= -0.177
log,& Jh,( LY)- (GF)z 11= - 3.597 log,,ll h,( cr) II = -0.24?
0.26.
log,,Il Jh,(c.u) - (GF), II = -3.518 cX= 1.2x
0.24.
10-4:
t=o
log,,llh,(a)ll
log,,tI Jh,( a, 1 - (GF)2 II = - 3.496 loglollhz(a)ll
0.28
0.30
0.32
= -0.177
= -0.248
0.34
Figure 14. The B-spline curve when noise is applied to equations (2a. b).
Comparison of these results and those of Example 2 in Ref. 1 confirm that if the discrete Picard condition is
Appl.
Math.
Modelling,
1997,
Vol. 21, October
661
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler Acknowledgment The author acknowledges the support of the UK Engineering and Physical Sciences Research Council (EPSRC) by an Advanced Research Fellowship, reference number B/93/AF/1703.
0.28. 0.27.
References 0.25. 0.24
0.29
0.30
0.31
0.32
0.33
1. Winkler, J. R. Polynomial basis conversion made stable by truncated singular value decomposition. Appl. Math. Mode&g, in press 2. Boehm, W., Farin, G., and Kahmann, J. A survey of curve and surface methods in CAGD. Compui. Aided Geo. Design 1984, 1, l-60 3. Lyche, T., Morken, K., and Strom, K. Conversion between Bspline bases using the generalized Oslo algorithm, in Knot Insertion and Delefion Algorithms for B-spline Curves and Surfaces. Eds. R. N. Goldman and T. Lyche, SIAM, Philadelphia, 1993 4. Hanke, M., Hansen, P. C. Regularization methods for large-scale problems. Survey Math. Ind. 1993, 3, 253-315 5. Hansen, P. C. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM J. Sci. Staf. Comput. 1990, 11, 503-518 6. Hansen, P. C. Numerical tools for analysis and solution of Fredholm integral equations of the first kind. Inverse Problems 1992,
0.34
(a)
8, 849-872 7. Golub,
G. H. and Van Loan, C. F. Matrix Computations. Johns Hopkins University Press, Baltimore, 1989 8. Hansen, P. C. The discrete Picard condition for discrete ill-posed problems. BIT 1990,30,658-672 9. Farouki, R. T. and Rajan, V. T. On the numerical condition of polynomials in Bernstein form. Comput. Aided Geo. Design 1987, 4, 191-216 10. Winkler, J. R. The sensitivity of linear algebraic equations. Proceedings
I (b) Figure 16. (a) The ELspline curve for a1 = CY~= 1.2 X 10m4. (b) The B-spline curve and convex hull of its control polygon for aI = (~z = 1.2 x 10-4.
of the Fifth SIAM
put. 1993,
satisfied,
truncated
662
Appl.
Math.
value decomposition yield very similar solutions
singular
Modelling,
1997,
and (see
Vol. 21, October
on Applied
Linear Algebra,
14, 1487-1503
14. Lachance,
Tikhonov regularisation Ref. 5, Theorem 3.2).
Conference
Snowbird, Utah, USA, June 15-18 1994, pp. 279-283 11. Winkler, J. R. The condition number of a matrix and the stability of linear algebraic equations. Department of Computer Science, The University of Sheffield, internal research report CS-97-05 12. Hansen, P. C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Reu. 1992,34, 561-580 13. Hansen, P. C. and O’Leary, D. P. The use of the L-curve in the regularization of discrete ill-posed problems, SIAM J. Sci. Com-
mial curves. Algorithms
Shrivenham, pp. 125-133
M. A. Piecewise
polynomial
approximation
of polyno-
Proceedings of the Second International Conference on for Approximation, Royal Military College of Science,
UK, July 1988, eds. J. C. Mason
and M. G. Cox,
UK
Polynomial basis conversion is an important industrial requirement for high-integn’ty data exchange between computer-aided design systems. These conversions may be numen’cally ill-conditioned, and this paper considers the use of Tikhonov regulatisation in standard form for their stabilisation. It is shown that the method yields a good approximation to the theoretically exact solution if it is known a priori that the desired solution satisfies the discrete Picard condition. An example of the use of Tikhonov regularisation for the stabilisation of the conversion between the Bernstein and B-spline polynomial bases is given. This basis transformation can also be performed by the Oslo algorithm, but it is difficult to determine whether this algotithm ir unstable for given knot vectors. However when this basis transformation is performed by Tikhonov regularisation, it is readily determined whether the linear algebraic equation is unstable, and if it is, regulan’sation follows immediately if the discrete Picard condition is satisfied. 0 1997 by Elsevier Science inc.
Keywords: computer-aided
design, polynomial
basis conversion, Tikhonov regularisation,
1. Introduction Polynomial basis conversion is an important industrial requirement because it allows the exchange of data between computer-aided design systems that use different mathematical representations of curves and surfaces. Even when the conversion is theoretically exact, numerical difficulties may occur when the linear algebraic equation that defines the transformation is ill-conditioned, in which case the solution is computationally unreliable. Conventional methods such as Gaussian elimination and Cholesky decomposition cannot be used, and special techniques, called regularisation methods, must be employed. One popular regularisation method is truncated singular value decomposition, and its effectiveness in regularising a specific class of ill-conditioned linear algebraic equation has been considered.’ Another regularisation method, known as Tikhonov regularisation in standard form (henceforth called Tikhonov regularisation), is considered in this paper, and it is shown that this method is closely related to truncated singular value decomposition. An example of the use of Tikhonov regularisation for the transformation of a Bezier curve segment into its B-spline’ form is given.
Address reprint requests to Dr. Joab R. Winkler at The University of Sheffield, Regent Court, Department of Computer Science, 211 Portobello Street, Sheffield S14DP, UK. Received
27 January
1997; revised
5 June
1997; accepted
Appl. Math. Modelling 1997, 21:651-662, October 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
3 July 1997.
B-splines
It is shown in Section 2 that the condition number of the coefficient matrix of the linear algebraic equation that defines the transformation of a Bezier curve segment into its B-spline form increases rapidly as the degree of the polynomial increases. Although the transformation between piecewise polynomial curves is of greater industrial relevance, it is advantageous to consider first the transformation of individual curve segments because attention can be restricted to the application of regularisation to the simplest class of polynomial basis conversion problem. This is an important point because Tikhonov regularisation and truncated singular value decomposition are closely related to the numerical solution of Fredholm integral equations of the first kind, and these equations, and their discretised equivalents, are known apriori to be ill-conditioned. However polynomial basis transformation problems are not associated with integral equations, and thus the linear algebraic equation that defines the basis transformation may be ill-conditioned. Tikhonov regularisation and the discrete Picard condition are introduced in Section 3, and it is shown in Section 4 that if this condition is satisfied, the equation is ill-conditioned but that the converse is not true. This leads to the generalised singular value decomposition, and this is considered in Section 5. Regularisation must, of necessity, introduce an error, and this error is quantified in Section 6. It is shown that if the discrete Picard condition is satisfied, this error is guaranteed to be small, that is, the regularised solution is a good approximation to the theoretically exact solution.
0307-904x/97/$17.00 PII s0307-904x(97)00081-4
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
An important issue in regularisation is the degree of smoothing that is imposed, because if too little smoothing is applied, the computed solution may be sensitive to small perturbations in the input data, but if too much regularisation is applied, the computed solution is too smooth, that is, some of the desired underlying solution has been removed from the computed solution. A graphical method for the determination of the degree of smoothing to be imposed is described in Section 7, and Section 8 contains an example of the conversion of a Bezier curve segment of degree 12 to its B-spline form.
2. Polynomial
basis transformation
The transformation between a Bezier curve segment its B-spline form, both of degree n, is’ JH=GF
and
(1)
where J, G E !Rmxm, H, FE !Rmmx2,m = n + 1, and rank J = rank G = m. Explicit forms for the elements of G and J are in Ref. 1. The matrices H and F store the coordinates of the vertices of the control polygons of the B-spline curve and Bezier curve, respectively. The conversion between the Bernstein and B-spline polynomial bases can be performed by the generalized Oslo algorithm, but it is noted in Ref. 3, page 152, that “. . . it is possible to determine that the generalized Oslo algorithm is stable in whether the algorithm is some cases.. . . To determine unstable for given knot vectors.. . seems like a more difficult task.” It is shown in Ref. 1 that the condition number of J increases rapidly with the degree of the polynomial, and it follows that there exist situations in which the basis conversion problem is ill-conditioned. It will be shown that when the transformation is performed numerically, these ill-conditioned situations are readily apparent, and that in some, but not all, situations, Tikhonov regularisation leads to a computationally reliable solution that is a good approximation to the theoretically exact answer. It is shown in Ref. 1 that it is necessary to rewrite equation (1) as Jh, = (GF)l
Jh, = (GF);?
where H = [h, h,] and GF = [(GF), these equations separately. It follows linear algebraic equation
(2a,b)
GF),l, and that
3. Tikhonov
regularisation
Since equation (3) is ill-conditioned, it is necessary to regularise the equation, that is, make it less sensitive to noise. Most numerical methods for the regularisation of a linear algebraic equation seek to overcome the problems of ill-conditioning by replacing it with a “nearby” wellconditioned problem whose solution approximates the desired solution and is more satisfactory than is the simple solution. The latter goal is achieved by the addition of more information about the sought solution, often that the computed solution should be smooth. These methods are called regularisation procedures, and they always include a parameter that controls the degree of smoothing that is imposed on the computed solution. For example, the regularisation parameter in truncated singular value decomposition is a discrete index that controls the number of nonzero singular values that are retained in the expression for the solution x. The solution that is obtained from a regularisation procedure must possess the following two properties: (i) It must be a good approximation to the theoretically exact solution; and (ii) it must be less sensitive to noise than is the simple (direct) solution of the equation that is being regularised. It will be shown that both properties are satisfied if it is known a prior-f that the theoretically exact solution satisfies the discrete Picard condition. The regularised solution of Tikhonov regularisation4-6 is defined by the unique solution of min {IlAx - bll’ + ~u~llx1l~l X
(4)
where (Y is the regularisation parameter and controls the degree with which the constraint on ((x1(is imposed on x = x( CY). 11(1 denotes the 2-norm. The solution x(a) of this minimisation problem is given by the solution of the equation 64% + (Y21)X(a) = ATb
(5)
This equation is most conveniently considered by using the singular value decomposition7 of A, given by USVT where U, S, V E sm”““, U and V are orthogonal matrices, and S is the diagonal matrix whose diagonal elements, sii = si, i = l..m, are the singular values of A, arranged in nonincreasing order. It follows that
solve
the general
x(cy)=v(s2+(y21)-1sc= 5
L-v.
ci
i=l S12+(Y2 Si ’
Ax=b
(3)
can be considered. A E smxm, x, b E %“‘, and A is nonsingular and ill-conditioned. Equation (2a) defines the transformation between the x-coordinates of the vertices of the control polygons of the Bezier and B-spline forms of a curve, and equation (2b) defines the transformation between the y-coordinates of the vertices of the control polygons of the curves.
652
Appl.
(6)
Math.
Modelling,
1997,
Vol. 21, October
where c = {ci} = UTb, and v, is the ith column of V. U and V are the matrices of the left and right singular vectors, respectively. The theoretically exact solution of equation (3) is
x(0) = vs-‘c
m ci = c ;vi i=, I
(7)
Tikhonov regularisation Equation
X(a)
(6) can be written
in the form
is solved. The solution
Y,(,,%‘~
= ~ i=l
is given by
2
1
&W=*=
2’
1
1+
i = l..m
(9)
” i s, 1
It is clear that
=
J. R. Winkler
is
ci + SC; 5 ~ v L s;
i
i=l
i
Since the ratios Ici + 6cil/si and Ic,l/si usually span several orders of magnitude, it is more convenient to consider the Picard ratios logic; + 6ci(/si and logIciJ/s, flog = log,,). Since it is assumed that the desired solution satisfies the discrete Picard condition, lcil + 0 as i + m. However the components lc, + 6ciI/si typically decrease to a minimum level that is determined by the noise and then increase, as shown in Figure 2. In particular, if p is the index at the minimum point of the curve in the figure, then
1
if (Y<< si
lcil B lSc,l
l/2
if a=si
lcil=lSc,I
if (YB si
lcil e 16cil
i 0
conversion:
of this equation
x(0) + 6x(O) =
(8)
I
where the filter factor g,(a)
g,(a)
for polynomial
for
i
for
i=p i >p
for
Figure 2 shows that The variation of gi( CY>with (Y/S; is shown in Figure 1. As cr increases, the filter factor g;(c~y) decreases, and the terms in equation (8) that are associated with the small singular values are removed from the expression for X((Y). It follows that if X((Y) is to be a good approximation to x(O), it is necessary that x(O) be dominated by the large singular values, which implies that the ratio lcil/si must decrease towards zero as the index i increases. Since the singular values are arranged in nonincreasing order, it is necessary that the components lcil decrease towards zero more rapidly than do the singular values. This is the discrete Picard condition.*
x(0) + 6x(O) =
+vm
(11)
m
and thus the simple solution of equation (10) is dominated by the noise and all knowledge of the exact solution is lost. However the regularisation parameter CY can be chosen such that g,(ff)
=
1 0 i
i=l..p-1 i=p..m
in which case Definition: The discrete Picard condition. The unperturbed right-hand side b of equation (3) satisfies the discrete Picard condition if for all numerically nonzero singular values si, the coefficients lcil decay towards zero, on average, more rapidly than do the singular values. In practical situations the right-hand side vector b is corrupted by noise such that the perturbed equation A(x+6x)=b+6b
(10)
This solution is a good approximation to x(O) if the decay of the ratio Icil/sj is sufficiently rapid. It is seen that the satisfaction of the discrete Picard condition guarantees that the regularised solution x(cu) + OX is a good approximation to the theoretically exact solution x(01.
Picard ratio Ii
filter factor
,
i
P
Figure
1.
The variation
of the filter
factor
g,(a)
with
a/s,
Figure 2.
Appl.
Math.
The variation
Modelling,
of the Picard
1997,
ratio logjc,+
Gc,Vs, with
Vol. 21, October
653
i.
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
The difference between truncated singular value decomposition and Tikhonov regularisation is the manner in which the singular values are removed; in the former method, they are removed in discrete steps by reducing the value of an integer, but in the latter method they are reduced towards zero gradually by increasing the value of the continuous parameter CY.The close relationship between these two regularisation methods will be shown by example in Section 8. Tikhonov regularisation and truncated singular value decomposition are based on the following heuristic:5
Heuristic. The number of oscillations in the left and right singular vectors ui and vi tends to increase with increasing i. The number of oscillations in a vector z is a measure of its spectral components; as the number of oscillations increases, the components of z exhibit a greater number of sign changes. The heuristic6 is clearly evident in Figure 3 in Hansen (1992j6, taken from an example of image restoration. It is clear that when the heuristic is satisfied, the regularised solution is smooth because it is a linear combination of nonoscillatory vectors. This oscillatory feature is satisfied when equation (3) is derived from the discretisation of a Fredholm integral equation of the first kind using orthonormal basis functions.6 However this property may not be satisfied for polynomial basis conversion problems. For example, Figure 3 shows the left singular matrix U for the transformation matrix P between the power and Bernstein polynomial bases for univariate polynomials of degree twelve, where Py = r and y and r store, respectively, the coefficients of the polynomial in the Bernstein and power bases.’ Similarly Figure 4 shows the right singular matrix V of P. Figures 5 and 6 show, respectively, the matrices U and V for the transformation matrix J in equation (1). It is clear from the figures that the heuristic is not satisfied for the transformation matrices P and J. It is interesting to note that only the left singular matrix of J has a strong pattern. Even though the
column Figure 4. matrix P.
index The right singular
matrix V for the transformation
heuristic may not be satisfied for polynomial basis conversion problems, it is shown in Section 8 that Tikhonov regularisation is effective in regularising a specific class of ill-conditioned linear algebraic equation.
4. The stability
of linear algebraic
equations
It was shown in the previous section that if the discrete Picard condition is satisfied, the solution of equation (3) is susceptible to the effects of random noise, from which it follows that the equation is ill-conditioned. This qualitative statement is quantified in this section, and it is shown that Tikhonov regularisation can only solve a particular class of ill-conditioned linear algebraic equation. The ratio relative error in x( (Y> relative error in b
U
column Figure 3. matrix P.
654
index
column
”
The left singular
Appl.
Math.
matrix
Modelling,
U for the transformation
1997,
Vol. 21, October
Figure 5. matrix J.
index The left singular
matrix U for the transformation
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
It is important to note that this measure cannot be used if the discrete Picard condition is satisfied. Consider the effective condition number of equation (10)
K,&,
IIC+ 6cll
b + 6 b) = L s,
Ils-‘(c
+ Sc)II
Since it is assumed that the discrete Picard condition is satisfied and that llbll x=-IlSbll, it follows that IIS-‘(c + 6c)ll = llS-'Sell = ISc,l/s, and Ilc + 6cll= llcll = Icll. Thus ~,,(-kb column Figure 6. ma’trix J.
index The right singular
matrix V for the transformation
is called the sensitivity function, and a formula for it is developed in Theorem 1. This ratio is important for the determination of the error magnification that may occur.
This shows that the effective condition number is dominated by the effects of noise, and it therefore provides totally misleading answers. However the measure in equation (15) is used for the derivation of the following important result. Further study of the discrete Picard condition is facilitated by postulating the following simple relationship between the magnitude of the components ci and singular values si
Theorem 1. Let A E ‘Smxm,x, b E S”, and rank A = m. The sensitivity function 5(x( a), b) for equation (5) is
Jcil =
sf,i = 1..m,
830
The parameter 0 determines ratio IciI/si with i:
S(x(a),b) =
ICll + ab) = ~sc,l
llW + a21)%cII IISCII
Ml ll(S2 + c2I)%ll
(12)
0 < 0 < 1 j IciI/si increases
(16) the dependence
of the
monotonically, i= l...m
where
SC = U% b.
O=l*(ci(/si=l,i=l...m
Proofi It follows from equation
1 < 0 * Icil/si decreases (6) that
8x((Y) =V(S2 + cr21)-1s6c The relative
error in X((Y) is given by ll(S2 + cY21)-k6cJI
IlSx((r)ll Ax(a)
= IlX(Q)II
since V is an orthogonal given by Ab = ll6bll -=llbll
=
ll(S2 + (y21)-1Scll
matrix. The relative
(13)
error in b is
IlScll (14)
llcll
Equations (13) and (14) lead directly to equation (12). This completes the proof. The effective condition number of equation (3) is obtained by setting CY= 0 in equation (12) and taking the maximum value over all perturbations SC.” This yields +&‘&b)
I
llcll
= sc~a& Sk(O), b) = y- m IIS-‘cll
(15)
monotonically,
i = 1. . . m
This model distinguishes clearly between the right-hand side vectors that do and do not satisfy the discrete Picard condition. It should be noted that although the model in equation (16) is simpler than is likely to occur in a real problem, it is adequate for the sequel. It is shown in Ref. 11 that when this model of c is substituted into equation (15), equation (3) is ill-conditioned if the discrete Picard condition is satisfied (0 > l), but that the converse is not true. This emphasises the importance of’establishing that the discrete Picard condition is satisfied and shows that Tikhonov regularisation can only solve a particular class of ill-conditioned linear algebraic equation. The satisfaction of the discrete Picard condition can be verified by plotting the terms Ici + &,I and si against i. If the coefficients Jci + &,I decay towards zero more rapidly than do the singular values si, it can be assumed that the discrete Picard condition is satisfied. However if the coefficients level off at a value that is sufficiently small, then this value is an estimate of the noise level in c. The imposition of mild statistical assumptions about the noise 8 b sheds further insight into the discrete Picard condition, and moreover, it reveals the degradation in the solution x(a) that may occur if this condition is not satisfied. The square of the sensitivity function of equa-
Appl.
Math.
Modelling,
1997,
Vol. 21, October
655
Tikhonov
regularisation
for polynomial
conversion:
J. R. Winkler
tion (5) is
into the expression
for $(cr) it follows that
s%(a),b) =
KS2 + a~I)-‘S6cl12
llCl12
l16cl12
i=l
(Sf+(Y2)2
ll(S2+ a21)- ‘Scl12 _
j=l
6cV
E
(17)
(Sf+a2)3
I I
5
i-1 (S,‘+ (Y2)3 j=l Theorem 2. The derivative S 2(x( (Y>,b) satisfies
with
dS(x(a ), b)
S(x(a),b)
da2
= -
respect
to
a2
of
which simplifies
ll~cl12
a’)’
to
Since this expression written as m-l
+ (Y~I)-~SSC
tic(Y)=
(19a, b)
m
c
i = j, it can
is zero when
6c~c;s,‘s;(
s’
-
be
si”)
c
i=l
j=i+l
(S,2+a2)3(S,~+CX2)3
+ a21)-3Sc
T3 = 6cTS(S2
+ a21)-3S6c
m-l
m
j=l
c i=j+l
+c
(19c, d) T4 = crS(S2
+
T42
where
T, = crS(S2
(Sf
11412U,T2 - TJ,) (18)
T, = 6cTS(S2
+:
6c~c;s;s;( (Sf
+
si’ -
c12)3(S;
+
ST) CT’)’
+ LU~I)-~SC
Proofi
It follows from equation (17) that it is necessary evaluate the derivative of the expression
to
On making the substitutions i + j, j + i in the second term on the right-hand side, this expression becomes
ll(S2 + (Y21)-1s6c112
ll(S2+ a2K1Scl12 and then multiply the result by (l~~~~/ll6cJl~. This completes the proof. Let 6 b be a random vector with a diagonal covariance matrix proportional to the identity matrix
If the statistical assumption about 6c in equation (21) is substituted into this expression, it follows that the expected value of $4 a 1 is m-l E{+((Y)}=E~
c i=l
E(6b 6bT) = ~~1
m
s~sj’(s~-q(c;-c~)
c j=i+l
(Sj
+
a2)3(S,!
+
Cf2)3
(20) (22)
where I E !RH”““,c2 is the variance of each component 6 b and E is the expectation operator. It follows that E{ 6c 6cT} = ~~1
of
(21)
Theorem 3. If the right-hand side vector b in equation (3) satisfies the discrete Picard condition and the perturbation 6 b satisfies equation (20), then Tikhonov regularisation reduces, on average, the effect of noise on the solution of equation (3). Proofi
It is required
that
+(a)=T,T2-T3T4<0
because this implies that the derivative of the sensitivity function is negative. Also, it follows from equation (21) that 116~11~ = me2. On substituting equations (19a)-Cd)
656
Appl.
Math.
Modelling,
1997,
Vol. 21, October
If the discrete Picard condition is satisfied, the average value of 1+4(a) is negative and Tikhonov regularisation reduces, on average, the effect of noise on the solution of equation (3). However if the discrete Picard condition is not satisfied, the term on the right-hand side of equation (22) may be positive, in which case Tikhonov regularisation degrades the solution of equation (3) because the effect of the noise 6 b is amplified. This completes the proof. Theorem 3 shows that if lcjl < lcil for j > i, the effect of Tikhonov regularisation is a reduction in the sensitivity of x(a) to noise. However the satisfaction of this condition does not guarantee that the discrete Picard condition is satisfied, and thus the error between X((Y) and the theoretically exact solution may be large, such that x(a) is of little use. It was noted in Section 1 that this paper is concerned in standard form, and it has with Tikhonov regularisation
Tikhonov regularisation been shown that this form of regularisation can only stabilise a particular class of ill-conditioned linear algebraic equation. It is shown in the next section that the generalisation of this form of regularisation enables other classes of ill-conditioned linear algebraic equation to be stabilised.
for polynomial
Theorem 4. The generalised singular value decomposition. Let the matrices A and L be defined in equations (31 and (23b), respectively. Then there exist matrices R E 8 mxm,T E %PpxPwith RTR = I,,,,TTT = I, (the identity matrices of orders m and p, respectively), and a nonsingular matrix X E Ytmxm such that A=RNX-’
5. The generalised
singular value decomposition
conversion: J. R. Winkler
L = T[M,O]X-’
(26a, b)
where
It is clear from equation (4) that Tikhonov regularisation places a constraint on Ilxll. It is possible to place constraints on the norm of other functions of x, such as the norm of its first derivative. Thus equation (4) can be written more generally as
N = diag[ r1 . . ap
1.. 11 E %I”““, (27a)
M=diag[~l..~plE%pPXP and
min{JIAx - b112+ a2~jLx~~2)
(23a)
x
where L E %Ppxm,m >p,rank(L)
=p, N(A) n N(L) = (0)
ai*+j$=l,
(27~)
i = l..p
(23b) The generalised IlLxll is a semi-norm and L is typically a discrete approximation to a derivative operator. N(X) denotes the null space of X, and it is shown in Theorem 5 that the condition that the intersection of the null spaces of A and L is the zero vector guarantees that equation (23a) has a unique minimiser. Examples of L include L, E MC” - ‘jxrn and L, E !@‘-2)xm, given by 1 0
-1
0
0 0
l-1
0 0
0 0
1
&,
0 1.00
0 00.
-1 L, =
0 01. 2
0 0 0 0
l-l
i = l..p
If L = I,, then X-’ = M-ITT and A = R(NM-‘)TT, and the generalised singular values y, reduce to the singular values of A (with pi # 0, but in reverse order. It is convenient to partition the matrices R, N, and X as
-1
0
0 0
2
-1
0
2 -1 .
-1
[
0 0 0
0 0 0
2.
0 -1 I
-1
L, and L, are, respectively, proportional to approximations of the first and second derivative operators on a uniform net. It follows that the minimiser of equation (23a) is given by the solution of the equation (ATA +
a2LTL)X(a)
=
ATb
I
NP O
N=
IO
0
0. I
0 0
-1
-1
0 0
01-l . 0
2
&,I
0
-1 0.
-1
values ‘y, are defined by
Pi
R=[R, L,=
singular
(24)
X=[X,
x,1
(28a, b, c)
where R, and X, have p columns, Np E %PxP and o = m -p. Equations (26a), (26b), (27a1, (27b), (27~1, (28a), (28b1, and (28~) enable the matrix A’ to be simplified to A’ = X,FN;R;
+ X,R;
where Nl is the pseudo-inverse the diagonal matrix
of Np, and F E %Pp*P is
It follows that the expression written as
in equation
(25) can be
and thus X(a) = (ATA + (y2LTL)-lATb
(25)
The matrix A’ = (ATA+ (Y*LTL>-’ AT is most conveniently analyzed by the generalised singular value decomposition’ of the matrix pair (A,L). This decomposition reduces to the singular value decomposition of A when L = I.
x((y)’
2
5 yi i=, ++a*
qTb -xi a,
+
E (rTb)xi i=p+l
(29)
where ri is the ith column of R. It is noted that the generalised singular values are arranged in nondecreasing order, which must be compared with the order in which
Appl. Math.
Modelling,
1997,
Vol. 21, October
657
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler
the singular values si are arranged. The second term on the right-hand side of equation (29) is independent of the regular&&on parameter ry, and the vectors xi, i =p + l..m lie in the null space of L. This follows immediately from equations (26b) and (28~) because LX=L[X,
X,1 =T[M
Y?
i = l..p
+Cr2’
‘yi
As the regularisation parameter (Y increases, the terms in equation (29) that are associated with the small generalised singular values are filtered out from the expression for x(a). It follows that the expression in equation (29) is a good approximation to the theoretically exact solution x(O) if the dominant components of x(0) are associated with the large generalised singular values yi, that is, the expression in equation (29) is dominated by the components that are associated with large values of the index i, and this implies that it is necessary that
lr’bl
-=-ai
1 IrTbl +O Pi
asi-
Yi
Since pi + 1 as i + 1, it follows approximation to x(O) if
lr,% -ro
rank[ t] where
01
and thus LX, = 0. Equation (29) shows that the solution x(a) is expressed as a linear combination of the vectors xi and that the filter factor g&a) in equation (9) is replaced by the filter factor
2
is y = 0, and this implies that the rank of the coefficient matrix of this homogeneous equation is m. Since
that
x(a)
is a good
= rank[ aAL]
(Y is an arbitrary rank1 iL]
and
nonzero
= raok{MT
rankA = rankATA
constant,
it follows that
aLTl[ :L]}
= rank(ATA + cy2LTL) = m
(31)
( = ) If rank(ATA + (Y2LTL) = m, it follows from equation (31) that the only solution of equation (30) is y = 0, and this implies that N(A) n N(L) = (0). This completes the proof. If the ratio Icil/si does not decrease towards zero as the index i increases but equation (3) is ill-conditioned because 8 = 1 in equation (16), it is necessary to select a matrix L such that the coefficients lrTb[ decay towards zero more rapidly than do the generalised singular values 7,. The geometric interpretation is that a set of linearly independent vectors xi is sought such that the dominant components of x(O) lie in the subspace of St” that is spanned by those vectors xi that are associated with the large generalised singular values ‘yi. An example of the use of the matrix L, for the regularisation of the equation that arises from the discretisation of a Fredholm integral equation of the first kind is in Ref. 8. It is concluded that the form of regularisation in equation (5) is a particular example of a more general class of regularisation methods, and that this larger class must be considered when solving ill-conditioned linear algebraic equations. The selection of the matrix L is a major research issue, and the incorrect choice yields an erroneous answer.
asi-+l 6. The regularisation
Yi
Since yi -+ 0 as i --) 1, it follows that the coefficients lr’bl must decay towards zero more rapidly than do the generalised singular values yi. This is the discrete Picard condition for equation (241, and its satisfaction is crucial for the correct use of Tikhonov regularisation in general form.
Theorem 5. Let the matrices A, L be defined (3) and (23b), respectively. Then
in equations
error
It is shown in Ref. 1 that truncated singular value decomposition introduces a bias because the regularised solution does not reduce to the theoretically exact solution in the absence of noise. It is shown in this section that a similar result is valid for Tikhonov regularisation. The total absolute error in the solution of equation (3) is x(0) -(x((Y)
+ 6x((Y))
=x(O) - (ATA + (r21)-‘AT(b + 6b) N(A) n N(L) = (0)
ts
rank(A’A + (Y’ L’L) = m
Proofi (a) If N(A) n N(L) = {O}, the only solution equation
[
658
I
A L Y=o
Appt.
and this simplifies
Modelling,
1997,
+ s b)
to
of the
(30)
Math.
= X(O) - (ADA+ n2r)-‘~T(Ax(0)
Vol. 21, October
The first term in this expression is the error, in the absence of noise, due to the regularisation procedure, and the second term includes the effects of noise and the
Tikhonov
regularisation procedure. It is seen that as 6 b -+ 0, the absolute error does not approach the zero vector, which implies that Tikhonov regularisation introduces a bias in the solution X((Y). The regularisation error e(a) is defined as the normalised value of the total absolute error in the absence of noise e(a)
=
=
regularisation
for polynomial
IA
B
I C
I
ll(I - (ADA+ (Y~I)-~A~A)x(O)II
IIx(O)II
Figure
(16) is
Pro08 See Ref. 5, Theorem
3.1. If (Y is sufficiently small compared to s1 and 8 > 1, X((Y) is guaranteed to approximate x(O), such that the regularisation error e(cr> is small, and the larger the value of 8, the better the approximation. However if the discrete Picard condition is not satisfied (0 < 0 Q l), x( LY)is a poor approximation to x(O). As the magnitude of the perturbation increases, the required value of (Y increases because a larger number of the small singular values must be removed from the expression in equation (6). It follows that 0 must increase in order to maintain a small value of e(a).
7. The L-curve The success of Tikhonov regularisation in regularising equation (3) is dependent upon the value of (Y,even if it is known that the theoretically exact solution satisfies the discrete Picard condition. If the chosen value of cr is too small, x(cw> may be corrupted by noise, but if it is too large, a significant portion of the desired solution may be removed from X((Y). It is shown in Refs. 12 and 13 that the L-curve is a good method for the calculation of the optimum value of LY. Equation (4) shows that it is desired to minimise IlAx - bll while constraining Ilx(cr>ll, from which it follows that a method for the optimum determination of a should measure both these quantities. The L-curve is a parametric plot of logllx(cy)II against logllAx(cy) - bll. It is shown in Refs. 12 and 13 that if (i) the unperturbed right-hand side vector b of equation (3) satisfies the discrete Picard condition, (ii) the perturbation 6 b is a zero mean random variable that satisfies equation (20), and (iii)
J. R. Winkler
log magnitude
IIx(0) -x(a)II IIx(O) II
Theorem 6. If the simple model for lcil in equation assumed, the regularisation error e( (Y> satisfies
conversion:
7.
A
schematic
D
E
log residual
diagram
of an L-curve
II6 bll < llbll, the sequence of points that define the L-curve assume the shape shown in Figure 7. The value of (Y increases as one progresses along the curve from A to E. In particular, the section ABC of the curve is characterised by solutions that are dominated by the effects of noise, and the portion CDE is characterised by solutions that have been smoothed too much. The optimum value of (Y is its value at the corner C because it provides the best value that minimises both IlAx - bll and Ilx(cy)II. It is shown in Refs. 12 and 13 that the curve in Figure 7 is a combination of two separate L-curves (Figure 8): (9 the L-curve when b = b, is defined by zero mean random noise that satisfies equation (20); and (ii> the L-curve when b = b, has no noise component and satisfies the discrete Picard condition. When (Y is small, the composite L-curve is dominated by the effects of b,, but as it increases, the curve is dominated by contributions from b,. The region between these two extremes defines the corner C, as shown in Figure 7.
8. Example This section contains an example of the conversion of a curve segment of degree 12 (m = 13) from the Bernstein basis to the B-spline basis. Since high-degree polynomials (up to degree 15 or 21) are used in some computer-aided design systems, l4 there is a need to obtain a computationally reliable answer to this ill-conditioned basis transformation problem. The notation of equation (3) that was used in the previous sections is replaced by the notation of
log magnitude _______ : [ (a)
(b)
I Figure
8.
log residual The
L-curves
that satisfies
the discrete
Math.
Modelling,
Appl.
for (a) only Picard
noise
and (b) only
signal
condition.
1997,
Vol. 21, October
659
Tikhonov regularisation
for polynomial
Bezier
conversion: J. R. Winkler
curve
0.28
0.27
\
-81 Figure (2a).
0.26
0.25
t=o 0.24
...
Figure 9. polygon.
...... 0.29 0.3
0.31
0.32
0.33
0.34
The Bezier curve and the convex hull of its control
equations (2a) and (2b). In order to facilitate comparison between Tikhonov regularisation and truncated singular value decomposition, the data of Example 2 in Ref. 1 is used. Figures 9 and 10, respectively, the Bezier and B-spline curves and the convex hulls of their control polygons. Although the curves are necessarily identical, their control polygons differ. The vertices of the control polygons in these figures, and Figures 13 and 16(b), are denoted by black dots. Zero mean Gaussian noise was added to the right-hand side of equations (2a) and (2b), and its variance was adjusted so that
ll(G~h II II 6(GF)l II
11 (GFh 11 = 2.12 x lo3
B-spline 0.5
II6(GF)2 II
11.
(a) loglc,I/s,
and (b) logic,+ 6c,I/s,
Figures 11 and 12 show the variation of the ratios loglcil/si and logic, + 6ciI/s, with i for equations (2a) and (2b), respectively. They show that the simple solutions of equations (2a) and (2b) is unsatisfactory. This is confirmed by Figure 13, which shows the convex hull of the control polygon of the B-spline curve when the random noise is applied to the right-hand side of the equations and no regularisation is used. Comparison of Figures 10 and 13 shows that all knowledge of the exact solution is lost and that the noise dominates the solution. Figure 14 shows the B-spline curve in this circumstance, computed from
p(t)
+
6p(t) = t(GF + 6(GF)) = tJ(H + 6H)
and comparison of Figure 9 and 14 shows that the change in the curve is negligible. Since 6 p(t) = t 6(GF) it follows that
11 Sp(t)ll,
< lltllxll S(GF>h
=
11 S(GF)lI=
= 1.07 x lo3 for 0
curve
F
4
Picard ratio
(b)
2
Figure 10.
The B-spline curve
and the convex
hull of its con-
Appl.
Math.
Modelling,
1997,
Vol. 21, October
,---,*
-8 t Figure 12.
trol polygon.
660
for equation
(a) loglc,I/sj
and (b) logic;+ Sc,I/s,
for equation G’b).
Tikhonov
B-spline
basis
regularisation
for polynomial
conversion:
J. R. Winkler
log magnitude
log magnitude
Figure 13. The convex hull of the control polygon of the B-spline curve when noise is applied to equations (2a. b).
is small. (Note that the infinity norm clidean norm is used.) For 0 < t Q 1 Ilp(t)ll,
= Il[x(t)
ywlll,
and not the Eu-
= 0.34
and this occurs at t = 0.8. It follows that a measure relative change in the curve due to the perturbation
of the is
II 6pCt) IL < 4.38 x 1O-4 IIp(t) IL which is small. Figures 13 and 14 show that it is necessary to distinguish between the integrity of the transformed curve and the integrity of the transformed convex hull. Figure 15(a) and (b), shows the discrete points that define the L-curves for equations (2a) and (2b), respec-
B-spline
(b) Figure 15. (a) The L curve for equation (2a). -, Tikhonov regularisation, * truncated singular value decomposition. (b) The L-curve for equation (2b). -, Tikhonov regularisation, . truncated singular value decomposition.
basis
tively. The results that are obtained when truncated singular value decomposition is used, reproduced from Ref. 1, are shown, and it is seen that there is excellent agreement between the L-curves for the two methods. Figure 16(a) and (b) show the B-spline curve, and the B-spline curve and convex hull of its control polygon, for (Y, = (Y*= 1.2 x 10e4, where (Y, and (Y* denote the optimum regularisation parameters, as determined from Figure 15(a) and (b), respectively. It is seen that Tikhonov regularisation is very effective in restoring the quality of the convex hull of the B-spline curve segment. Very good answers were also obtained for (Ye= CY*= 1.2 x lo-“. The magnitudes of the residual and solution at these two values of (Y are
0.301
log,,II Jh,( a) - (GFh 11= -3.841 Ly= 1.2 x lo-“:
0.28.
log,,llh,(a)ll
= -0.177
log,& Jh,( LY)- (GF)z 11= - 3.597 log,,ll h,( cr) II = -0.24?
0.26.
log,,Il Jh,(c.u) - (GF), II = -3.518 cX= 1.2x
0.24.
10-4:
t=o
log,,llh,(a)ll
log,,tI Jh,( a, 1 - (GF)2 II = - 3.496 loglollhz(a)ll
0.28
0.30
0.32
= -0.177
= -0.248
0.34
Figure 14. The B-spline curve when noise is applied to equations (2a. b).
Comparison of these results and those of Example 2 in Ref. 1 confirm that if the discrete Picard condition is
Appl.
Math.
Modelling,
1997,
Vol. 21, October
661
Tikhonov regularisation
for polynomial
conversion: J. R. Winkler Acknowledgment The author acknowledges the support of the UK Engineering and Physical Sciences Research Council (EPSRC) by an Advanced Research Fellowship, reference number B/93/AF/1703.
0.28. 0.27.
References 0.25. 0.24
0.29
0.30
0.31
0.32
0.33
1. Winkler, J. R. Polynomial basis conversion made stable by truncated singular value decomposition. Appl. Math. Mode&g, in press 2. Boehm, W., Farin, G., and Kahmann, J. A survey of curve and surface methods in CAGD. Compui. Aided Geo. Design 1984, 1, l-60 3. Lyche, T., Morken, K., and Strom, K. Conversion between Bspline bases using the generalized Oslo algorithm, in Knot Insertion and Delefion Algorithms for B-spline Curves and Surfaces. Eds. R. N. Goldman and T. Lyche, SIAM, Philadelphia, 1993 4. Hanke, M., Hansen, P. C. Regularization methods for large-scale problems. Survey Math. Ind. 1993, 3, 253-315 5. Hansen, P. C. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM J. Sci. Staf. Comput. 1990, 11, 503-518 6. Hansen, P. C. Numerical tools for analysis and solution of Fredholm integral equations of the first kind. Inverse Problems 1992,
0.34
(a)
8, 849-872 7. Golub,
G. H. and Van Loan, C. F. Matrix Computations. Johns Hopkins University Press, Baltimore, 1989 8. Hansen, P. C. The discrete Picard condition for discrete ill-posed problems. BIT 1990,30,658-672 9. Farouki, R. T. and Rajan, V. T. On the numerical condition of polynomials in Bernstein form. Comput. Aided Geo. Design 1987, 4, 191-216 10. Winkler, J. R. The sensitivity of linear algebraic equations. Proceedings
I (b) Figure 16. (a) The ELspline curve for a1 = CY~= 1.2 X 10m4. (b) The B-spline curve and convex hull of its control polygon for aI = (~z = 1.2 x 10-4.
of the Fifth SIAM
put. 1993,
satisfied,
truncated
662
Appl.
Math.
value decomposition yield very similar solutions
singular
Modelling,
1997,
and (see
Vol. 21, October
on Applied
Linear Algebra,
14, 1487-1503
14. Lachance,
Tikhonov regularisation Ref. 5, Theorem 3.2).
Conference
Snowbird, Utah, USA, June 15-18 1994, pp. 279-283 11. Winkler, J. R. The condition number of a matrix and the stability of linear algebraic equations. Department of Computer Science, The University of Sheffield, internal research report CS-97-05 12. Hansen, P. C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Reu. 1992,34, 561-580 13. Hansen, P. C. and O’Leary, D. P. The use of the L-curve in the regularization of discrete ill-posed problems, SIAM J. Sci. Com-
mial curves. Algorithms
Shrivenham, pp. 125-133
M. A. Piecewise
polynomial
approximation
of polyno-
Proceedings of the Second International Conference on for Approximation, Royal Military College of Science,
UK, July 1988, eds. J. C. Mason
and M. G. Cox,