- Email: [email protected]
Polynomial basis conversion made stable by truncated singular value decomposition
Polynomial basis conversion made stable by truncated singular value decomposition Joab R. Winkler
The transformations between the power and Bernstein polynomial bases, and the Bernstein and B-spline polynomial bases, are dejked by a linear algebraic equation whose coeficient matrix is square and nonsingular. Direct solution of the equation is not always recommended because it may be ill-conditioned. This paper considers the use of truncated singular value decomposition for the regularisation of the equation. An important consideration in this method is the deletion of the correct number of singular values of the coeficient matrix, and a method that enables this number to be determined is described. Examples of the use of truncated singular value decomposition for both transformations are presented. 0 1997 by Elsevier Science fnc.
Keywords: computer-aided value decomposition
design, power basis, Bernstein
1. Introduction The widespread use of computer-aided design (CAD) systems has led to an increased interest in the problem of high integrity data exchange between different systems. Difficulties arise because different systems use different mathematical representations of the same geometric entities, and thus the compatibility of these systems is not guaranteed. Even when the conversion between two representations is theoretically exact, probfems may occur because of the iIl-conditioned nature of the coefficient matrix of the linear algebraic equation that defines the transformation between the representations. These problems are particularly likely to manifest themselves when derived data, such as 6) the computed curve of intersection of two surfaces or (ii) a computed surface that blends two or more other surfaces, is transformed, in which case a minor change in the representation in the sending system may have a substantial effect on the representation in the receiving system. This paper considers the problem of the conversion between different representations for which computation~ly reliable answers may be difficuit to obtain. Specifically, the transformations between the power and Bernstein poiynomiat bases, and B-spline and Bernstein polynomial bases, are investigated. It is as-
basis, Bezier curve, B-spline curve, singular
sumed that the knots in the B-spline spaced.
2. Basis transformation
basis are uniformly
matrices
The transformation of a univariate polynomial of degree n between its power and Bernstein forms is defined by the equation Py=r
(1)
where P E RmXrn, y, r E R”, P is nonsingular, and m = n + 1. The vectors y and r store the coefficients of the polynomial in the two bases. Formulas for the elements of P and its inverse are in Ref. 1. A single segment of a B-spline curve of degree n on a uniform knot set can be written as p(t) =
kdt)y(t)l
= tJH, t E LO,11
(2)
where J, H, and t are of order m X m, m X 2, and 1 x m, respectively, and rank J = m. Their explicit forms are (Cohen and Riesenfeld’; Rogers and Adams,-i p. 334; Yamaguchi: pp. 327-329) t =: [t” t*- 1 1.. t1]
Address reprint requests to Dr. J. R. Winkler at the University of Sheffield, Regent Court, Department of Computer Science, 211 Portobella Street, Sheffield Sl 4DP, UK. Received
16 August
1996; accepted
18 May 1997.
Appl. Math. Modelling 1997, 21, 557-568, September 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
J={jLj}~j=o=&(~)
~(n-k)‘(-l)k-i k-j
H = [h, h, *a*h,_ Ih,17
0307-904X/97/$1 7.00 PII SO307-904X(97~0052-8
Pdynomial
basis conversion: J. R. Winkler
Since (n - k)’ = 1 for k =j . . _n, when i = 0, it follows that
The position vectors hi = [xi yilT, i = 0.. . n, store the coordinates of the vertices of the control polygon of the B-spline curve. The form of J when the knots are not uniformly spaced is considered in Ref. 5. A Bezier curve of degree n can be written in the form p(r) = [x(t) y(t)1 = tGF
(3)
where G and F are of order m X m and m X 2, respectively, and rank G = WZ.Explicit forms for G and F are (Cohen and Riesenfeld’; Rogers and Adams,3 p. 297; Yamaguchi: pp. 195-196) G = {gij];j=,,
t
0
otherwise
F=[fOff,*-fn_,f,IT The position vectors fi = [xi yi]*, i = 0.. . n, store the coordinates of the vertices of the control polygon of the Bezier curve. On equating the coefficients of t”, i = 0.. . n, in equations (2) and (3), it follows that JH=GF
(4)
The Bernstein basis is an example of the B-spline basis that is obtained for one particular knot vector. ~thou~ the conversion between these two bases can be accomplished by the generalized Oslo algorithm6 it is noted (on page 152 of Ref. 6) that “ . . . it is possible to determine that the generalized OS/Oalgorithm is stable in some cases. . . . To determine whether the algorithm is unstable for given knot vectors . . . seems like a more diflicult task.” It is shown in
this paper that some, but not all, of the situations in which the algorithm is unstable are revealed when the basis conversion is performed by truncated singular value decomposition. Furthermore, remedial action, that is, regularisation of the solution, follows immediately. The variation of logI K(P) and log,, K(J) against R, where K(X) denotes the condition number of X with respect to the 2-norm, is shown in Figure 1. It is seen that both condition numbers increase rapidly with the degree of the polynomial. Some CAD systems use polynomials of degree fifteen or twenty-one,’ and Figure 1 shows that substantial numerical problems may occur when transforming these high degree polynomials between different representations. The ill-conditioned nature of P has led to the conclusion that the transformation between the power and Bernstein bases is not to be recommended because of the amplification of errors that may occur.* Since K(P) and
558
Appl. Math. dwelling,
1997, Vol. 21, September
t
lO.log
condition
nuber
6%
polynomial
degree
Figure 1. The variationof (al log,, K(P) end (b) log,, K(J) with the polynomial degree.
K(J) may be poor measures of the stability of equations (1) and (41, respectively, it is necessary to derive criteria that establish the conditions when these equations are stable and when they are not stable. This topic is addressed in Section 3. Since the transformation matrix P(J) is ill-conditioned, small changes in dGF) may give rise to large changes in y(H), in which case standard methods cannot be used to solve equations (1) and (4% In particular equations (1) and (4) must be regularised, that is, made less sensitive to perturbations in r and GF, respectively, in order to filter out the effects of noise. This problem has been addressedgel for equation (11, but it is usually assumed that this equation is derived from the discretisation of a Fredholm integral equation of the first kind. This approach is not satisfactory for polynomial basis conversion, and it is therefore necessary to consider regularisation from an alternative aspect. This topic is addressed in Section 4, and it will be seen that there is a close connection between the results of Sections 3 and 4. An ~~rtant issue in regul~~tion is the degree of smoothing that must be applied, if too little smoothing is used, then the reguiarised solution may still be corrupted by noise, but if too much smoothing is imposed, a significant portion of the desired solution may be filtered out. A simple graphical technique for the determination of the optimum amount of smoothing is described in Section 5. Section 6 contains two examples that demonstrate the theory. One method of reducing the condition numbers of P and J is by a sequence of row and column scalings. However this procedure was not used for the work described in this paper because there does not exist a satisfactory algorithm for scaling a general matrix (Stewart,‘4 pp. 157-158).
3. The stability of linear algebraic equations The regularisation of equations (1) and (4) by truncated singular value decomposition requires the determination of the solutions that minimise, respectively, IlPy- rllz and IUH- (GF)IIF when P and J become singular because of the deletion of one or more of their singular values. The
Polynomial
subscripts 2 and F in these expressions denote the Euclidean and Frobenius norms, respectively. Since the Frobenius norm is a generalisation of the Euclidean vector norm, that is, IDH- GFII; = 11Jh, - (GF)lIl; + 11Jh, - (GF)z 11; (5) where H = [h, h2] and GF = [(GF),(GF),], it follows that IlJH- GFIIF is minimised when each of the terms on the right-hand side of equation (5) is minimised. Since (GF), and (GF), are independent, it follows that when considering the family of solutions that minimise 11~ - rll2 and IlJH- GFIIF, rank P, J < m, one can consider, without loss of generality, the least squares solutions of the general linear algebraic equation Ax=b
Jh, = (GFh
(74
Jh, = (GF)z
(7b)
k (UTb)i x(k) = VS:UTb = c ----v. si
i=l
(81
’
The notation x(k) denotes that the solution is obtained when rank A = k. (U’b), is the ith element of Urb and vi is the ith column of V. The integer k is called the truncation parameter, and
separately, that is, solve equation (6) with A = J twice: once for (x, b) = (h,, (GF),) and once for (x, b) = (h,,(GF),). Note that equation (7a) (7b) considers the transformation between the n(y)-coordinates of the vertices of the control polygons of the Bezier and B-spline forms of a curve. The theoretically exact and unique solutions of the polynomial basis conversion problems that are addressed in this paper occur when rank A = m. It follows that the development of the theory can be expressed in terms of the generic quantities A,x, b, of equation (6). However in Section 6, where two examples are considered, the symbols in equations (11, (7a) and C7b) will be used. The usual method of measuring the stability of equation (6) is by the condition number K(A) of A. This is a worst-case upper bound of the ratio r=
J. R. Winkler
matrices, and S is the diagonal matrix whose diagonal elements are called the singular values of A. These diagonal elements sli = sj, i = 1.. . m, are arranged in nonincreasing order, si 2 s2 2 ... 2 s, > 0. U and V are the matrices of left and right singular vectors, respectively. It is easily seen that the columns of U and V are the eigenvectors of AAT and ATA, respectively. When one or more of the singular values of A is set to zero (deleted), A becomes singular and there exists an infinite number of vectors x, called the least squares solutions, that minimise the residual IlAx - b/l, where 1)* 1) denotes the 2-norm. In particular, when the last m - k singular values of A are deleted, the minimum norm least squares solution of equation (6) is given by
(6)
where A E Rm xm, x, b E R”, and rank A < m. It will be shown that it is necessary to solve the equations
basis conversion:
0
o,,,:;;rk Sk,, is the k X k diagonal matrix that contains the first k singular values of A. The subscripts on the submatrices indicate their size. The solution of equation (6) when rank A = m is
x(m) =VS-‘UTb
m (UTb)i = c ----v, i=
1
si
’
3.2 The sensitivity function of a linear algebraic equation
The main result of this section is established in the following theorem. It will be used in later sections.
relative error in x relative error in b
In particular K(A) is only a function of A, and so for a given vector b, it may be an unrealistic measure of this ratio. It follows that there exist situations in which equation (6) is stable even when K(A) is large. In order to measure the stability of equation (6), it is necessary to extend the formula for the ratio r such that it becomes a function of A and b, with the consequence that it provides a true measure of the stability of equation (6). This new formula for r is called the sensitivity function.
Theorem 1-_(Winkleri6). Let A E R”‘“, x, b E R”, and rank A = m. If the last m -k singular values of A are deleted, the sensitivity function of the minimum norm least squares solution of equation (6) is given by
S(x(k),b)
= !$$
(10)
&
where c and 6c are given by
3.1 The singular ualue decomposition and least squares solutions The singular value decomposition14T’5 of A is given by USV’, where U, S,V E Rmxm. U and V are orthogonal
Appl.
c=UTb
(lla)
6c = UT8 b
(lib)
Math.
Modelling,
1997,
Vol. 21, September
559
Polynomial basis conversion: J. R. Winkler Proof-It
follows from equation
(8) that
= VS&
6x(k)
The relative
error in x(k) is given by
IISx(k)ll llS~6cll =Ax(k) = IIx(k) II IlSlcll since V is an orthogonal given by
ll6bll
(12)
matrix. The relative
IIScll
Ab=lbll=llcll
error in b is
(13)
Equations (12) and (13) lead directly to an expression for the ratio r, and equation (10) follows. Since the perturbation 6c is not known, it is advantageous to determine the lower and upper bounds of the first term on the right-hand side of equation (10) as 6c ranges over all possible values. It is seen that for k = m
min
S(x(m)
SCERm
2
b) = ?l
maxS(x(m),b)
= ?- &-s,
IIS
1
(14b) ‘cl1
= -
max
S(x(m),b)
= K(A)
t
ficiv,
(
0
(15) of c and the filter factor fi
for ik
It is seen that fi filters out the terms that correspond to the small singular values in the expression for x(k). The importance of the discrete Picard condition can be motivated by considering the following simple relationship between the magnitude of the components ci and singular values si Jcil = se, i = 1.. . m, 0 > 0
1 S(x(m),b)
c,6crRm
=
where ci is the ith element satisfies
fi=
and thus the ratio of the maximum value of S(x(m), b) to the minimum value of S(x(m), b) for a given vector c is which is equal to K(A). It follows from given by s,/s,, equations (14a) and (14b) that min
x(k)
““’
IF-‘cll
s1
GCERrn
ill-conditioned linear algebraic equation. More precisely, it will be shown that if it is known a prioti that the theoretically exact solution, that is, the solution that is specified in equation (9), satisfies the discrete Picard solution,” then the regularised solution (equation (8) with k < m) is a good approximation to the theoretically exact solution. The derivation of the discrete Picard condition requires the determination of the conditions for which a Fredholm integral equation of the first kind possesses a square integrable solution. However, since the underlying problem in polynomial basis conversion is not governed by this type of equation, it is desirable to motivate the discrete Picard condition directly from equation (6). This topic is now considered. In TSVD the smallest singular values are set equal to zero and equation (8) is rewritten as
The parameter lcil/si with i:
K (-4)
0 determines
(16)
the dependence
of the ratio
c.8cER”’
0 < 0 < 1 = Jc,\/si increases
which shows that the minimum and maximum values of the sensitivity function, taken over all vectors c and SC, are K-~(A) and K(A), respectively. Equations (14a) and (14b) show that as the dominant components of b = UC move from the directions that are defined by the eigenvectors of AA* which are associated with its large eigenvalues, to those that are defined by its small eigenvalues, the expressions in these equations decrease. In particular they attain their maximum (minimum) values when b is aligned along the eigenvector of AA* that is associated with the largest (smallest) eigenvalue of this matrix.
4. The discrete Picard condition Singular singular method algebraic be used
560
value decomposition forms the basis of truncated value decomposition (TSVD), which is a popular for the regularisation of ill-conditioned linear equations. It will be shown that TSVD can only for the regularisation of a particular class of
Appl.
Math.
Modelling,
1997,
Vol. 21, September
monotonically, i = l(l)m
0 = 1 * lcil/si = 1, i = l(l)m 1 < 8 =j lcil/si decreases
monotonically,
i = l(l)m
Since the ratio Icil/si may span several orders of magnitude, it is more convenient to consider the ratio log,Olcil/si, and this will be called the Picard ratio. Although the model in equation (16) is simpler than is likely to occur in a real problem, it is adequate for the motivation of the discrete Picard condition. Figure 2 shows three different forms of the Picard ratio, and Figure 3 shows the variation of the filter factor f, with i. It follows from equation (15) that the component of x(k) along the direction vi is equal to the product of the ordinate values of the graphs in Figures 2 and 3 at the abscissa value i =j. It is seen from Figure 262) that for 0 G 0 < 1, a large portion of the solution x(m) lies in the region of R” that is spanned by the right singular vectors that are associated with the large values of i, that is, the effects of the small singular values dominate the expres-
Polynomial
Picard
---_ ---_ -----)__ -----_____!fl_________.
i
j Figure
2.
ratio with
Three
different
the index
J. R. Winkler
sarily monotonically. However it is clear that the more rapid the decay of the ratio Icil/si, the better the approximation of x(k) to x(m) for k
ratio
--._
basis conversion:
forms
for the variation
of the
i, (a) 0 d 0 < 1, (b) tI= 1, and (c) fl>
Picard 1.
sion for x(m). This implies that the truncation of the series in equation (8) at an index k 1, the removal of the small singular values causes a relatively small error because x(m) is dominated by the large singular values. This analysis suggests that TSVD is guaranteed to yield a solution that is a good approximation to the solution in equation (9) only if the ratio Ic,l/s, decays towards zero as the index i increases. Since the singular values are arranged in non-increasing order, it follows that the absolute values of the components ci must decrease towards zero more rapidly than do the singular values si. This is the discrete Picard condition.17
Definition: The discrete Picard condition-The unperturbed right-hand side b of equation (6) satisfies the discrete Picard condition if for all numerically non-zero singular values si, the absolute values of the components ci decay towards zero, on average, more rapidly than do the singular values. It should be noted that the discrete Picard condition need only be satisfied approximately, that is, the ratio IciI/sj must decrease “on average” to zero and not neces-
filter factor 1.
A(x+iSx)=b+6b
(17)
is solved, and this problem does not satisfy the discrete Picard condition. However if b satisfies the discrete Picard condition, then k can be chosen so that b + 6b also satisfies the discrete Picard condition, from which it follows that a useful solution of equation (17) can be obtained. In the presence of noise, the components Ici + 6ciI/si typically decrease gradually until they reach a minimum at a level determined by the errors, and then increase. (This will be shown in the examples in Section 6.) It follows that x(m) is likely to be dominated by the noise corrupted terms that are associated with the small singular values. These noise corrupted terms must be filtered out if a useful solution to equation (17) is to be obtained. The correct choice of k allows the ratio filci + 6ciI/si to satisfy the discrete Picard condition, that is, the filter factor can be chosen so that the components of x(k) along the directions vi, i = 1.. . m, that are associated with the small singular values are zero. It is shown in Ref. 18 that when the discrete Picard condition is satisfied, equation (6) is ill-conditioned, but the converse is not true because of the situation that prevails when 0 = 1. This result, which is obtained by substituting equation (16) into equation (14b), emphasises the importance of establishing that the discrete Picard condition is satisfied and clarifies the statement that TSVD can only solve a particular class of ill-conditioned linear algebraic equation. 4.1 Statistical analysis The discussion of the discrete Picard condition above was adequate to suggest the importance of this condition, but more analysis must be performed in order to assert that it must be satisfied for a useful regularised solution to exist. This section and the subsequent section consider the discrete Picard condition in more detail. Let 6 b be a random vector with a diagonal covariance matrix proportional to the identity matrix
E{6b6bT}= (~‘1
(18)
where I E R”‘“, (T’ is the variance of each component of 6 b and E is the expectation operator. It follows from equation (lib) that i
j Figure
3.
The variation
of the filter
factor
with
the index
i.
E{ScGcT} = (T21
Appt. Math.
Modelling,
1997,
(19)
Vol. 21,
September
561
Polynomial
basis
conversion:
J. R. Winkler
Theorem 2-If the right-hand side vector b of equation (6) satisfies the discrete Picard condition and the noise vector 6 b satisfies equation (18), then TSVD reduces, on average, the effect of noise on the solution of equation (6).
The first term is identically becomes
k:l<-c;~$ j=l
Proof-If it is required (6), it is necessary that S(x(k),
that TSVD
regularise
b) - Stx(k - l), b)
IL$W =--IISCII
llcll
IlSX_ 1Wl
Ikll
Ils:cll
IISCII
IlSd- 14
be greater than zero. The term on the right-hand this expression simplifies to llcll
lIs:6cll
IIS~CII IIs:_ 141
and it is therefore
necessary
= IIs:~cl1211sx-
This expression
side of
I
that
1c11*- lls:cl1211s~_,6cl12
'j
1
(20)
(21)
It is clear that if the discrete Picard condition is satisfied, the expression in equation (21) is necessarily positive and TSVD reduces, on average, the sensitivity function of the minimum norm least squares solution of equation (6). However if the discrete Picard condition is not satisfied, the expression in equation (21) may be positive or negative. 4.2 The regularisation error There are two errors that must be considered when analysing TSVD. They are (i) the error in x(k) that arises from the error in b, .and (ii) the regularisation error that arises because x(m) # x(k), k # m, in the absence of noise. These two errors must not be confused; the first considers the effect of noise in b on the solution x(k), and the second is a measure of the bias that is introduced in x(k) in the absence of noise by the smoothing properties of TSVD. It has already been shown that if the discrete Picard condition is satisfied, then TSVD is effective in reducing the effect of noise on the minimum norm least squares solution x(k). The second type of error is now considered. The total error in this solution of equation (6) is
IIs:_ ICI1 - IIS%_ ,6cll IIS~CII
II WI [
5(k)
equation
zero, and thus equation
> 0
becomes
x(m) - (x(k) + 6x(k)) If the statistical assumption about 6c in equation (19) is substituted into this expression, the expected value of 5(k) is
Using equation written as
=x(m)
-At(b
+ 6b)
(6), it follows that this expression
(22) can be
x(m) - A+(Ax(m) + 6 b) = (I - A+A)x(m) - A%b
(23) E{c$(k)I
= a* (20)
The expression as
inside the square brackets
can be written
The bias that is introduced in x(k) by TSVD follows immediately because the expression in equation (23) does not approach the zero vector as S b + 0. The term (I A+A)x(m) is the error due to the regularisation procedure, and the second term is the error due to noise. The regularisation error e(k) is defined as e(k) = 11x(m) - x(k)11 IIx(m>II when 6 b = 0. It follows from equations
(22) and (23) that
e(k) = I~(I- A+A)x(m)ll II x(m) and this simplifies
II
The analysis of this expression is most easily performed assuming the simple model for (cj in equation (16).
to
by
Theorem 3-If the simple model for lcil in equation (16) is assumed, the regularisation error of the minimum norm
562
Appl.
Math.
Modelling,
1997,
Vol. 21, September
Polynomial
least squares solution
of equation
(6) satisfies
log
Proof-See (Hansen,” theorem 3.1). This theorem shows that if 8 > 1 and ski 1 is sufficiently small compared to sr, then x(k) is guaranteed to approximate x(m). The truncation parameter k is usually determined such that the larger the error in b,the smaller the value of k, that is, more singular values must be deleted. It follows that if there is an error in b, 0 must increase to offset the smaller value of k, so that a small regularisation error is retained.
discrete
magnitude
6.
-I 4.
A schematic
log diagram
I The
equation Picard
log
L-curves (18)
for (a) zero and
(b) only
mean signal
residual
random that
noise
that
satisfies
the
condition.
condition, and the other is the L-curve when b represents zero mean random noise that satisfies equation (18). These two separate L-curves are shown in Figure 5. The graph in Figure 4 is obtained by adding the ordinate values of the graphs in Figure 5 at each abscissa value. When very few singular values are deleted (segment AB in Figure 4) the L-curve is almost horizontal, but the gradient decreases rapidly as more singular values are deleted until the corner C is reached. In the vertical portion BC, x(k) is dominated by the effects of errors in b, magnified by the division by the small singular values. The deletion of these singular values causes a large change in log,,llx(k)ll, but the change in log,,l~(k) - bll is small. The gradient of the curve then increases rapidly, and the L-curve is approximately horizontal (segment CD in Figure 4). The change in logrollx(k)ll in this segment is small, but the change in log,,l!Ax(k) - blI is large. The gradient of the L-curve then decreases as more singular values are deleted, corresponding to the segment DE in Figure 4. The optimum value of k is its value at the corner C, because this provides the correct balance in keeping both Ilx(k)ll and I!Axtk) - bll small; if fewer singular values were to be deleted, x(k) would be dominated by the effects of noise, but if more singular values were to be deleted, x(k) would be smoothed too much. The log-log scale is used for two reasons: (i) The singular values span several orders of magnitude, and (ii) logarithmic axes emphasise portions of the L-curve where the variation in one variable is small compared to the variation in the other variable, and this has the effect of emphasising the corner.
An important issue that must be considered when using TSVD is the choice of the truncation parameter k; if it is too large, then x(k) may still be corrupted by noise, but if it is too small, a significant portion of the desired solution may be removed from x(k). This section describes the L-curve method,‘9,20 which is a simple graphical technique for the determination of the optimum amount of smoothing. Other methods,‘9*20 including the discrepancy principle and generalised cross-validation, exist, but they suffer from disadvantages. The motivation for the regularisation of equation (6) is to limit Ilx(k)ll, and thus a criterion for the optimum choice of k must include information on I!dk)ll and ILAx - bll. The L-curve is a parametric plot of log,,lMk1ll against log,,l~(k) - bll. It is shown in Refs. 19 and 20 that if (i) the unperturbed right-hand side vector b of equation (6) satisfies the discrete Picard condition, (ii) the perturbation S b is a zero mean random vector that satisfies equation (181, and (iii) llsbll +K llbll, then the L-curve assumes the shape shown in Figure 4. It consists of a sequence of discrete points that may be joined up to form a continuous curve. This curve is a combination of two separate L-curves, one of which is the L-curve when b is not corrupted by noise and satisfies the discrete Picard
Figure
5.
satisfies
5. The L-curve
log
J. R. Winkler
magnitude
I Figure
basis conversion:
Examples
This section contains two examples that demonstrate the theory of the previous sections. Example 1 considers the conversion from the power basis to the Bernstein basis, and Example 2 considers the conversion from the Bernstein basis to the B-spline basis. The notation of equation (6) that was used in the previous sections is replaced by the notation of equation (1) in Example 1 and by the notation of equations (7a) and (7b) in Example 2.
residual
of an L-curve.
Appl.
Math.
Modelling,
1997,
Vol. 21, September
563
Polynomial
basis conversion:
J. R. Winkle/
Example 1. Let the vector r in equation (1) store the coefficients of a univariate polynomial when it is expressed in the power basis, and it is required to compute the vector y, which stores the coefficients of the polynomial when it is expressed in the Bernstein basis. The degree of the polynomial is 12, which implies that m = 13. The vector r was chosen to satisfy the discrete Picard condition by specifying 0 = 2. Zero mean Gaussian noise was added to r, and its variance was such that Ilr II/11 6rll= 105.Figures 6(a) and 6(b) show the variation of the ratios log,,lcil/si and log,,lci + 6cil/si with i. It is seen that the ratio Icil/si decreases monotonically with i, but that the noise causes the ratio lci + &,l/s, to increase, albeit not monotonically, for i > 6. It is clear that the direct solution of equation (1) with no regularisation would yield an unsatisfactory answer because of the significant effect of the noise. Figzue 7(a) shows the polynomial that is being transformed, in the range 0 < t < 1, and Figure 7(b) shows the transformed polynomial when equation (1) is solved directly and r is corrupted by noise. The curves in these two figures should be similar, but it is clear that the noise in r severely degrades the solution of equation (1). Figure 8 shows the discrete points that define the L-curve, and it is seen that the corner is rounded, i.e., not well defined, thus emphasising the need for a rigorous definition of the corner of the L-curve. The figure suggests that k = 6,7,8, or 9 because these indices define the points that are closest to the corner of the L-curve. Figzue 6(b) suggests that the summation in equation (8) should be truncated at k = 6, corresponding to the deletion of seven singular values. Figures 9(a), 9(b), and 9(c) show the (original) power basis polynomial and the (transformed) Bernstein basis polynomial for k = 5, 6, and 7, respectively. It is clear that the transformed polynomials for k = 5, 6 are similar and that they differ from the transformed polynomial that is obtained for k = 7. It is seen that the optimum value of k is 6, but k = 5 also yields an acceptable solution. The transformed polynomial is improvement in the clear-compare Figures 7(b) and 9(b).
Picard
Exact
solution
(a) Noise
corrupted
solution
(4 Figure 7.
(a) The
power
basis
to the
formed
polynomial
polynomial
whose
Bernstein
basis
when
r is corrupted
transformation is desired.
from
(b) The
the
trans-
by noise.
Figure 10 shows a sequence of L-curves for four different values of the ratio of the signal power to the noise power for a given vector r that satisfies the discrete Picard condition. As this ratio decreases, the corner of the L-curve becomes more pronounced and the optimum value of the truncation parameter decreases. This is because the portion of the ratio jci + 6cil/si that decreases terminates at a lower value of i, since the noise 6ci corrupts a larger number of the components ci, i = 13,12,. . . .
ratio
log magnitude
':::I'"'"L;log
1,
i
4.99
4.96 Figure 6.
The
variation
log,,lc,+Gcj~sj
with
564
Math.
Appl.
of the
ratios
(a)
loglo~ci~si
and
(b)
4.91
Figure 8.
i.
Modelling,
1997,
Vol. 21, September
I
-0
.
k=4
I
The L-curve.
residual .
k:2
.
Polynomial basis conversion: J. R. Winkler
k=5
log magnitude 7.06.5, (c)
6.0 5.5, .:I
-.. I
0 Figure
10.
(a)
(b) ) z
2
4
The L-curves
to (a) 30, (b) 10”.
t 6
8
for (signal
log
residual
power)/(noise
power)
equal
(cl 1Oa. and fd) 10’.
(4 k=6 Bezier
curve
0.29 0.2
u-6
1
-100 0.28
0.27
0.26
0.25
t=o 0.24 Figure
0.29 11.
0.3
The Bezier
0.31
curve
0.32
and the convex
0.33
0.34
hull of its control
polygon.
Figure
9.
(a) The
exact
for k= 5. (b) The exact k=6.
fc) The
exact
solution,
solution,
solution,
-, -,
-,
and TSVD
and TSVD and
TSVD
solution,
solution, solution,
----,
B-sgline curve
- - - -, for ----,
0.5
for
f
k=7. 0.4.
Example2. Consider the transformation of a Bezier curve of degree 12 (m = 13) to a B-spline curve. The matrix F in equation (4) stores the coordinates of the vertices of the control polygon of the Bezier curve, and it is required to compute the coordinates of the vertices of the control polygon of the equivalent B-spline curve. These coordinates are stored in the matrix H. Figures 11 and 12 show, respectively, the Bezier and B-spline curves and the convex hulls of their control polygons. The - in these (and in Figures 14 and 17(b))
Figure
The
12.
B-spline
curve
and the convex
hull of its con-
trol polygon.
Awl.
Math.
Modelling,
1997,
Vol. 21, September
565
Polynomial
basis conversion:
J. R. Winkler
denote the vertices of the control polygons. It is seen that although the curves are necessarily identically equal, the convex hulls are different. In particular, the convex hull of the control polygon of the B-spline curve is substantially larger than is the convex hull of the control polygon of the Bezier curve.21 Zero mean Gaussian noise was added to the right-hand side of equations (7a) and Ub), and its variance was adjusted so that
II II = 2.12 x 103 11 i?(GF)lll
II(GF)~II 116(GF)z II
B-spline
basis
= 1.07 x 103
Figures 13(a) and 13(b) show the variation of the ratios log,,lcjl/si and log,,lc, + tk,l/s, with i for equations (7a) and (7b), respectively. These figures show that the direct (simple) solution of equations (7a) and (7b) will yield unsatisfactory answers, because the largest values of the terms lci + 6cil/si are Figure 14. The convex hull of the control B-spline curve in the presence of noise.
polygon of the
4
Picard ratio
W
a
-6 * -8A Figure 13.
(A) The ratios (a) for equation (7a).
lOgloICiVS/
and (b) log,olci+acivsi
approximately equal to 16cil/si. This is confirmed in Figure 14, which shows the computed convex hull of the control polygon of the B-spline curve when no regularisation is applied to the equations. The cluster of points near the origin of this diagram are some of the vertices of the control polygon of the curve. Comparison of Figures 12 and 14 shows that the effect of the noise on the convex hull is substantial. Figure 15 shows the B-spline curve in this circumstance, computed from p(t) + Sp(t) = tJ(H+ 8H),and a comparison of Figures 11 and 15 shows that the change in the curve is negligible. Although the fidelity of the curve has been retained in the transformation, the gross magnification and distortion of the convex hull implies that the transformation is unsatisfactory. Figures 16(a) and 16(b) show the discrete points that define the L-curves for equations (7a) and Ob), respec-
4-
Picard ratio a,
B-spline
0.301
basis
0.26
0.24~
-8.
t=0
B Figure 13. (B) The ratios (a) for equation (7b).
566
Appl.
lOg,,ICiVS;
Math. Modelling,
0.28
and (b) log,cIci+aciVs; Figure 15.
1997, Vol. 21, September
0.30
0.32
0.34
The B-spline curve in the presence of noise.
Polynomial basis conversion: J. R. Winkler log
magnitude
k=6
J.Or l
2.0
1.0 0.0
0.29
0
0.28
a l
-6.0
-4W
46.0
l
0.
&og
0.21’
residual
0.26,
-1.01 (a)
0.25.
log
magnitude
t=o 0.24’
0.29
0.3
0.31
0.32
0.33
0.34
W
03 Figure 16.
(a) The L-curve for equation (7a). (b) The L-curve for equation (7b).
0.3
tively. Detailed examination of the curves showed that the values of k nearest the corner of the L-curve in Figure 16(u) are 5, 6, 7, and 8, and that the values of k nearest the corner of the L-curve in Figure 16(b) are 4, 5, 6, 7, and 8. Figures 17(a) and 17(b) show the B-spline curve, and the B-spline curve and convex hull of its control polygon, for k = 6. It is seen that very good answers are obtained for this value of k. Very similar answers were obtained for k = 5. The solutions for k = 7, 8 were underregularised. In general equations (7a) and 0%) may require different amounts of regularisation, that is, different values of k, or one of the equations may require regularisation and the other one may not. It is for this reason that equation (4) cannot be treated as a single equation with a multiple right-hand side. Rather, it must be decomposed into the form given in equations (7a) and G’b), and each of these equations must be solved individually. The error measure in the L-curve in Example 1 is a function of the transformation matrix P and coefficients r of the polynomial when it is expressed in the power basis. It may be desirable to measure the error by a property that is intrinsic to a curve because (i) this has more geometric meaning, and (ii) it is well known that small changes in the coefficients of a polynomial that is expressed in the power basis may alter it substantially-the Wilkinson polynomial is a common example. It follows that a small residual may not guarantee that the polynomial in the receiving CAD system is a good approximation to the polynomial in the sending CAD system. A similar situation arises for the problem in Example 2; the residual is a function of the transformation matrix J and the coordinates of vertices of the control polygon of the Bezier curve, but it may be advantageous to define the error by the proximity of the curve to its convex hull.
’
0.1
0.3
0.2
0.4
0.5
(b) Figure 17. (a) The B-spline curve for k=6. (b) The B-spline curve and convex hull of its control polygon for k=6.
Acknowledgment The author acknowledges the support of the UK EPSRC by an Advanced Research Fellowship, reference number B/93/AF/1703.
References 1. Farouki, R. T. and Rajan, V. T. On the numerical condition of polynomials in Bernstein form. Comput. Aided Geo. Design 1987, 4, 191-216 2. Cohen, E. and Riesenfeld, R. F. General matrix representations for Bezier and B-spline curves. Compuf. Indust. 1982, 3, 9-1.5 3. Rogers, D. F. and Adams, J. A. Mathematical Elements for Computer Graphics, McGraw-Hill, New York, 1990 4. Yamaguchi, F. Curves and Surfaces in Computer Aided Geomettic Design, Springer-Verlag, Berlin, 1988 5. Grabowski, H. and Li, X. Coefficient formula and matrix of nonuniform B-spline functions. Comput. Aided Design 1992, 24, 637-642 6. Lyche, T., Morken, K. and Strom, K. Conversion between B-spline bases using the generalized Oslo algorithm. Knot Insertion and
Appl.
Math.
Modelling,
1997,
Vol. 21, September
567
Polynomial
7.
8.
9. 10. 11.
12.
13.
basis conversion:
J. R. Winkler
Deletion Algorithms for B-Spline Curues and Surfaces, eds. R. N. Goldman and T. Lyche, SIAM, 1993, pp. 135-153 Lachance, M. A. Piecewise polynomiai approximation of poiynomiai curves. 2nd inte~ational Conference on Algorithms for Approximation, Royaf Military College of Science, Shrivenham, UK, July 1988, eds. J. C. Mason and M. G. Cox, pp. 125133 Daniel M. and Daubisse, J. C. The numerical problem of using Bezier curves and surfaces in the power basis. Comput. Aided Geo. Design 1989, 6, 121-128 Hanke, M. and Hansen, P. C. Regularization methods for largescale problems. Surti. Mark. Ind. 1993, 3, 253-315 Hansen, P. C. The truncated SVD as a method for regularization. BIT 1987, 27, 534-553 Hansen, P. C. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM J. Sci. Stat. Comput. 1990, 11, 503-518 Hansen, P. C. Numerical tools for analysis and solution of Fredhoim integral equations of the first kind. InrJerse Problems 1992. 8, 849-872 Varah, J. M. A practical examination of some numerical methods for linear discrete ill-posed problems. SIAM ReL). 1979, 21, 100-111
568
Appt.
Math. Modelling,
1997, Vol. 21, September
14. Stewart, G. W. Introduction to Matrix Computations. Academic Press, New York, 1973 15. Golub, G. H. and Van Loan, C. F. auto Computu?ions. Johns Hopkins University Press, Baltimore, MD, 1989. 16. Winkler, J. R. The sensitivity of linear algebraic SIAM Conference on Applied Linear Algebra, USA, June 15-18 1994, pp. 279-283 17. Hansen, P. C. The discrete Picard problems. RIT 1990,30,658-672
condition
equations, Snowbird,
for discrete
Fifth Utah,
ill-posed
18. Winkler, J. R. The condition number of a matrix and the stability of linear algebraic equations. Department of Computer Science, The University of Sheffield, memoranda CS-97-05 19. Hansen, P. C. Analysis of discrete ill-posed the L-curve. SIAM Reu. 1992, 34, 561-580
problems
by means of
20. Hansen, P. C. and O’Leary, D. P. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM .I. Sci. Comput. 1993.14, 1487-1503 21. Goult, R. J. Standards for curves and surface data exchange, Fourth IMA Conference on The Mathematics of Surfaces, Bath, England, September 1990, ed. A. Bowyer, 1990, pp. 1-15
The transformations between the power and Bernstein polynomial bases, and the Bernstein and B-spline polynomial bases, are dejked by a linear algebraic equation whose coeficient matrix is square and nonsingular. Direct solution of the equation is not always recommended because it may be ill-conditioned. This paper considers the use of truncated singular value decomposition for the regularisation of the equation. An important consideration in this method is the deletion of the correct number of singular values of the coeficient matrix, and a method that enables this number to be determined is described. Examples of the use of truncated singular value decomposition for both transformations are presented. 0 1997 by Elsevier Science fnc.
Keywords: computer-aided value decomposition
design, power basis, Bernstein
1. Introduction The widespread use of computer-aided design (CAD) systems has led to an increased interest in the problem of high integrity data exchange between different systems. Difficulties arise because different systems use different mathematical representations of the same geometric entities, and thus the compatibility of these systems is not guaranteed. Even when the conversion between two representations is theoretically exact, probfems may occur because of the iIl-conditioned nature of the coefficient matrix of the linear algebraic equation that defines the transformation between the representations. These problems are particularly likely to manifest themselves when derived data, such as 6) the computed curve of intersection of two surfaces or (ii) a computed surface that blends two or more other surfaces, is transformed, in which case a minor change in the representation in the sending system may have a substantial effect on the representation in the receiving system. This paper considers the problem of the conversion between different representations for which computation~ly reliable answers may be difficuit to obtain. Specifically, the transformations between the power and Bernstein poiynomiat bases, and B-spline and Bernstein polynomial bases, are investigated. It is as-
basis, Bezier curve, B-spline curve, singular
sumed that the knots in the B-spline spaced.
2. Basis transformation
basis are uniformly
matrices
The transformation of a univariate polynomial of degree n between its power and Bernstein forms is defined by the equation Py=r
(1)
where P E RmXrn, y, r E R”, P is nonsingular, and m = n + 1. The vectors y and r store the coefficients of the polynomial in the two bases. Formulas for the elements of P and its inverse are in Ref. 1. A single segment of a B-spline curve of degree n on a uniform knot set can be written as p(t) =
kdt)y(t)l
= tJH, t E LO,11
(2)
where J, H, and t are of order m X m, m X 2, and 1 x m, respectively, and rank J = m. Their explicit forms are (Cohen and Riesenfeld’; Rogers and Adams,-i p. 334; Yamaguchi: pp. 327-329) t =: [t” t*- 1 1.. t1]
Address reprint requests to Dr. J. R. Winkler at the University of Sheffield, Regent Court, Department of Computer Science, 211 Portobella Street, Sheffield Sl 4DP, UK. Received
16 August
1996; accepted
18 May 1997.
Appl. Math. Modelling 1997, 21, 557-568, September 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
J={jLj}~j=o=&(~)
~(n-k)‘(-l)k-i k-j
H = [h, h, *a*h,_ Ih,17
0307-904X/97/$1 7.00 PII SO307-904X(97~0052-8
Pdynomial
basis conversion: J. R. Winkler
Since (n - k)’ = 1 for k =j . . _n, when i = 0, it follows that
The position vectors hi = [xi yilT, i = 0.. . n, store the coordinates of the vertices of the control polygon of the B-spline curve. The form of J when the knots are not uniformly spaced is considered in Ref. 5. A Bezier curve of degree n can be written in the form p(r) = [x(t) y(t)1 = tGF
(3)
where G and F are of order m X m and m X 2, respectively, and rank G = WZ.Explicit forms for G and F are (Cohen and Riesenfeld’; Rogers and Adams,3 p. 297; Yamaguchi: pp. 195-196) G = {gij];j=,,
t
0
otherwise
F=[fOff,*-fn_,f,IT The position vectors fi = [xi yi]*, i = 0.. . n, store the coordinates of the vertices of the control polygon of the Bezier curve. On equating the coefficients of t”, i = 0.. . n, in equations (2) and (3), it follows that JH=GF
(4)
The Bernstein basis is an example of the B-spline basis that is obtained for one particular knot vector. ~thou~ the conversion between these two bases can be accomplished by the generalized Oslo algorithm6 it is noted (on page 152 of Ref. 6) that “ . . . it is possible to determine that the generalized OS/Oalgorithm is stable in some cases. . . . To determine whether the algorithm is unstable for given knot vectors . . . seems like a more diflicult task.” It is shown in
this paper that some, but not all, of the situations in which the algorithm is unstable are revealed when the basis conversion is performed by truncated singular value decomposition. Furthermore, remedial action, that is, regularisation of the solution, follows immediately. The variation of logI K(P) and log,, K(J) against R, where K(X) denotes the condition number of X with respect to the 2-norm, is shown in Figure 1. It is seen that both condition numbers increase rapidly with the degree of the polynomial. Some CAD systems use polynomials of degree fifteen or twenty-one,’ and Figure 1 shows that substantial numerical problems may occur when transforming these high degree polynomials between different representations. The ill-conditioned nature of P has led to the conclusion that the transformation between the power and Bernstein bases is not to be recommended because of the amplification of errors that may occur.* Since K(P) and
558
Appl. Math. dwelling,
1997, Vol. 21, September
t
lO.log
condition
nuber
6%
polynomial
degree
Figure 1. The variationof (al log,, K(P) end (b) log,, K(J) with the polynomial degree.
K(J) may be poor measures of the stability of equations (1) and (41, respectively, it is necessary to derive criteria that establish the conditions when these equations are stable and when they are not stable. This topic is addressed in Section 3. Since the transformation matrix P(J) is ill-conditioned, small changes in dGF) may give rise to large changes in y(H), in which case standard methods cannot be used to solve equations (1) and (4% In particular equations (1) and (4) must be regularised, that is, made less sensitive to perturbations in r and GF, respectively, in order to filter out the effects of noise. This problem has been addressedgel for equation (11, but it is usually assumed that this equation is derived from the discretisation of a Fredholm integral equation of the first kind. This approach is not satisfactory for polynomial basis conversion, and it is therefore necessary to consider regularisation from an alternative aspect. This topic is addressed in Section 4, and it will be seen that there is a close connection between the results of Sections 3 and 4. An ~~rtant issue in regul~~tion is the degree of smoothing that must be applied, if too little smoothing is used, then the reguiarised solution may still be corrupted by noise, but if too much smoothing is imposed, a significant portion of the desired solution may be filtered out. A simple graphical technique for the determination of the optimum amount of smoothing is described in Section 5. Section 6 contains two examples that demonstrate the theory. One method of reducing the condition numbers of P and J is by a sequence of row and column scalings. However this procedure was not used for the work described in this paper because there does not exist a satisfactory algorithm for scaling a general matrix (Stewart,‘4 pp. 157-158).
3. The stability of linear algebraic equations The regularisation of equations (1) and (4) by truncated singular value decomposition requires the determination of the solutions that minimise, respectively, IlPy- rllz and IUH- (GF)IIF when P and J become singular because of the deletion of one or more of their singular values. The
Polynomial
subscripts 2 and F in these expressions denote the Euclidean and Frobenius norms, respectively. Since the Frobenius norm is a generalisation of the Euclidean vector norm, that is, IDH- GFII; = 11Jh, - (GF)lIl; + 11Jh, - (GF)z 11; (5) where H = [h, h2] and GF = [(GF),(GF),], it follows that IlJH- GFIIF is minimised when each of the terms on the right-hand side of equation (5) is minimised. Since (GF), and (GF), are independent, it follows that when considering the family of solutions that minimise 11~ - rll2 and IlJH- GFIIF, rank P, J < m, one can consider, without loss of generality, the least squares solutions of the general linear algebraic equation Ax=b
Jh, = (GFh
(74
Jh, = (GF)z
(7b)
k (UTb)i x(k) = VS:UTb = c ----v. si
i=l
(81
’
The notation x(k) denotes that the solution is obtained when rank A = k. (U’b), is the ith element of Urb and vi is the ith column of V. The integer k is called the truncation parameter, and
separately, that is, solve equation (6) with A = J twice: once for (x, b) = (h,, (GF),) and once for (x, b) = (h,,(GF),). Note that equation (7a) (7b) considers the transformation between the n(y)-coordinates of the vertices of the control polygons of the Bezier and B-spline forms of a curve. The theoretically exact and unique solutions of the polynomial basis conversion problems that are addressed in this paper occur when rank A = m. It follows that the development of the theory can be expressed in terms of the generic quantities A,x, b, of equation (6). However in Section 6, where two examples are considered, the symbols in equations (11, (7a) and C7b) will be used. The usual method of measuring the stability of equation (6) is by the condition number K(A) of A. This is a worst-case upper bound of the ratio r=
J. R. Winkler
matrices, and S is the diagonal matrix whose diagonal elements are called the singular values of A. These diagonal elements sli = sj, i = 1.. . m, are arranged in nonincreasing order, si 2 s2 2 ... 2 s, > 0. U and V are the matrices of left and right singular vectors, respectively. It is easily seen that the columns of U and V are the eigenvectors of AAT and ATA, respectively. When one or more of the singular values of A is set to zero (deleted), A becomes singular and there exists an infinite number of vectors x, called the least squares solutions, that minimise the residual IlAx - b/l, where 1)* 1) denotes the 2-norm. In particular, when the last m - k singular values of A are deleted, the minimum norm least squares solution of equation (6) is given by
(6)
where A E Rm xm, x, b E R”, and rank A < m. It will be shown that it is necessary to solve the equations
basis conversion:
0
o,,,:;;rk Sk,, is the k X k diagonal matrix that contains the first k singular values of A. The subscripts on the submatrices indicate their size. The solution of equation (6) when rank A = m is
x(m) =VS-‘UTb
m (UTb)i = c ----v, i=
1
si
’
3.2 The sensitivity function of a linear algebraic equation
The main result of this section is established in the following theorem. It will be used in later sections.
relative error in x relative error in b
In particular K(A) is only a function of A, and so for a given vector b, it may be an unrealistic measure of this ratio. It follows that there exist situations in which equation (6) is stable even when K(A) is large. In order to measure the stability of equation (6), it is necessary to extend the formula for the ratio r such that it becomes a function of A and b, with the consequence that it provides a true measure of the stability of equation (6). This new formula for r is called the sensitivity function.
Theorem 1-_(Winkleri6). Let A E R”‘“, x, b E R”, and rank A = m. If the last m -k singular values of A are deleted, the sensitivity function of the minimum norm least squares solution of equation (6) is given by
S(x(k),b)
= !$$
(10)
&
where c and 6c are given by
3.1 The singular ualue decomposition and least squares solutions The singular value decomposition14T’5 of A is given by USV’, where U, S,V E Rmxm. U and V are orthogonal
Appl.
c=UTb
(lla)
6c = UT8 b
(lib)
Math.
Modelling,
1997,
Vol. 21, September
559
Polynomial basis conversion: J. R. Winkler Proof-It
follows from equation
(8) that
= VS&
6x(k)
The relative
error in x(k) is given by
IISx(k)ll llS~6cll =Ax(k) = IIx(k) II IlSlcll since V is an orthogonal given by
ll6bll
(12)
matrix. The relative
IIScll
Ab=lbll=llcll
error in b is
(13)
Equations (12) and (13) lead directly to an expression for the ratio r, and equation (10) follows. Since the perturbation 6c is not known, it is advantageous to determine the lower and upper bounds of the first term on the right-hand side of equation (10) as 6c ranges over all possible values. It is seen that for k = m
min
S(x(m)
SCERm
2
b) = ?l
maxS(x(m),b)
= ?- &-s,
IIS
1
(14b) ‘cl1
= -
max
S(x(m),b)
= K(A)
t
ficiv,
(
0
(15) of c and the filter factor fi
for i
It is seen that fi filters out the terms that correspond to the small singular values in the expression for x(k). The importance of the discrete Picard condition can be motivated by considering the following simple relationship between the magnitude of the components ci and singular values si Jcil = se, i = 1.. . m, 0 > 0
1 S(x(m),b)
c,6crRm
=
where ci is the ith element satisfies
fi=
and thus the ratio of the maximum value of S(x(m), b) to the minimum value of S(x(m), b) for a given vector c is which is equal to K(A). It follows from given by s,/s,, equations (14a) and (14b) that min
x(k)
““’
IF-‘cll
s1
GCERrn
ill-conditioned linear algebraic equation. More precisely, it will be shown that if it is known a prioti that the theoretically exact solution, that is, the solution that is specified in equation (9), satisfies the discrete Picard solution,” then the regularised solution (equation (8) with k < m) is a good approximation to the theoretically exact solution. The derivation of the discrete Picard condition requires the determination of the conditions for which a Fredholm integral equation of the first kind possesses a square integrable solution. However, since the underlying problem in polynomial basis conversion is not governed by this type of equation, it is desirable to motivate the discrete Picard condition directly from equation (6). This topic is now considered. In TSVD the smallest singular values are set equal to zero and equation (8) is rewritten as
The parameter lcil/si with i:
K (-4)
0 determines
(16)
the dependence
of the ratio
c.8cER”’
0 < 0 < 1 = Jc,\/si increases
which shows that the minimum and maximum values of the sensitivity function, taken over all vectors c and SC, are K-~(A) and K(A), respectively. Equations (14a) and (14b) show that as the dominant components of b = UC move from the directions that are defined by the eigenvectors of AA* which are associated with its large eigenvalues, to those that are defined by its small eigenvalues, the expressions in these equations decrease. In particular they attain their maximum (minimum) values when b is aligned along the eigenvector of AA* that is associated with the largest (smallest) eigenvalue of this matrix.
4. The discrete Picard condition Singular singular method algebraic be used
560
value decomposition forms the basis of truncated value decomposition (TSVD), which is a popular for the regularisation of ill-conditioned linear equations. It will be shown that TSVD can only for the regularisation of a particular class of
Appl.
Math.
Modelling,
1997,
Vol. 21, September
monotonically, i = l(l)m
0 = 1 * lcil/si = 1, i = l(l)m 1 < 8 =j lcil/si decreases
monotonically,
i = l(l)m
Since the ratio Icil/si may span several orders of magnitude, it is more convenient to consider the ratio log,Olcil/si, and this will be called the Picard ratio. Although the model in equation (16) is simpler than is likely to occur in a real problem, it is adequate for the motivation of the discrete Picard condition. Figure 2 shows three different forms of the Picard ratio, and Figure 3 shows the variation of the filter factor f, with i. It follows from equation (15) that the component of x(k) along the direction vi is equal to the product of the ordinate values of the graphs in Figures 2 and 3 at the abscissa value i =j. It is seen from Figure 262) that for 0 G 0 < 1, a large portion of the solution x(m) lies in the region of R” that is spanned by the right singular vectors that are associated with the large values of i, that is, the effects of the small singular values dominate the expres-
Polynomial
Picard
---_ ---_ -----)__ -----_____!fl_________.
i
j Figure
2.
ratio with
Three
different
the index
J. R. Winkler
sarily monotonically. However it is clear that the more rapid the decay of the ratio Icil/si, the better the approximation of x(k) to x(m) for k
ratio
--._
basis conversion:
forms
for the variation
of the
i, (a) 0 d 0 < 1, (b) tI= 1, and (c) fl>
Picard 1.
sion for x(m). This implies that the truncation of the series in equation (8) at an index k
Definition: The discrete Picard condition-The unperturbed right-hand side b of equation (6) satisfies the discrete Picard condition if for all numerically non-zero singular values si, the absolute values of the components ci decay towards zero, on average, more rapidly than do the singular values. It should be noted that the discrete Picard condition need only be satisfied approximately, that is, the ratio IciI/sj must decrease “on average” to zero and not neces-
filter factor 1.
A(x+iSx)=b+6b
(17)
is solved, and this problem does not satisfy the discrete Picard condition. However if b satisfies the discrete Picard condition, then k can be chosen so that b + 6b also satisfies the discrete Picard condition, from which it follows that a useful solution of equation (17) can be obtained. In the presence of noise, the components Ici + 6ciI/si typically decrease gradually until they reach a minimum at a level determined by the errors, and then increase. (This will be shown in the examples in Section 6.) It follows that x(m) is likely to be dominated by the noise corrupted terms that are associated with the small singular values. These noise corrupted terms must be filtered out if a useful solution to equation (17) is to be obtained. The correct choice of k allows the ratio filci + 6ciI/si to satisfy the discrete Picard condition, that is, the filter factor can be chosen so that the components of x(k) along the directions vi, i = 1.. . m, that are associated with the small singular values are zero. It is shown in Ref. 18 that when the discrete Picard condition is satisfied, equation (6) is ill-conditioned, but the converse is not true because of the situation that prevails when 0 = 1. This result, which is obtained by substituting equation (16) into equation (14b), emphasises the importance of establishing that the discrete Picard condition is satisfied and clarifies the statement that TSVD can only solve a particular class of ill-conditioned linear algebraic equation. 4.1 Statistical analysis The discussion of the discrete Picard condition above was adequate to suggest the importance of this condition, but more analysis must be performed in order to assert that it must be satisfied for a useful regularised solution to exist. This section and the subsequent section consider the discrete Picard condition in more detail. Let 6 b be a random vector with a diagonal covariance matrix proportional to the identity matrix
E{6b6bT}= (~‘1
(18)
where I E R”‘“, (T’ is the variance of each component of 6 b and E is the expectation operator. It follows from equation (lib) that i
j Figure
3.
The variation
of the filter
factor
with
the index
i.
E{ScGcT} = (T21
Appt. Math.
Modelling,
1997,
(19)
Vol. 21,
September
561
Polynomial
basis
conversion:
J. R. Winkler
Theorem 2-If the right-hand side vector b of equation (6) satisfies the discrete Picard condition and the noise vector 6 b satisfies equation (18), then TSVD reduces, on average, the effect of noise on the solution of equation (6).
The first term is identically becomes
k:l<-c;~$ j=l
Proof-If it is required (6), it is necessary that S(x(k),
that TSVD
regularise
b) - Stx(k - l), b)
IL$W =--IISCII
llcll
IlSX_ 1Wl
Ikll
Ils:cll
IISCII
IlSd- 14
be greater than zero. The term on the right-hand this expression simplifies to llcll
lIs:6cll
IIS~CII IIs:_ 141
and it is therefore
necessary
= IIs:~cl1211sx-
This expression
side of
I
that
1c11*- lls:cl1211s~_,6cl12
'j
1
(20)
(21)
It is clear that if the discrete Picard condition is satisfied, the expression in equation (21) is necessarily positive and TSVD reduces, on average, the sensitivity function of the minimum norm least squares solution of equation (6). However if the discrete Picard condition is not satisfied, the expression in equation (21) may be positive or negative. 4.2 The regularisation error There are two errors that must be considered when analysing TSVD. They are (i) the error in x(k) that arises from the error in b, .and (ii) the regularisation error that arises because x(m) # x(k), k # m, in the absence of noise. These two errors must not be confused; the first considers the effect of noise in b on the solution x(k), and the second is a measure of the bias that is introduced in x(k) in the absence of noise by the smoothing properties of TSVD. It has already been shown that if the discrete Picard condition is satisfied, then TSVD is effective in reducing the effect of noise on the minimum norm least squares solution x(k). The second type of error is now considered. The total error in this solution of equation (6) is
IIs:_ ICI1 - IIS%_ ,6cll IIS~CII
II WI [
5(k)
equation
zero, and thus equation
> 0
becomes
x(m) - (x(k) + 6x(k)) If the statistical assumption about 6c in equation (19) is substituted into this expression, the expected value of 5(k) is
Using equation written as
=x(m)
-At(b
+ 6b)
(6), it follows that this expression
(22) can be
x(m) - A+(Ax(m) + 6 b) = (I - A+A)x(m) - A%b
(23) E{c$(k)I
= a* (20)
The expression as
inside the square brackets
can be written
The bias that is introduced in x(k) by TSVD follows immediately because the expression in equation (23) does not approach the zero vector as S b + 0. The term (I A+A)x(m) is the error due to the regularisation procedure, and the second term is the error due to noise. The regularisation error e(k) is defined as e(k) = 11x(m) - x(k)11 IIx(m>II when 6 b = 0. It follows from equations
(22) and (23) that
e(k) = I~(I- A+A)x(m)ll II x(m) and this simplifies
II
The analysis of this expression is most easily performed assuming the simple model for (cj in equation (16).
to
by
Theorem 3-If the simple model for lcil in equation (16) is assumed, the regularisation error of the minimum norm
562
Appl.
Math.
Modelling,
1997,
Vol. 21, September
Polynomial
least squares solution
of equation
(6) satisfies
log
Proof-See (Hansen,” theorem 3.1). This theorem shows that if 8 > 1 and ski 1 is sufficiently small compared to sr, then x(k) is guaranteed to approximate x(m). The truncation parameter k is usually determined such that the larger the error in b,the smaller the value of k, that is, more singular values must be deleted. It follows that if there is an error in b, 0 must increase to offset the smaller value of k, so that a small regularisation error is retained.
discrete
magnitude
6.
-I 4.
A schematic
log diagram
I The
equation Picard
log
L-curves (18)
for (a) zero and
(b) only
mean signal
residual
random that
noise
that
satisfies
the
condition.
condition, and the other is the L-curve when b represents zero mean random noise that satisfies equation (18). These two separate L-curves are shown in Figure 5. The graph in Figure 4 is obtained by adding the ordinate values of the graphs in Figure 5 at each abscissa value. When very few singular values are deleted (segment AB in Figure 4) the L-curve is almost horizontal, but the gradient decreases rapidly as more singular values are deleted until the corner C is reached. In the vertical portion BC, x(k) is dominated by the effects of errors in b, magnified by the division by the small singular values. The deletion of these singular values causes a large change in log,,llx(k)ll, but the change in log,,l~(k) - bll is small. The gradient of the curve then increases rapidly, and the L-curve is approximately horizontal (segment CD in Figure 4). The change in logrollx(k)ll in this segment is small, but the change in log,,l!Ax(k) - blI is large. The gradient of the L-curve then decreases as more singular values are deleted, corresponding to the segment DE in Figure 4. The optimum value of k is its value at the corner C, because this provides the correct balance in keeping both Ilx(k)ll and I!Axtk) - bll small; if fewer singular values were to be deleted, x(k) would be dominated by the effects of noise, but if more singular values were to be deleted, x(k) would be smoothed too much. The log-log scale is used for two reasons: (i) The singular values span several orders of magnitude, and (ii) logarithmic axes emphasise portions of the L-curve where the variation in one variable is small compared to the variation in the other variable, and this has the effect of emphasising the corner.
An important issue that must be considered when using TSVD is the choice of the truncation parameter k; if it is too large, then x(k) may still be corrupted by noise, but if it is too small, a significant portion of the desired solution may be removed from x(k). This section describes the L-curve method,‘9,20 which is a simple graphical technique for the determination of the optimum amount of smoothing. Other methods,‘9*20 including the discrepancy principle and generalised cross-validation, exist, but they suffer from disadvantages. The motivation for the regularisation of equation (6) is to limit Ilx(k)ll, and thus a criterion for the optimum choice of k must include information on I!dk)ll and ILAx - bll. The L-curve is a parametric plot of log,,lMk1ll against log,,l~(k) - bll. It is shown in Refs. 19 and 20 that if (i) the unperturbed right-hand side vector b of equation (6) satisfies the discrete Picard condition, (ii) the perturbation S b is a zero mean random vector that satisfies equation (181, and (iii) llsbll +K llbll, then the L-curve assumes the shape shown in Figure 4. It consists of a sequence of discrete points that may be joined up to form a continuous curve. This curve is a combination of two separate L-curves, one of which is the L-curve when b is not corrupted by noise and satisfies the discrete Picard
Figure
5.
satisfies
5. The L-curve
log
J. R. Winkler
magnitude
I Figure
basis conversion:
Examples
This section contains two examples that demonstrate the theory of the previous sections. Example 1 considers the conversion from the power basis to the Bernstein basis, and Example 2 considers the conversion from the Bernstein basis to the B-spline basis. The notation of equation (6) that was used in the previous sections is replaced by the notation of equation (1) in Example 1 and by the notation of equations (7a) and (7b) in Example 2.
residual
of an L-curve.
Appl.
Math.
Modelling,
1997,
Vol. 21, September
563
Polynomial
basis conversion:
J. R. Winkle/
Example 1. Let the vector r in equation (1) store the coefficients of a univariate polynomial when it is expressed in the power basis, and it is required to compute the vector y, which stores the coefficients of the polynomial when it is expressed in the Bernstein basis. The degree of the polynomial is 12, which implies that m = 13. The vector r was chosen to satisfy the discrete Picard condition by specifying 0 = 2. Zero mean Gaussian noise was added to r, and its variance was such that Ilr II/11 6rll= 105.Figures 6(a) and 6(b) show the variation of the ratios log,,lcil/si and log,,lci + 6cil/si with i. It is seen that the ratio Icil/si decreases monotonically with i, but that the noise causes the ratio lci + &,l/s, to increase, albeit not monotonically, for i > 6. It is clear that the direct solution of equation (1) with no regularisation would yield an unsatisfactory answer because of the significant effect of the noise. Figzue 7(a) shows the polynomial that is being transformed, in the range 0 < t < 1, and Figure 7(b) shows the transformed polynomial when equation (1) is solved directly and r is corrupted by noise. The curves in these two figures should be similar, but it is clear that the noise in r severely degrades the solution of equation (1). Figure 8 shows the discrete points that define the L-curve, and it is seen that the corner is rounded, i.e., not well defined, thus emphasising the need for a rigorous definition of the corner of the L-curve. The figure suggests that k = 6,7,8, or 9 because these indices define the points that are closest to the corner of the L-curve. Figzue 6(b) suggests that the summation in equation (8) should be truncated at k = 6, corresponding to the deletion of seven singular values. Figures 9(a), 9(b), and 9(c) show the (original) power basis polynomial and the (transformed) Bernstein basis polynomial for k = 5, 6, and 7, respectively. It is clear that the transformed polynomials for k = 5, 6 are similar and that they differ from the transformed polynomial that is obtained for k = 7. It is seen that the optimum value of k is 6, but k = 5 also yields an acceptable solution. The transformed polynomial is improvement in the clear-compare Figures 7(b) and 9(b).
Picard
Exact
solution
(a) Noise
corrupted
solution
(4 Figure 7.
(a) The
power
basis
to the
formed
polynomial
polynomial
whose
Bernstein
basis
when
r is corrupted
transformation is desired.
from
(b) The
the
trans-
by noise.
Figure 10 shows a sequence of L-curves for four different values of the ratio of the signal power to the noise power for a given vector r that satisfies the discrete Picard condition. As this ratio decreases, the corner of the L-curve becomes more pronounced and the optimum value of the truncation parameter decreases. This is because the portion of the ratio jci + 6cil/si that decreases terminates at a lower value of i, since the noise 6ci corrupts a larger number of the components ci, i = 13,12,. . . .
ratio
log magnitude
':::I'"'"L;log
1,
i
4.99
4.96 Figure 6.
The
variation
log,,lc,+Gcj~sj
with
564
Math.
Appl.
of the
ratios
(a)
loglo~ci~si
and
(b)
4.91
Figure 8.
i.
Modelling,
1997,
Vol. 21, September
I
-0
.
k=4
I
The L-curve.
residual .
k:2
.
Polynomial basis conversion: J. R. Winkler
k=5
log magnitude 7.06.5, (c)
6.0 5.5, .:I
-.. I
0 Figure
10.
(a)
(b) ) z
2
4
The L-curves
to (a) 30, (b) 10”.
t 6
8
for (signal
log
residual
power)/(noise
power)
equal
(cl 1Oa. and fd) 10’.
(4 k=6 Bezier
curve
0.29 0.2
u-6
1
-100 0.28
0.27
0.26
0.25
t=o 0.24 Figure
0.29 11.
0.3
The Bezier
0.31
curve
0.32
and the convex
0.33
0.34
hull of its control
polygon.
Figure
9.
(a) The
exact
for k= 5. (b) The exact k=6.
fc) The
exact
solution,
solution,
solution,
-, -,
-,
and TSVD
and TSVD and
TSVD
solution,
solution, solution,
----,
B-sgline curve
- - - -, for ----,
0.5
for
f
k=7. 0.4.
Example2. Consider the transformation of a Bezier curve of degree 12 (m = 13) to a B-spline curve. The matrix F in equation (4) stores the coordinates of the vertices of the control polygon of the Bezier curve, and it is required to compute the coordinates of the vertices of the control polygon of the equivalent B-spline curve. These coordinates are stored in the matrix H. Figures 11 and 12 show, respectively, the Bezier and B-spline curves and the convex hulls of their control polygons. The - in these (and in Figures 14 and 17(b))
Figure
The
12.
B-spline
curve
and the convex
hull of its con-
trol polygon.
Awl.
Math.
Modelling,
1997,
Vol. 21, September
565
Polynomial
basis conversion:
J. R. Winkler
denote the vertices of the control polygons. It is seen that although the curves are necessarily identically equal, the convex hulls are different. In particular, the convex hull of the control polygon of the B-spline curve is substantially larger than is the convex hull of the control polygon of the Bezier curve.21 Zero mean Gaussian noise was added to the right-hand side of equations (7a) and Ub), and its variance was adjusted so that
II II = 2.12 x 103 11 i?(GF)lll
II(GF)~II 116(GF)z II
B-spline
basis
= 1.07 x 103
Figures 13(a) and 13(b) show the variation of the ratios log,,lcjl/si and log,,lc, + tk,l/s, with i for equations (7a) and (7b), respectively. These figures show that the direct (simple) solution of equations (7a) and (7b) will yield unsatisfactory answers, because the largest values of the terms lci + 6cil/si are Figure 14. The convex hull of the control B-spline curve in the presence of noise.
polygon of the
4
Picard ratio
W
a
-6 * -8A Figure 13.
(A) The ratios (a) for equation (7a).
lOgloICiVS/
and (b) log,olci+acivsi
approximately equal to 16cil/si. This is confirmed in Figure 14, which shows the computed convex hull of the control polygon of the B-spline curve when no regularisation is applied to the equations. The cluster of points near the origin of this diagram are some of the vertices of the control polygon of the curve. Comparison of Figures 12 and 14 shows that the effect of the noise on the convex hull is substantial. Figure 15 shows the B-spline curve in this circumstance, computed from p(t) + Sp(t) = tJ(H+ 8H),and a comparison of Figures 11 and 15 shows that the change in the curve is negligible. Although the fidelity of the curve has been retained in the transformation, the gross magnification and distortion of the convex hull implies that the transformation is unsatisfactory. Figures 16(a) and 16(b) show the discrete points that define the L-curves for equations (7a) and Ob), respec-
4-
Picard ratio a,
B-spline
0.301
basis
0.26
0.24~
-8.
t=0
B Figure 13. (B) The ratios (a) for equation (7b).
566
Appl.
lOg,,ICiVS;
Math. Modelling,
0.28
and (b) log,cIci+aciVs; Figure 15.
1997, Vol. 21, September
0.30
0.32
0.34
The B-spline curve in the presence of noise.
Polynomial basis conversion: J. R. Winkler log
magnitude
k=6
J.Or l
2.0
1.0 0.0
0.29
0
0.28
a l
-6.0
-4W
46.0
l
0.
&og
0.21’
residual
0.26,
-1.01 (a)
0.25.
log
magnitude
t=o 0.24’
0.29
0.3
0.31
0.32
0.33
0.34
W
03 Figure 16.
(a) The L-curve for equation (7a). (b) The L-curve for equation (7b).
0.3
tively. Detailed examination of the curves showed that the values of k nearest the corner of the L-curve in Figure 16(u) are 5, 6, 7, and 8, and that the values of k nearest the corner of the L-curve in Figure 16(b) are 4, 5, 6, 7, and 8. Figures 17(a) and 17(b) show the B-spline curve, and the B-spline curve and convex hull of its control polygon, for k = 6. It is seen that very good answers are obtained for this value of k. Very similar answers were obtained for k = 5. The solutions for k = 7, 8 were underregularised. In general equations (7a) and 0%) may require different amounts of regularisation, that is, different values of k, or one of the equations may require regularisation and the other one may not. It is for this reason that equation (4) cannot be treated as a single equation with a multiple right-hand side. Rather, it must be decomposed into the form given in equations (7a) and G’b), and each of these equations must be solved individually. The error measure in the L-curve in Example 1 is a function of the transformation matrix P and coefficients r of the polynomial when it is expressed in the power basis. It may be desirable to measure the error by a property that is intrinsic to a curve because (i) this has more geometric meaning, and (ii) it is well known that small changes in the coefficients of a polynomial that is expressed in the power basis may alter it substantially-the Wilkinson polynomial is a common example. It follows that a small residual may not guarantee that the polynomial in the receiving CAD system is a good approximation to the polynomial in the sending CAD system. A similar situation arises for the problem in Example 2; the residual is a function of the transformation matrix J and the coordinates of vertices of the control polygon of the Bezier curve, but it may be advantageous to define the error by the proximity of the curve to its convex hull.
’
0.1
0.3
0.2
0.4
0.5
(b) Figure 17. (a) The B-spline curve for k=6. (b) The B-spline curve and convex hull of its control polygon for k=6.
Acknowledgment The author acknowledges the support of the UK EPSRC by an Advanced Research Fellowship, reference number B/93/AF/1703.
References 1. Farouki, R. T. and Rajan, V. T. On the numerical condition of polynomials in Bernstein form. Comput. Aided Geo. Design 1987, 4, 191-216 2. Cohen, E. and Riesenfeld, R. F. General matrix representations for Bezier and B-spline curves. Compuf. Indust. 1982, 3, 9-1.5 3. Rogers, D. F. and Adams, J. A. Mathematical Elements for Computer Graphics, McGraw-Hill, New York, 1990 4. Yamaguchi, F. Curves and Surfaces in Computer Aided Geomettic Design, Springer-Verlag, Berlin, 1988 5. Grabowski, H. and Li, X. Coefficient formula and matrix of nonuniform B-spline functions. Comput. Aided Design 1992, 24, 637-642 6. Lyche, T., Morken, K. and Strom, K. Conversion between B-spline bases using the generalized Oslo algorithm. Knot Insertion and
Appl.
Math.
Modelling,
1997,
Vol. 21, September
567
Polynomial
7.
8.
9. 10. 11.
12.
13.
basis conversion:
J. R. Winkler
Deletion Algorithms for B-Spline Curues and Surfaces, eds. R. N. Goldman and T. Lyche, SIAM, 1993, pp. 135-153 Lachance, M. A. Piecewise polynomiai approximation of poiynomiai curves. 2nd inte~ational Conference on Algorithms for Approximation, Royaf Military College of Science, Shrivenham, UK, July 1988, eds. J. C. Mason and M. G. Cox, pp. 125133 Daniel M. and Daubisse, J. C. The numerical problem of using Bezier curves and surfaces in the power basis. Comput. Aided Geo. Design 1989, 6, 121-128 Hanke, M. and Hansen, P. C. Regularization methods for largescale problems. Surti. Mark. Ind. 1993, 3, 253-315 Hansen, P. C. The truncated SVD as a method for regularization. BIT 1987, 27, 534-553 Hansen, P. C. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM J. Sci. Stat. Comput. 1990, 11, 503-518 Hansen, P. C. Numerical tools for analysis and solution of Fredhoim integral equations of the first kind. InrJerse Problems 1992. 8, 849-872 Varah, J. M. A practical examination of some numerical methods for linear discrete ill-posed problems. SIAM ReL). 1979, 21, 100-111
568
Appt.
Math. Modelling,
1997, Vol. 21, September
14. Stewart, G. W. Introduction to Matrix Computations. Academic Press, New York, 1973 15. Golub, G. H. and Van Loan, C. F. auto Computu?ions. Johns Hopkins University Press, Baltimore, MD, 1989. 16. Winkler, J. R. The sensitivity of linear algebraic SIAM Conference on Applied Linear Algebra, USA, June 15-18 1994, pp. 279-283 17. Hansen, P. C. The discrete Picard problems. RIT 1990,30,658-672
condition
equations, Snowbird,
for discrete
Fifth Utah,
ill-posed
18. Winkler, J. R. The condition number of a matrix and the stability of linear algebraic equations. Department of Computer Science, The University of Sheffield, memoranda CS-97-05 19. Hansen, P. C. Analysis of discrete ill-posed the L-curve. SIAM Reu. 1992, 34, 561-580
problems
by means of
20. Hansen, P. C. and O’Leary, D. P. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM .I. Sci. Comput. 1993.14, 1487-1503 21. Goult, R. J. Standards for curves and surface data exchange, Fourth IMA Conference on The Mathematics of Surfaces, Bath, England, September 1990, ed. A. Bowyer, 1990, pp. 1-15