- Email: [email protected]
Determination of shape preserving spline interpolants with minimal curvature via dual programs
JOURNAL
OF APPROXIMATION
52, 43-57 (1988)
THEORY
Determination of Shape Preserving Spline Interpolants with Minimal Curvature via Dual Programs SIEGFRIED
DIETZE
W.
AND JOCHEN
SCHMIDT
Deparrment of Mathematics, Technical University of‘ Dresden, Mommsenstraje 13, Dresden, 8027, German Democratic Republic Communicated by G. Meinardus Received May 1, 1985; revised February
1.
13, 1987
INTRODUCTION
In this paper programming problems that arise in the convexity or monotonicity preserving spline interpolation will be studied from a numerical point of view. The problems appearing there are of the following special structure. PROGRAM
PA. Minimize subject to
(i== 1, ...) n),
(nz,- , , mi) E W,
where Fi: R2 + R is a strictly convex function and Wi c R2 a closed convex set. For solving PA computationally here it is proposed to pass to the dual PROGRAM
DA. Maximize with
(1.2)
PO= P,!= 0.
The function H,*: R2 -+ i? is therein to be defined by H,*([,
~/)=sup{tx+qy-F,(x>
Y):
(x>
y)E
f+‘i).
(1.3)
In DA no constraints occur and, moreover, the gradient of the objective function is tridiagonal. These properties are indeed very convenient for 43 002 I -9045/88 $3.00 Copyright r 1988 by Academic Press. Inc All rights of reproduction in any farm reserved.
44
DIETZE AND SCHMIDT
numerical purposes and should be utilized by practical algorithms in order to be efficient. Attacking PA directly, e.g., by an active set strategy (see, e.g., Fletcher [Z]) is in general more expansive than starting from DA. Further, if a solution of DA is known then, under suitable assumptions, a solution of PA is explicitly given by the return-formula mr--I=alH*(Pi-l, mi=
-Pi)
a2HT(Pi-
12
(i = 1, ...) n).
-PiI
(1.4)
In addition, for problems arising in the shape preserving spline interpolation, the program DA can be explicitly stated becausethe subprograms (1.3) are solvable by formulas, and strong duality theorems are valid. Following this line the results of Burmeister et al. [l] concerning the convex spline interpolation are confirmed deductively. Moreover, also for the monotone spline interpolation problem, analogous results can be offered now that likewise are obtained in a deductive manner by applying the above-mentioned general concept. For further applications see [ 111, C131, C141, Cl51, and C161. 2. FORMULATION OF THE SPLINE PROBLEMS For describing the needed facts from the spline interpolation denote by d a partition of the unit interval [0, 11, A:x,=O
...
and let Sp(3, A) be the set of cubic Cl-splines on A. The points (x,, JJ,,),.... (x,, y,,) associated with A are said to be in convex position if T,
...
(2.1)
where si = ( yi - y, ,)/hi, hi = xi - x,+ , They are in monotone position if T, 2 0, ...) ?,,> 0.
(2.2)
Spline interpolants s E Sp(3, A) should satisfy the concrete interpolation condition 4x,) = Y,
(i = 0, ...) n).
(2.3)
Further, denote by m, = s’(x;) the first derivative of s in the node x,. Now, a result of Neuman [S] states that under (2.1) a spline interpolant s E Sp(3, A) is convex on [0, l] if and only if (m, I) mi) E W
(i= 1, .... n),
(2.4)
SHAPE
PRESERVING
SPLINE
INTERPOLANTS
45
where Wi = K, is the cone Ki={(X,y)ER2:2~+y~3~,,x+2y~3tj}.
(2.5)
An analogous result for the monotone spline interpolation, under the assumption (2.2), is due to Fritsch and Carlson [3]. Here Wj= Eiu D, is taken as the union of the ellipse E;= ((x, ,v)ER? x2+xy+y2-6q(x+y)+9tf
(2.6)
with the square D,= ((x, y)~R2:06x~32,,0~~~~35,}.
(2.7)
Naturally for a given set PPZ,,, .... m,, of first derivatives the spline interpolant s E Sp(3, A) is uniquely determined and can be directly computed. In the convex case the problem (2.4), with Wi= Ki, has in general no solution. An effective algorithm by which the solvability can be decided is given by Schmidt and Hess [7]; see also [9], [lo], and [12]. If at all, (2.4) has either one solution or an infinite number of solutions. In the monotone case the problem (2.4), with W, = E, u Di, is always solvable but in general not uniquely. Therefore, both in the convex and in the monotone spline interpolation case the question arises which solution of (2.4) i.e., which of the corresponding splines shall be chosen. A good strategy is to select that spline which has minimal mean curvature
I 5 s”(xJ2 dx = 0
i ;=,
i { rnf- 1+ mi_ ,m, + rnf - 3T,(m,-, + mi) + 3rf). I
X ?
FIG.
1.
The cone K.
(2.8)
46
DIETZE AND SCHMIDT
FIG. 2.
The set E,u D,
Thus one is led to a programming function Fj(X,
y)=R
I
X
37i
T
problem
of the type PA, with the
{X2+,uJ’+4.2-3ti(-\-+y)+32~)
(2.9)
and the sets Wi=K,
or
Wi=E,uDi
introduced by means of (2.5)-(2.7). The Hessian of the function (2.9) is positive implies
(2.10)
definite. This property
PROPOSITION 1. The program PA and the auxiliary programs (1.3) are uniquely solvable if speczfied by (2.9) and (2.10) and {f the feasible domains are nonempty. The optimal values are finite.
3. APPLICATIONS
OF FENCHEL’S
THEORY
In order to verify the duality of the above-mentioned pair PA and DA of programs and to derive duality statements one can refer to Fenchel’s theory and to extensions of it due to Rockafellar and Stoer, see, e.g., [4, 6, 81, but also Lagrange’s duality theory may be applied. Letf,g: R”-+a=Ru{+a, - cc } be proper convex functions, that is, f and g are convex functions greater than - cx, whose effective domains dom f and dom g are nonempty. Moreover, let ,f and g be closed. As usual,
47
SHAPE PRESERVING SPLINE INTERPOLANTS
the subdifferential off is denoted by af while the conjugate function will be called f*, f*(u)=sup{Uru-f(U)XER”‘}.
(3.1)
Further remember that for the indicator function I of a set W is I(u; W) = 0 for u E W and Z(u; W) = +ar, for u $ W. The subsequently used concepts of stability are defined, e.g., in [S]. Now, in Fenchel’s theory, to the primal convex program P: inf{ f( u) + g(u): u E R”’ }
(3.2)
belongs the dual concave program D: sup{ -f*(r)-g*(-u):
UER”j.
(3.3)
A corresponding strong duality theorem of a generality necessary here reads as follows (see [8, Section 5.31 for the first part and see,e.g., [6, Section 311 for the second part). THEOREM
2.
Let f and g be closed proper convex functions. Assume that
inf{f(u)+g(u):
UE R”} E R
(3.4)
and that there are vectors w, z E dom f n dom g such that f is w-stable and g is z-stable. Then, program D has an optimal solution, and the optimal values of P and D are equal:
inf{f(u)+g(u):
MER”} = maxi-f*(u)-g*(-u):oER”).
In addition, u and v are optimal for
(3.5)
P and D, respectively, jf and only if
u~df*(v)ndg*(-v)
(3.6)
holds.
Next, the programs PA and DA are shown to be special casesof P and D while the verification of the assumptions of the preceding theorem including the application of formula (3.6) will be undertaken later and only for the concrete spline problems. By introducing additional variables, Program PA can be reformulated as PA’: Minimize
i I=
Fi(mip , , Mi) 1
subject to (m,-,, M;)c
W,(i=
1, .... n), Mi=m,
while M, = m, is only a change of notation. 640’52’1.4
(i= 1, .... n- l),
48
DIETZE AND SCHMIDT
Naturally this is a separable program treated by several authors, see,e.g., [6, Section 311. Using the function Hi: R2 -+ R, Hib, y) = Fib, 4’) + 4x, y; W,), set, with u=(m,,M,,m
,,..., M,-,,m,-,, f(u)
=
f
i= I
H,(m,-
(3.7)
M,,)ER~“, I,
(3.8)
Mi).
Further let
g(u)=4%w,
W={u~R~“:M~=m~(i=l,...,n-1)).
(3.9)
By these definitions problem PA’ and therefore likewise PA become a program P. The occurring functions f and g are proper convex and closed. In order to state the corresponding dual program D the conjugate functions f* and g* to (3.8) and (3.9) are to be determined. Let u = (PO,P,,p,, .... Pnp,,pnm ,, P,*)E R’“. Then, in view of (3.8) it follows immediately that f*(U)=
i
,=I
H*(pi-
(3.10)
1, pi)
holds if H,* is given by (1.3). Observe further that g*(u) = sup{ urzX2.4 E W}
=o
for po=O,P,+pi=O(i=l
=+cc
otherwise
,..., n),p,,=O
(3.11)
and therefore g*( -u) = g*(u). Hence, problem D reduces to program DA. Now, the formula (3.6) shall be specified. To this end, it is additionally assumed that the auxiliary programs (1.3) are solvable for all ([, q) E R2 with finite optimal values HT([, q). Becauseof the strict convexity of F, the maximizer
in (1.3) is unique. Further, since (xi, yi) E i3H,+(<, q) if and only if (xi, yi) is a maximizer of (1.3) one gets the differentiability of HF and (3.12)
49
SHAPE PRESERVING SPLINE INTERPOLANTS
In addition, the derivatives a, H,? and & H,f+ are continuous. For this, see, e.g., [6, Sections 23 and 251. Thus, f* defined by (3.10) is continuously differentiable. Further, let v be optimal for D and let the optimal value be finite. This implies g*( -v) = 0, and the Kuhn-Tucker conditions for program D, i.e., for inf{f*(v):p,=O, P,+p,=O (i=l,..., n-l), P,,=O}, yield ?,H,*,,(p,, vector u by
P,+,)=8,H~(p,~,,
m,~ I=~lHiY~,-,,Pj),
P,), i= 1, .... n- 1. Now, define the
Mi=&H,%-
i = 1, .... n.
,, P,),
Then one gets g(u) = 0, v*u = 0, and hence g*( -v) + g(u) = -uTu. Thus it follows that UE 8g*( -u), see, e.g., [6, Section 231. Further, because of u E $f*(c), relation (3.6) is valid. Thus, the mentioned assumptions imply PROPOSITION 3. Let (m,, m, , .... m,,) with
(p,,, p,, .... p,,)
he
m,- l = 8, K+YP,-, , -pi) = fitpi-
a
solution
of program
DA
I) -PA
mi = azHI*(p;- 1) -P;) = YAPi- 13-P;) is the only solution continuous.
of’
(i= 1, .... n)
then
(3.13)
PA. The derivatives 2, H,*, 8, H,* are
4. CONVEX SPLINE INTERPOLATION As said above, for determining the convex spline intepolant with minimal mean curvature one has to solve the PROGRAM
PC. Minimize subject to
(m, _, , mi) c K,
(i= 1, .... n)
(4.1)
with Fi and Ki according to (2.9) and (2.5), respectively. By Proposition 1 this program is uniquely solvable if the feasible domain K= {(m,,m,,
.... m,)ER”+‘:
(mip,,mi)EK,(i=
is not empty. Note that K # @ implies (2.1).
l,..., n)}
(4.2)
50
DIETZE AND SCHMIDT
For numerical purposes the corresponding dual program DC shall be stated. To this end, the subprogram (1.3) now maximize d(x, u)=sx+qy-i
{(x-T)~+(x-~)(J--T)+(J~--)~) (4.3
subject to 2x + JJ< 32, x + 2~133r,
where h = h, and 7 = z,, is to be treated for all (t, q) E R2. In view of (1.21 and (1.4) especially, the optimal value d,,, = H,*([, q) is of interest. Proposition 3 assures that H,* is continuous differentiable. PROPOSITION 4.
The optimal value of program (4.3) is equal to
2 =r,(t+q)+k
6-i (
=7;(5+q)
)
.for q d 0, 21”6 ~1
.for 5 2 2q,21”3
FIG. 3. Pieces of the domain of function (4.4)
q.
(4.4)
SHAPE PRESERVING SPLINE INTERPOLANTS
51
The proof of (4.4) is elementary and is therefore only sketched. In the first case, the center of the ellipse d(x, y) = c, I=r+;(zc-a),
4$-w
(4.5)
is the maximizer of the problem (4.3). This is possible if and only if 5 ,< 0, q > 0, and the optimal value HT(t, q) = d(.Z, J) is as given in (4.4). In the second case, the maximizer is located on the line 2x +J = 3~, x < 5. This is valid only if ir > 0. The maximizer is the solution of 0,
Y) = c,
2x+y=3~,
where the constant c is to be chosen such that exactly one solution exists. One now gets
whenever X d z; that is, r - 2~ d 0. And HT(& q) = d(Z, j) with H,*(t, q) from (4.4) can be confirmed. The third case is analogous to the second. The maximizer is now lying on x + 2y = 32, y > z. The last case occurs for 4 - 2~ 2 0, 2t -q 2 0. Then the maximizer is in the vertex of the cone K,, X=j=r, and H*(5,11)=d(x,y)=t;(~+rl). Thus the proof is complete. Using Proposition 4 and (1.2) one is led to the following program DC which is dual to PC. PROGRAM
DC.
Maximize - i i=
I
HT(p,+ , , -Pi)
with
p. = pn = 0,
(4.6)
where H,? is given by (4.4). Next, by means of Theorem 2 the solvability of DC is shown. For this, let the feasible domain K of PC not be empty. Obviously, for (m,, m,, .... m,) E K it follows that z=(m,,m,,m,,
.... m,,-,,m,~,,m,)Edomfndomg
if, with Wi= K,,fand g are given by (3.8) and (3.9), respectively. Now, by stability statements obtainable from [IS, Sections 4.8, 5.3, and 5.51 the functions are proven to be z-stable. Indeed, F, considered as a function defined in R2” is z-stable, since z is an interior point of dom Fi. Further, in view of It-x, Y; K,) =1(x, y; K,,) + 4x, I’; K,,)
52
DIETZE AND SCHMIDT
with Kn={(x,y):2x+y<3r,}, K,j={(x,y):x+2y>3~,} and, since I( ., K,,), I( ., K,,) in spite of the affine constraints are polyhedral, the function Z( ., K,) is z-stable. Thus, this property holds also forf: In the same manner g is shown to be z-stable. Therefore, by the Theorem 2, program D possessesan optimal solution and thus likewise Program DC. Summarizing the main result of this section concerning the convex spline interpolation reads of follows. THEOREM 5. Both Programs PC and DC are solvable if and only if the feasible domain K of PC is nonempty, and then the optimal values are equal. Furthermore, if a solution (po,pI, .... p,) of DC is known, the solution (m,, ml, .... m,) of PC is given by formula (3.13) where HT is to he defined by (4.4).
As mentioned in the Introduction this theorem can be found already in paper [ 11. There it is derived in a more direct way be starting from the KuhnTucker conditions. Numerical tests presented in [l] show the strategy of solving PC via DC to be very efficient even if DC is handled by the ordinary Newton method; see also [ 111, [ 141, and [ 161.
5. MONOTONE SPLINE INTERPOLATION
As stated before, the monotone spline interpolation leads to PROGRAM PM.
Minimize
i
F,(m,- 134)
,=l
subject to (m,- , , mi) E Ei u D;
(i= 1, .... n),
(5.1)
where F;, E;, and Di are given by (2.9), (2.6), and (2.7). Here the feasible domain is always nonempty. Therefore, in view of Proposition 1, PM is always uniquely solvable but, as in the convex case, the solution can be computed more effectively by means of the corresponding dual program DM. For formulating this program the subprogram (1.3), now maximize d(x, y)=[x+qy-i{(x--T)2+(~--i)(y-T)+(y-r)2} (5.2) subject to (x, y)~ E u D E = E,, and D = D,, is to be considered for all with h=hi, T-T,, (5, v)eR*.
SHAPE
PRESERVING
SPLINE
53
INTERPOLANTS
PROPOSITION 6. The optimal value d,,, = HF(<, q) of program (5.2) is, with ai= 12zj/hi, given by
for (5,rl)EAiuBj for (5, v) E S,,
for (i", vIESI,
for (5,V)ETi +2a,Jt2
- 5yl+ q2 - ai(t + q) + af
for
(5,
VIE
ci
(5.3) where A,= {(5,rl)ER2:-ai~2r-)1~2a,,
-2a,<[-2qda,},
Bj= {(tl,~)ER2:52-5~+~2-ai(5+vl)~0}, Sli=
((5,rl)ER2:t-2v>a,,
S2;= {((,ty)ER2:2[-q6
-a,65da,}, -ai,
-a,Gq
Tj = {(L q) E R2: 5 < -ai, q 6 -a,},
C~=R'\(A~UB~US,~US,,U
T,).
%i Ti
FIG. 4.
Pieces of the domain
of function
(5.3).
54
DIETZE
AND
SCHMIDT
The proof of (5.3) is elementary. Only a sketch will be given. Compare also the proof of (4.4). In the first case the maximizer is equal to the center (4.5) of the ellipse d(x, y) = c, valid if and only if (t, q) eAiu Bi. In the second case the maximizer (X, j) is lying on ,v=O, 0 bx 6 35. Then one gets
if (5, yl) E S,;. The third case is symmetric to the second. In the fourth case the maximizer equals X = j = 0 true for (5, q) E Ti. The last case occurs if the maximizer is located on the “outer” boundary of the ellipse e(x, y)=i
{x~+x~+~~-~~(x+,v)+~~~}=O.
Then, the maximizer (X, j) solves d(x, y) = c, e(x, y) = 0 if c is chosen to be the greatest value such that exactly one solution exists. Hence, by adding the two equations one gets
with c1= 5 - 127/h, p = q - 12r/h, y = c - 24T2/h. Eliminating J in e(x, y) = 0, a quadratic equation in x arises. This has only one solution if y2 - 4r( LX+ p) y + 12&/l = 0. The greater root is
Now, because of d,,, = d,,,,, + emax= y + 24T2/h the value of H,*(t, q) given in (5.3) immediately results if (5, 9) E C,. Thus the proof is complete. By means of Proposition 6 and (1.2) the following program DM is obtained which is dual to PM.
PROGRAM DM. Maximize - i H,+(p, _ , , -P,) i= 1 where H,* is defined by (5.3).
with
p. = p,! = 0,
(5.4)
SHAPE PRESERVING SPLINE INTERPOLANTS
55
In order to verify the solvability of DM suppose that, in sharpening (2.21, T1 > 0, ...) Tn > 0
(5.5)
holds. Then, e.g., the vector z = (m,, M,, m,, .... M,.. ,, m,, M,) E R2”
with m,=M;=E/2, s=min{r ,, .... r,} > 0 is an interior point of dom H, even if the function H,, defined by (3.7) with W, = E, u D,, is considered as a function given on R2”. Thus F, is z-stable. In view of (3.8) one gets z ~domf and the z-stability off: Further, for the function g defined by (3.9) obviously z E dom g. And g is z-stable because g can be represented by a sum of polyhedral functions. Therefore, Theorem 2 assures the existence of an optimal solution for Program DM. Summarizing, the main result of this paper concerning the monotone spline interpolation is the following. THEOREM 7.
Both programs PM and DM possess optimal solutions if (5.5) holds, and then the optimal values are equal. Moreover, if a solution (PO,PI, ...?p,) of DM is given, the solution (m,, ml, .... m,) of PM can he computed by means of formula (3.13), where HI is to be defined by (5.3).
6. PROGRAMS WITH GIVEN BOUNDARY VALUES In the monotone spline interpolation ti = 0
for some jE { 1, .... n}
(6.1)
may occur. Then, mjp, = mj= 0 for the solution of PM follows, and the optimization problem PM can be decomposed into subproblems of lower dimension. Analogously, if in the convex spline interpolation 5j=T,+
I
for someJo ( 1, .... n - 1 }
(6.2)
appears one gets m, ~ 1= mj = mi + , = TV, and Program PC divides into lower dimensional subprograms. Of course, from a numerical point of view such a decomposition should be undertaken if possible. In the resulting subproblems, one or both of the boundary values, say m, and m,, are prescribed. For example, let m, = z and m, = j be given. Then one has to consider the slightly modified
56
DIETZE AND SCHMIDT
PR~CRAM
PB.
Minimize
i
F,(m,- , , mi) I subject to (m,+, , mi) E Wj(i = 1, .... n),
(6.3)
i=
mo= M,m,,= fl.
Two ways are offered to derive a corresponding dual program. In the first case, define the functions fand g by (3.8) and (3.9) and let the set W now be W= {zLER~~: Mi=mj(i=l,...,n),m,=rx,m,=~}.
(6.4)
Then, it is immediately seen that for ~~+P,=O(i=l,...,n)
g*(u) = vo + PPll = +m
otherwise.
(6.5)
Using the functions HT given by (1.3) the dual program now obtained reads as follows. PROGRAM
DB,. Maximize - f
H,?(p,
,, -P,) + wo - BP,,.
(6.6)
r=l
In a second case, the variables m, and m, are eliminated, and thus PB is considered to be a problem in the space RZnp2. Following the line described in Section 3 one is led to the dual PROGRAM
DB,.
Maximize -G:( -p,)-
II- I 1 NT@,- ,, -pi)-G,*(p,_,).
(6.7)
,=2
Here H,* are defined by (1.3) and G:, G,* by G:(v) = sup{~ -P,(a, Y): (~9Y) E W, 1,
(6.8)
G,*(O=sup{tx-F,(x,
(6.9)
8): (x, B)E W,}.
In the shape preserving spline interpolation, if (6.1) or (6.2) occur, it is recommended to decompose the program for determining the interpolant with minimal curvature into lower dimensional subprograms. These are programs with boundary conditions of the type PB, to be solved via DB, or DB,. For the convex spline problem, Theorem 5 holds verbatim if PB is
SHAPE PRESERVING SPLINE INTERPOLANTS
51
coupled with DB, or with DB,, while for the monotone spline problem, an analog to Theorem 7 can be assured if PB is coupled with DB,. The details are omitted because their verification is now straightforward.
REFERENCES I. W. BURMEISTER,W. HESS, AND J. W. SCHMIDT, Convex spline interpolants with minimal curvature, Computing 35 (1985), 219-229. 2. R. FLETCHER, “Practical Methods of Optimization. Vol. 2. Constrained Optimization,” Wiley, New York/Chichester/Brisbane/Toronto/Singapore, 1981. 3. F. N. FKKSCH AND R. E. CARLSON, Monotone piecewise cubic interpolation, SIAM J. Numer. Anal. 17 (1980), 238-246. 4. P.-J. LAURENT, “Approximation et Optimisation,” Hermann, Paris, 1972. 5. E. NEUMAN, Uniform approximation by some Hermite interpolating splines, J. Compuf. Appl. Math. 4 (1978), 7-9. 6. R. T. ROCKAFELLAR, “Convex Analysis,” Princeton Univ. Press, Princeton, NJ, 1970. 7. J. W. SCHMIDT UND W. HESS, Schwach verkoppelte Ungleichungssysteme und konvexe Spline-Interpolation, E/em. Math. 39 (1984), 85-95. and Optimization in Finite Dimensions,” 8. J. STOER AND C. WITZGALL, “Convexity Springer-Verlag, Berlin/Heidelberg/New York, 1970. 9. P. COSTANTINI, On monotone and convex spline interpolation, Mafh. Compuf. 46 (1986), 203-214. 10. P. COSTAN~IN~ AND R. MORANDI, Monotone and convex cubic spline interpolation, Calcolo 21 (1984), 281-294. 1 I. W. HESSAND J. W. SCHMIDT, Convexity preserving interpolation with exponential splines, Computing 36 (1986), 335-342. 12. J. W. SCHMIDT, Convex interval interpolation with cubic splines, BIT 26 (1986), 377-387. 13. J. W. SCHMIDT, On shape preserving spline interpolation: Existence theorems and determination of optimal splines, Banach Center Publ., Semester XXVII, Polish Scient. Publishers, Warsaw, in press. 14. J. W. SCHMIDT, On tridiagonal linear complementarity problems, Numer. Mafh. 51 (1987), 11-21. 15. J. W. SCHMIDT, An unconstrained dual program for computing convex C’-spline approximants, Compufinz, in press. 16. J. W. SCHMIDT AND W. HESS, Positivity of polynomials on intervals and positive spline interpolation, submitted.
OF APPROXIMATION
52, 43-57 (1988)
THEORY
Determination of Shape Preserving Spline Interpolants with Minimal Curvature via Dual Programs SIEGFRIED
DIETZE
W.
AND JOCHEN
SCHMIDT
Deparrment of Mathematics, Technical University of‘ Dresden, Mommsenstraje 13, Dresden, 8027, German Democratic Republic Communicated by G. Meinardus Received May 1, 1985; revised February
1.
13, 1987
INTRODUCTION
In this paper programming problems that arise in the convexity or monotonicity preserving spline interpolation will be studied from a numerical point of view. The problems appearing there are of the following special structure. PROGRAM
PA. Minimize subject to
(i== 1, ...) n),
(nz,- , , mi) E W,
where Fi: R2 + R is a strictly convex function and Wi c R2 a closed convex set. For solving PA computationally here it is proposed to pass to the dual PROGRAM
DA. Maximize with
(1.2)
PO= P,!= 0.
The function H,*: R2 -+ i? is therein to be defined by H,*([,
~/)=sup{tx+qy-F,(x>
Y):
(x>
y)E
f+‘i).
(1.3)
In DA no constraints occur and, moreover, the gradient of the objective function is tridiagonal. These properties are indeed very convenient for 43 002 I -9045/88 $3.00 Copyright r 1988 by Academic Press. Inc All rights of reproduction in any farm reserved.
44
DIETZE AND SCHMIDT
numerical purposes and should be utilized by practical algorithms in order to be efficient. Attacking PA directly, e.g., by an active set strategy (see, e.g., Fletcher [Z]) is in general more expansive than starting from DA. Further, if a solution of DA is known then, under suitable assumptions, a solution of PA is explicitly given by the return-formula mr--I=alH*(Pi-l, mi=
-Pi)
a2HT(Pi-
12
(i = 1, ...) n).
-PiI
(1.4)
In addition, for problems arising in the shape preserving spline interpolation, the program DA can be explicitly stated becausethe subprograms (1.3) are solvable by formulas, and strong duality theorems are valid. Following this line the results of Burmeister et al. [l] concerning the convex spline interpolation are confirmed deductively. Moreover, also for the monotone spline interpolation problem, analogous results can be offered now that likewise are obtained in a deductive manner by applying the above-mentioned general concept. For further applications see [ 111, C131, C141, Cl51, and C161. 2. FORMULATION OF THE SPLINE PROBLEMS For describing the needed facts from the spline interpolation denote by d a partition of the unit interval [0, 11, A:x,=O
...
and let Sp(3, A) be the set of cubic Cl-splines on A. The points (x,, JJ,,),.... (x,, y,,) associated with A are said to be in convex position if T,
...
(2.1)
where si = ( yi - y, ,)/hi, hi = xi - x,+ , They are in monotone position if T, 2 0, ...) ?,,> 0.
(2.2)
Spline interpolants s E Sp(3, A) should satisfy the concrete interpolation condition 4x,) = Y,
(i = 0, ...) n).
(2.3)
Further, denote by m, = s’(x;) the first derivative of s in the node x,. Now, a result of Neuman [S] states that under (2.1) a spline interpolant s E Sp(3, A) is convex on [0, l] if and only if (m, I) mi) E W
(i= 1, .... n),
(2.4)
SHAPE
PRESERVING
SPLINE
INTERPOLANTS
45
where Wi = K, is the cone Ki={(X,y)ER2:2~+y~3~,,x+2y~3tj}.
(2.5)
An analogous result for the monotone spline interpolation, under the assumption (2.2), is due to Fritsch and Carlson [3]. Here Wj= Eiu D, is taken as the union of the ellipse E;= ((x, ,v)ER? x2+xy+y2-6q(x+y)+9tf
(2.6)
with the square D,= ((x, y)~R2:06x~32,,0~~~~35,}.
(2.7)
Naturally for a given set PPZ,,, .... m,, of first derivatives the spline interpolant s E Sp(3, A) is uniquely determined and can be directly computed. In the convex case the problem (2.4), with Wi= Ki, has in general no solution. An effective algorithm by which the solvability can be decided is given by Schmidt and Hess [7]; see also [9], [lo], and [12]. If at all, (2.4) has either one solution or an infinite number of solutions. In the monotone case the problem (2.4), with W, = E, u Di, is always solvable but in general not uniquely. Therefore, both in the convex and in the monotone spline interpolation case the question arises which solution of (2.4) i.e., which of the corresponding splines shall be chosen. A good strategy is to select that spline which has minimal mean curvature
I 5 s”(xJ2 dx = 0
i ;=,
i { rnf- 1+ mi_ ,m, + rnf - 3T,(m,-, + mi) + 3rf). I
X ?
FIG.
1.
The cone K.
(2.8)
46
DIETZE AND SCHMIDT
FIG. 2.
The set E,u D,
Thus one is led to a programming function Fj(X,
y)=R
I
X
37i
T
problem
of the type PA, with the
{X2+,uJ’+4.2-3ti(-\-+y)+32~)
(2.9)
and the sets Wi=K,
or
Wi=E,uDi
introduced by means of (2.5)-(2.7). The Hessian of the function (2.9) is positive implies
(2.10)
definite. This property
PROPOSITION 1. The program PA and the auxiliary programs (1.3) are uniquely solvable if speczfied by (2.9) and (2.10) and {f the feasible domains are nonempty. The optimal values are finite.
3. APPLICATIONS
OF FENCHEL’S
THEORY
In order to verify the duality of the above-mentioned pair PA and DA of programs and to derive duality statements one can refer to Fenchel’s theory and to extensions of it due to Rockafellar and Stoer, see, e.g., [4, 6, 81, but also Lagrange’s duality theory may be applied. Letf,g: R”-+a=Ru{+a, - cc } be proper convex functions, that is, f and g are convex functions greater than - cx, whose effective domains dom f and dom g are nonempty. Moreover, let ,f and g be closed. As usual,
47
SHAPE PRESERVING SPLINE INTERPOLANTS
the subdifferential off is denoted by af while the conjugate function will be called f*, f*(u)=sup{Uru-f(U)XER”‘}.
(3.1)
Further remember that for the indicator function I of a set W is I(u; W) = 0 for u E W and Z(u; W) = +ar, for u $ W. The subsequently used concepts of stability are defined, e.g., in [S]. Now, in Fenchel’s theory, to the primal convex program P: inf{ f( u) + g(u): u E R”’ }
(3.2)
belongs the dual concave program D: sup{ -f*(r)-g*(-u):
UER”j.
(3.3)
A corresponding strong duality theorem of a generality necessary here reads as follows (see [8, Section 5.31 for the first part and see,e.g., [6, Section 311 for the second part). THEOREM
2.
Let f and g be closed proper convex functions. Assume that
inf{f(u)+g(u):
UE R”} E R
(3.4)
and that there are vectors w, z E dom f n dom g such that f is w-stable and g is z-stable. Then, program D has an optimal solution, and the optimal values of P and D are equal:
inf{f(u)+g(u):
MER”} = maxi-f*(u)-g*(-u):oER”).
In addition, u and v are optimal for
(3.5)
P and D, respectively, jf and only if
u~df*(v)ndg*(-v)
(3.6)
holds.
Next, the programs PA and DA are shown to be special casesof P and D while the verification of the assumptions of the preceding theorem including the application of formula (3.6) will be undertaken later and only for the concrete spline problems. By introducing additional variables, Program PA can be reformulated as PA’: Minimize
i I=
Fi(mip , , Mi) 1
subject to (m,-,, M;)c
W,(i=
1, .... n), Mi=m,
while M, = m, is only a change of notation. 640’52’1.4
(i= 1, .... n- l),
48
DIETZE AND SCHMIDT
Naturally this is a separable program treated by several authors, see,e.g., [6, Section 311. Using the function Hi: R2 -+ R, Hib, y) = Fib, 4’) + 4x, y; W,), set, with u=(m,,M,,m
,,..., M,-,,m,-,, f(u)
=
f
i= I
H,(m,-
(3.7)
M,,)ER~“, I,
(3.8)
Mi).
Further let
g(u)=4%w,
W={u~R~“:M~=m~(i=l,...,n-1)).
(3.9)
By these definitions problem PA’ and therefore likewise PA become a program P. The occurring functions f and g are proper convex and closed. In order to state the corresponding dual program D the conjugate functions f* and g* to (3.8) and (3.9) are to be determined. Let u = (PO,P,,p,, .... Pnp,,pnm ,, P,*)E R’“. Then, in view of (3.8) it follows immediately that f*(U)=
i
,=I
H*(pi-
(3.10)
1, pi)
holds if H,* is given by (1.3). Observe further that g*(u) = sup{ urzX2.4 E W}
=o
for po=O,P,+pi=O(i=l
=+cc
otherwise
,..., n),p,,=O
(3.11)
and therefore g*( -u) = g*(u). Hence, problem D reduces to program DA. Now, the formula (3.6) shall be specified. To this end, it is additionally assumed that the auxiliary programs (1.3) are solvable for all ([, q) E R2 with finite optimal values HT([, q). Becauseof the strict convexity of F, the maximizer
in (1.3) is unique. Further, since (xi, yi) E i3H,+(<, q) if and only if (xi, yi) is a maximizer of (1.3) one gets the differentiability of HF and (3.12)
49
SHAPE PRESERVING SPLINE INTERPOLANTS
In addition, the derivatives a, H,? and & H,f+ are continuous. For this, see, e.g., [6, Sections 23 and 251. Thus, f* defined by (3.10) is continuously differentiable. Further, let v be optimal for D and let the optimal value be finite. This implies g*( -v) = 0, and the Kuhn-Tucker conditions for program D, i.e., for inf{f*(v):p,=O, P,+p,=O (i=l,..., n-l), P,,=O}, yield ?,H,*,,(p,, vector u by
P,+,)=8,H~(p,~,,
m,~ I=~lHiY~,-,,Pj),
P,), i= 1, .... n- 1. Now, define the
Mi=&H,%-
i = 1, .... n.
,, P,),
Then one gets g(u) = 0, v*u = 0, and hence g*( -v) + g(u) = -uTu. Thus it follows that UE 8g*( -u), see, e.g., [6, Section 231. Further, because of u E $f*(c), relation (3.6) is valid. Thus, the mentioned assumptions imply PROPOSITION 3. Let (m,, m, , .... m,,) with
(p,,, p,, .... p,,)
he
m,- l = 8, K+YP,-, , -pi) = fitpi-
a
solution
of program
DA
I) -PA
mi = azHI*(p;- 1) -P;) = YAPi- 13-P;) is the only solution continuous.
of’
(i= 1, .... n)
then
(3.13)
PA. The derivatives 2, H,*, 8, H,* are
4. CONVEX SPLINE INTERPOLATION As said above, for determining the convex spline intepolant with minimal mean curvature one has to solve the PROGRAM
PC. Minimize subject to
(m, _, , mi) c K,
(i= 1, .... n)
(4.1)
with Fi and Ki according to (2.9) and (2.5), respectively. By Proposition 1 this program is uniquely solvable if the feasible domain K= {(m,,m,,
.... m,)ER”+‘:
(mip,,mi)EK,(i=
is not empty. Note that K # @ implies (2.1).
l,..., n)}
(4.2)
50
DIETZE AND SCHMIDT
For numerical purposes the corresponding dual program DC shall be stated. To this end, the subprogram (1.3) now maximize d(x, u)=sx+qy-i
{(x-T)~+(x-~)(J--T)+(J~--)~) (4.3
subject to 2x + JJ< 32, x + 2~133r,
where h = h, and 7 = z,, is to be treated for all (t, q) E R2. In view of (1.21 and (1.4) especially, the optimal value d,,, = H,*([, q) is of interest. Proposition 3 assures that H,* is continuous differentiable. PROPOSITION 4.
The optimal value of program (4.3) is equal to
2 =r,(t+q)+k
6-i (
=7;(5+q)
)
.for q d 0, 21”6 ~1
.for 5 2 2q,21”3
FIG. 3. Pieces of the domain of function (4.4)
q.
(4.4)
SHAPE PRESERVING SPLINE INTERPOLANTS
51
The proof of (4.4) is elementary and is therefore only sketched. In the first case, the center of the ellipse d(x, y) = c, I=r+;(zc-a),
4$-w
(4.5)
is the maximizer of the problem (4.3). This is possible if and only if 5 ,< 0, q > 0, and the optimal value HT(t, q) = d(.Z, J) is as given in (4.4). In the second case, the maximizer is located on the line 2x +J = 3~, x < 5. This is valid only if ir > 0. The maximizer is the solution of 0,
Y) = c,
2x+y=3~,
where the constant c is to be chosen such that exactly one solution exists. One now gets
whenever X d z; that is, r - 2~ d 0. And HT(& q) = d(Z, j) with H,*(t, q) from (4.4) can be confirmed. The third case is analogous to the second. The maximizer is now lying on x + 2y = 32, y > z. The last case occurs for 4 - 2~ 2 0, 2t -q 2 0. Then the maximizer is in the vertex of the cone K,, X=j=r, and H*(5,11)=d(x,y)=t;(~+rl). Thus the proof is complete. Using Proposition 4 and (1.2) one is led to the following program DC which is dual to PC. PROGRAM
DC.
Maximize - i i=
I
HT(p,+ , , -Pi)
with
p. = pn = 0,
(4.6)
where H,? is given by (4.4). Next, by means of Theorem 2 the solvability of DC is shown. For this, let the feasible domain K of PC not be empty. Obviously, for (m,, m,, .... m,) E K it follows that z=(m,,m,,m,,
.... m,,-,,m,~,,m,)Edomfndomg
if, with Wi= K,,fand g are given by (3.8) and (3.9), respectively. Now, by stability statements obtainable from [IS, Sections 4.8, 5.3, and 5.51 the functions are proven to be z-stable. Indeed, F, considered as a function defined in R2” is z-stable, since z is an interior point of dom Fi. Further, in view of It-x, Y; K,) =1(x, y; K,,) + 4x, I’; K,,)
52
DIETZE AND SCHMIDT
with Kn={(x,y):2x+y<3r,}, K,j={(x,y):x+2y>3~,} and, since I( ., K,,), I( ., K,,) in spite of the affine constraints are polyhedral, the function Z( ., K,) is z-stable. Thus, this property holds also forf: In the same manner g is shown to be z-stable. Therefore, by the Theorem 2, program D possessesan optimal solution and thus likewise Program DC. Summarizing the main result of this section concerning the convex spline interpolation reads of follows. THEOREM 5. Both Programs PC and DC are solvable if and only if the feasible domain K of PC is nonempty, and then the optimal values are equal. Furthermore, if a solution (po,pI, .... p,) of DC is known, the solution (m,, ml, .... m,) of PC is given by formula (3.13) where HT is to he defined by (4.4).
As mentioned in the Introduction this theorem can be found already in paper [ 11. There it is derived in a more direct way be starting from the KuhnTucker conditions. Numerical tests presented in [l] show the strategy of solving PC via DC to be very efficient even if DC is handled by the ordinary Newton method; see also [ 111, [ 141, and [ 161.
5. MONOTONE SPLINE INTERPOLATION
As stated before, the monotone spline interpolation leads to PROGRAM PM.
Minimize
i
F,(m,- 134)
,=l
subject to (m,- , , mi) E Ei u D;
(i= 1, .... n),
(5.1)
where F;, E;, and Di are given by (2.9), (2.6), and (2.7). Here the feasible domain is always nonempty. Therefore, in view of Proposition 1, PM is always uniquely solvable but, as in the convex case, the solution can be computed more effectively by means of the corresponding dual program DM. For formulating this program the subprogram (1.3), now maximize d(x, y)=[x+qy-i{(x--T)2+(~--i)(y-T)+(y-r)2} (5.2) subject to (x, y)~ E u D E = E,, and D = D,, is to be considered for all with h=hi, T-T,, (5, v)eR*.
SHAPE
PRESERVING
SPLINE
53
INTERPOLANTS
PROPOSITION 6. The optimal value d,,, = HF(<, q) of program (5.2) is, with ai= 12zj/hi, given by
for (5,rl)EAiuBj for (5, v) E S,,
for (i", vIESI,
for (5,V)ETi +2a,Jt2
- 5yl+ q2 - ai(t + q) + af
for
(5,
VIE
ci
(5.3) where A,= {(5,rl)ER2:-ai~2r-)1~2a,,
-2a,<[-2qda,},
Bj= {(tl,~)ER2:52-5~+~2-ai(5+vl)~0}, Sli=
((5,rl)ER2:t-2v>a,,
S2;= {((,ty)ER2:2[-q6
-a,65da,}, -ai,
-a,Gq
Tj = {(L q) E R2: 5 < -ai, q 6 -a,},
C~=R'\(A~UB~US,~US,,U
T,).
%i Ti
FIG. 4.
Pieces of the domain
of function
(5.3).
54
DIETZE
AND
SCHMIDT
The proof of (5.3) is elementary. Only a sketch will be given. Compare also the proof of (4.4). In the first case the maximizer is equal to the center (4.5) of the ellipse d(x, y) = c, valid if and only if (t, q) eAiu Bi. In the second case the maximizer (X, j) is lying on ,v=O, 0 bx 6 35. Then one gets
if (5, yl) E S,;. The third case is symmetric to the second. In the fourth case the maximizer equals X = j = 0 true for (5, q) E Ti. The last case occurs if the maximizer is located on the “outer” boundary of the ellipse e(x, y)=i
{x~+x~+~~-~~(x+,v)+~~~}=O.
Then, the maximizer (X, j) solves d(x, y) = c, e(x, y) = 0 if c is chosen to be the greatest value such that exactly one solution exists. Hence, by adding the two equations one gets
with c1= 5 - 127/h, p = q - 12r/h, y = c - 24T2/h. Eliminating J in e(x, y) = 0, a quadratic equation in x arises. This has only one solution if y2 - 4r( LX+ p) y + 12&/l = 0. The greater root is
Now, because of d,,, = d,,,,, + emax= y + 24T2/h the value of H,*(t, q) given in (5.3) immediately results if (5, 9) E C,. Thus the proof is complete. By means of Proposition 6 and (1.2) the following program DM is obtained which is dual to PM.
PROGRAM DM. Maximize - i H,+(p, _ , , -P,) i= 1 where H,* is defined by (5.3).
with
p. = p,! = 0,
(5.4)
SHAPE PRESERVING SPLINE INTERPOLANTS
55
In order to verify the solvability of DM suppose that, in sharpening (2.21, T1 > 0, ...) Tn > 0
(5.5)
holds. Then, e.g., the vector z = (m,, M,, m,, .... M,.. ,, m,, M,) E R2”
with m,=M;=E/2, s=min{r ,, .... r,} > 0 is an interior point of dom H, even if the function H,, defined by (3.7) with W, = E, u D,, is considered as a function given on R2”. Thus F, is z-stable. In view of (3.8) one gets z ~domf and the z-stability off: Further, for the function g defined by (3.9) obviously z E dom g. And g is z-stable because g can be represented by a sum of polyhedral functions. Therefore, Theorem 2 assures the existence of an optimal solution for Program DM. Summarizing, the main result of this paper concerning the monotone spline interpolation is the following. THEOREM 7.
Both programs PM and DM possess optimal solutions if (5.5) holds, and then the optimal values are equal. Moreover, if a solution (PO,PI, ...?p,) of DM is given, the solution (m,, ml, .... m,) of PM can he computed by means of formula (3.13), where HI is to be defined by (5.3).
6. PROGRAMS WITH GIVEN BOUNDARY VALUES In the monotone spline interpolation ti = 0
for some jE { 1, .... n}
(6.1)
may occur. Then, mjp, = mj= 0 for the solution of PM follows, and the optimization problem PM can be decomposed into subproblems of lower dimension. Analogously, if in the convex spline interpolation 5j=T,+
I
for someJo ( 1, .... n - 1 }
(6.2)
appears one gets m, ~ 1= mj = mi + , = TV, and Program PC divides into lower dimensional subprograms. Of course, from a numerical point of view such a decomposition should be undertaken if possible. In the resulting subproblems, one or both of the boundary values, say m, and m,, are prescribed. For example, let m, = z and m, = j be given. Then one has to consider the slightly modified
56
DIETZE AND SCHMIDT
PR~CRAM
PB.
Minimize
i
F,(m,- , , mi) I subject to (m,+, , mi) E Wj(i = 1, .... n),
(6.3)
i=
mo= M,m,,= fl.
Two ways are offered to derive a corresponding dual program. In the first case, define the functions fand g by (3.8) and (3.9) and let the set W now be W= {zLER~~: Mi=mj(i=l,...,n),m,=rx,m,=~}.
(6.4)
Then, it is immediately seen that for ~~+P,=O(i=l,...,n)
g*(u) = vo + PPll = +m
otherwise.
(6.5)
Using the functions HT given by (1.3) the dual program now obtained reads as follows. PROGRAM
DB,. Maximize - f
H,?(p,
,, -P,) + wo - BP,,.
(6.6)
r=l
In a second case, the variables m, and m, are eliminated, and thus PB is considered to be a problem in the space RZnp2. Following the line described in Section 3 one is led to the dual PROGRAM
DB,.
Maximize -G:( -p,)-
II- I 1 NT@,- ,, -pi)-G,*(p,_,).
(6.7)
,=2
Here H,* are defined by (1.3) and G:, G,* by G:(v) = sup{~ -P,(a, Y): (~9Y) E W, 1,
(6.8)
G,*(O=sup{tx-F,(x,
(6.9)
8): (x, B)E W,}.
In the shape preserving spline interpolation, if (6.1) or (6.2) occur, it is recommended to decompose the program for determining the interpolant with minimal curvature into lower dimensional subprograms. These are programs with boundary conditions of the type PB, to be solved via DB, or DB,. For the convex spline problem, Theorem 5 holds verbatim if PB is
SHAPE PRESERVING SPLINE INTERPOLANTS
51
coupled with DB, or with DB,, while for the monotone spline problem, an analog to Theorem 7 can be assured if PB is coupled with DB,. The details are omitted because their verification is now straightforward.
REFERENCES I. W. BURMEISTER,W. HESS, AND J. W. SCHMIDT, Convex spline interpolants with minimal curvature, Computing 35 (1985), 219-229. 2. R. FLETCHER, “Practical Methods of Optimization. Vol. 2. Constrained Optimization,” Wiley, New York/Chichester/Brisbane/Toronto/Singapore, 1981. 3. F. N. FKKSCH AND R. E. CARLSON, Monotone piecewise cubic interpolation, SIAM J. Numer. Anal. 17 (1980), 238-246. 4. P.-J. LAURENT, “Approximation et Optimisation,” Hermann, Paris, 1972. 5. E. NEUMAN, Uniform approximation by some Hermite interpolating splines, J. Compuf. Appl. Math. 4 (1978), 7-9. 6. R. T. ROCKAFELLAR, “Convex Analysis,” Princeton Univ. Press, Princeton, NJ, 1970. 7. J. W. SCHMIDT UND W. HESS, Schwach verkoppelte Ungleichungssysteme und konvexe Spline-Interpolation, E/em. Math. 39 (1984), 85-95. and Optimization in Finite Dimensions,” 8. J. STOER AND C. WITZGALL, “Convexity Springer-Verlag, Berlin/Heidelberg/New York, 1970. 9. P. COSTANTINI, On monotone and convex spline interpolation, Mafh. Compuf. 46 (1986), 203-214. 10. P. COSTAN~IN~ AND R. MORANDI, Monotone and convex cubic spline interpolation, Calcolo 21 (1984), 281-294. 1 I. W. HESSAND J. W. SCHMIDT, Convexity preserving interpolation with exponential splines, Computing 36 (1986), 335-342. 12. J. W. SCHMIDT, Convex interval interpolation with cubic splines, BIT 26 (1986), 377-387. 13. J. W. SCHMIDT, On shape preserving spline interpolation: Existence theorems and determination of optimal splines, Banach Center Publ., Semester XXVII, Polish Scient. Publishers, Warsaw, in press. 14. J. W. SCHMIDT, On tridiagonal linear complementarity problems, Numer. Mafh. 51 (1987), 11-21. 15. J. W. SCHMIDT, An unconstrained dual program for computing convex C’-spline approximants, Compufinz, in press. 16. J. W. SCHMIDT AND W. HESS, Positivity of polynomials on intervals and positive spline interpolation, submitted.