# A Spline Smoothing Newton Method forL∞ Distance Regression with Bound Constraints

Authors: Li Dong; Bo Yu
Journal: ISRN Operations Research
 (2013)

---

## Abstract

Orthogonal distance regression is arguably the most common criterion for fitting a model to data with errors in the observations. It is not appropriate to force the distances to be orthogonal, when angular information is available about the measured data points. We consider here a natural generalization of a particular formulation of that problem which involves the replacement ofl2 norm by l∞ norm. This criterion may be a more appropriate one in the context of accept/reject decisions for manufacture parts. For l∞ distance regression with bound constraints, we give a smoothing Newton method which uses cubic spline and aggregate function, to smooth max function. The main spline smoothing technique uses a smooth cubic spline instead of max function and only few components in the max function are computed; hence it acts also as an active set technique, so it is more efficient for the problem with large amounts of measured data. Numerical tests in comparison to some other methods show that the new method is very efficient.


---

## Body

## 1. Introduction

For fitting curves or surfaces to observed or measured data, a common criterion is orthogonal distance regression (ODR). Given the data pairs(uj,wj), j=1,…,m, where uj∈Rp1 is the independent variable and wj∈R1 is the dependent variable; suppose
(1)wj=g(d,uj),
where d∈Rp2 is a vector of parameters to be determined. We assumed that δj is the random error associated with uj and that εj is the random error associated with wj. To be precise, we relate the quantities uj, wj, δj, and εj to
(2)wj=g(d,uj+δj)+εj,j=1,2,…,m.
As shown in Boggs et al. [1] this gives rise to the ODR problem given by
(3)mind,δ12∑j=1m‍[(wj-g(d,uj+δj))2+δj2].
The ODR problem can be solved by the Gauss-Newton or Levenberg-Marquardt methods (see [1, 2]). The general form of the bounded constrained ODR problem can be expressed by
(4)mind,δ12∑j=1m‍[(wj-g(d,uj+δj))2+δj2],s.t.Ld≤d≤Ud,
where Ld and Ud are vectors of length p2 that provide the lower and upper bounds on d, respectively. Zwolak et al. give the algorithm to handle (4) in [3].It is not appropriate to force the distances to be orthogonal, when angular information is available about the measured data points, such as the rotated cone fitting problem and rotated paraboloid fitting problem. Then, (4) becomes
(5)minx,δ12∑j=1m‍[(w-j(x)-g(x,δj))2+δj2],s.t.Lx≤x≤Ux,
where (u-j,w-j)=((uj+a),(wj+b))T(μ),  T(μ) is an orthogonal matrix, a∈Rp1, b∈R1, μ∈Rp3, x=(a,b,μ,d), and Lx and Ux are vectors of length p=p1+1+p3+p2.When the least squares norm is not appropriate, problem (5) can be generalized to use other measures in a variety of ways. Most generalizations have been of formulations (5), with the l2 norms replaced with other norms. We consider here l∞ norms. It may be a more appropriate one in the context of accept/reject decisions for manufacture parts (see [4]). In this paper, we consider the following l∞ distance regression with bound constraints.Let(6)f-j(x,δ)={w-j(x)-g(x,δj)j=1,…,mδj-mj=m+1,…,2m,(7)minx,δf(x,δ)=max1≤j≤2m|f-j(x,δ)|,s.t.x∈Λ={x∣Lx≤x≤Ux}.We know (7) is a minimax problem. There are several different algorithms that have been taken to solve (7), such as subgradient methods (see [5]), SQP methods (see [6–8]), bundle-type methods (see [9–12]), and smooth approximation methods (see [13–21]). For (7), |f-i(x)| is nonsmooth function including the absolute value function. Moreover, when large amounts of measured data are to be fitted to a model, the number of components in the maximum function is very large. It is necessary to develop efficient solution methods for problem (7).In this paper, we consider to uniformly approximatef-(x,δ) by the smooth splines introduced in [22].Let us first recall the formulation of multivariate splines. LetD be a polyhedral domain of Rm which is partitioned with irreducible algebraic surfaces into cells Δ={Δi∣i=1,…,N}. A function s(z) defined on D is called a k-spline function with rth order smoothness, expressed for short as s(z)∈Skr(D,Δ), if s(z)∈Cr(D) and s(z)∣Δi=pi∈Pk, where Pk is the set of all polynomials of degree k or less in m variables. Similar to the smooth splines which uniformly approximate min{z1,z2,…,zm} given in [22], we can construct a spline function s32(z;ε)∈S32(Rm,ΔMS2) to uniformly approximate max{z1,z2,…,zm} (as ε→+0), where ΔMS2 is the homogenous Morgan-Scott partition of type two in [22], as follows:
(8)s32(z1,z2,…,zm;ε)=zi1+∑l=1k-1‍cl(lzil+1-∑j=1l‍zij+ε)3,forz∈Δi1⋯ik(ε),
where c1=1/(6ε2), ck/ck+1=(k+2)/k, 1≤k≤m, and the cell Δi1⋯ik(ε) is the region defined by the following inequalities:
(9)zil-zil+1≥0,when1≤l<k,(k-1)zik-∑j=1k-1‍zij+ε≥0,kzil-∑j=1k‍zij+ε≤0,whenk+1≤l≤m.
The spline smoothing technique uses a smooth cubic spline instead of max function, and only few components in the max function are computed; hence it acts also as an active set technique, so it is more efficient for the minimax problems with nonsmoothness and large numbers of components.For that (7) is a minimax problem with bound constraints, and we cannot utilize SSN algorithm in [23] directly to solve (7). Here, we try to extend the idea of SSN algorithm to solve it. At first, we use penalty function to transform (7) into an unconstrained minimax problem. Then, using the smooth approximation, a smoothing Newton method (SN) can be used to solve the l∞ distance regression with bound constraints.


## 2. The SN Algorithm forl∞ Distance Regression with Bound Constraints

Firstly, some deformations for (7) are necessary. Due to |f-j(x,δ)|=max{f-j(x,δ),-f-j(x,δ)}, then (7) is equivalent to
(10)minx,δf(x,δ)=max1≤j≤4m{fj(x,δ)},s.t.x∈Λ,
where
(11)fj(x,δ)={f-j(x,δ)j=1,…,2m-f-j-2m(x,δ)j=2m+1,…,4m.Assumption 1.
We assume that the functionsfj(x,δ), j=1,…,4m, are twice continuously differentiable.Letϕ(x)=max{ϕ1(x),…,ϕ2p(x)}, where ϕ1(x)=x1-Ux,1,…, and ϕp(x)=xp-Ux,p, ϕp+1(x)=-x1-Lx,1,…, ϕ2p(x)=-xp-Lx,p. Denote the unknown variables x, δ to be x; then x∈Rn, where n=p+p1. Use penalty function with penalty parameter C>0 to transform (10) into the following unconstrained programming problem:
(12)min{ψ(x)={max1≤j≤4mfj(x)+Cϕ(x)}}.The following proposition concerning Theorem4.2.8 in [24] gives the first-order optimality condition for (12).Proposition 2.
Suppose that Assumption1 holds, and then if (12) attains the extremum at x*, then
(13)0∈∂ψ(x*)≜convj∈p(x*){∂fj(x*)}+Cconvi∈q(x*){∂ϕi(x*)},
where s(x*)={j∈s={1,…,4m}∣fj(x*)=f(x*)}, q(x*)={i∈q={1,…,2p}∣ϕi(x*)=ϕ(x*)}.We use the following cubic spline to smoothf(x) and the aggregate function to smooth ϕ(x) in (12).Consider(14)Ψt(x)=Ft(x)+CΦt(x),
where
(15)Ft(x)=s32(f1,t(x),f2,t(x),…,f4m,t(x);t),Φt(x)=tln(∑i=12p‍exp(ϕi(x)t)).Remark 3.
Under Assumption1, Ft(x) and Φt(x) are twice continuously differentiable for arbitrary t>0.From Lemma  1.1 in [23], Proposition  3.3 in [25], and Proposition  2.4 in [26], we have the following proposition.Proposition 4.
(
1
) For any x∈Λ, Ft(x) and Φt(x) are monotonically increasing with respect to t>0.
(
2
) Suppose that Assumption 1 holds. Then, for any t>0, Ψt(x) is twice continuously differentiable, and
(16)∇Ft(x)=∑i=14m‍λi,t(x)∇fi,t(x)=∑j=1k‍λtij(x)∇fij,t(x),
where
(17)λtij(x)={1-3∑l=1k-1‍cl(hl(x,t))2j=13(j-1)cj-1(hj-1(x,t))2-3∑l=jk-1‍cl(hl(x,t))22≤j<k3(k-1)ck-1(hk-1(x,t))2j=k0k<j≤4m,
and hl(x,t)=lfil+1,t(x)-∑r=1lfir,t(x)+t,
(18)∇Φt(x)=∑i∈q‍ξi,t(x)∇ϕi(x),
where
(19)ξi,t(x)=exp(ϕi(x)/t)∑i∈q‍exp(ϕi(x)/t)∈(0,1],∑i∈q‍ξi,t(x)=1.Lemma 5.
Suppose that Assumption1 holds, Then, for every bounded set S⊂Rn, there exists an L<∞ such that
(20)〈y,∇2Ψt(x)y〉≤1tL∥y∥2,
for all x∈S, y∈Rn, and 0<t≤1.Proof.
From Proposition4, Lemma  1.2 in [23], and Lemma  2.2 in [17], we know that for every bounded set S⊂Rn, there exist LF<∞ and LΦ<∞ such that 〈y,∇2Ft(x)y〉≤(1/t)LF∥y∥2 and 〈y,∇2Φt(x)y〉≤(1/t)LΦ∥y∥2 for all x∈S, y∈Rn, and 0<t≤1. We also know ∇2Ψt(x)=∇2Ft(x)+C∇2Φt(x). Let L=LF+CLΦ; then 〈y,∇2Ψt(x)y〉≤(1/t)L∥y∥2.Algorithm 6 (The smoothing Newton algorithm).
S
t
e
p
0. x0∈Λ; t0>0, t^≫1, and α, β, and κ1∈(0,1); κ2≫1, 0<κ3≪1 and δ>0; functions ϵa(t), ϵb(t) and τ(t):(0,∞)→(0,∞), satisfying ϵb(t)≥ϵa(t)>τ(t) for all t>0. Set i=0, k=0, s=1, and xk,i=x0.
Search the CellS
t
e
p
1. Let I-={j∣max{f1,tk(x(k,i)),…,f4m,tk(x(k,i))}-fj,tk(x(k,i))<tk}; k- is the cardinality of I-, and I-={i1,i2,…,ik-}; range {fij,tk(x(k,i))}j=1k- according to fi1,tk(x(k,i))≥fi2,tk(x(k,i))≥⋯≥fik-,tk(x(k,i)). If k-=1, the cell is Δi1(tk). Else, for every k~∈{k-,k--1,…,2}, if (k~-1)fik~,tk(x(k,i))-∑j=1k~-1‍fij,tk(x(k,i))+tk≥0, we have k~∈I⊆{k-,k--1,…,2}. Let k^ be the maximum element of I, then the cell is Δi1⋯ik^(tk).The Stabilized Newton-Armijo AlgorithmS
t
e
p
2. Go to Step 1, and compute ∇Ψtk(xk,i)=∇Ft(xk,i)+C∇Φt(xk,i). If ∥∇Ψtk(xk,i)∥2>τ(tk), go to Step 3. Else go to Step 8.S
t
e
p
3. Compute Btk(xk,i)(21)Btk(xk,i)=θ(xk,i)I+∇2Ψtk(xk,i),
where θ(x)=max{0,δ-e(x)} with e(x) denoting the minimum eigenvalue of ∇2Ψtk(x); then compute the Cholesky factor R such that Btk(xk,i)=RRT and the reciprocal condition number c(R) of R. If c(R)≥κ1 and pk≥κ3, go to Step 4. Else, if c(R)≥κ1 and the largest eigenvalue σpk,max(xk,i) of Btk(xk,i) satisfies σtk,max(xk,i)≤κ2, go to Step 4; else go to Step 5.S
t
e
p
4. Compute the search direction
(22)hk,i=-Btk(xk,i)-1∇Ψtk(xk,i),
and go to Step 6.S
t
e
p
5. Compute the search direction
(23)hk,i=-∇Ψtk(xk,i).S
t
e
p
6. Compute the step length λk,i=βl, where l≥0 is the smallest integer satisfying
(24)Ψtk(xk,i+βlhk,i)-Ψtk(xk,i)≤αβl〈∇Ψtk(xk,i),hk,i〉.S
t
e
p
7. Set xk,i+1=xk,i+λk,ihk,i, i=i+1. Go to Step 1, and compute ∇Ψtk(xk,i). If
(25)∥∇Ψtk(xk,i)∥2≤τ(tk),
go to Step 8; else go to Step 3.Adaptive Parameter DecreaseS
t
e
p
8. If s=1, compute t* such that
(26)ϵa(tk)≤∥∇Ψt*(xk,i)∥2≤ϵb(tk),
go to Step 9, else set tk+1=1/s(k+2), k=k+1, and i=0, and go to Step 2.S
t
e
p
9. If t*≥t^, set tk+1=min{t*,tk/(tk+1)}, k=k+1, and i=0, and go to Step 2; else set s=max{2,((1/t^)+2)/(k+1)}, tk+1=1/s(k+2), k=k+1, and i=0, and go to Step 2.Lemma 7.
Suppose that Assumption1 holds and that sequences {tk} and {x1,i}, {x2,i},…,{xk,i},…, are generated by Algorithm 6. Then the following properties hold: (1) the sequences {tk} is decreasing; (2) If {x1,i},{x2,i},…,{xk,i},… has an accumulation point, then tk→0 and ∑k=0∞tk=+∞.Lemma 8.
Suppose that Assumption1 holds. Then, for every bounded set S⊂Rn and parameters α,β∈(0,1), there exists a Ks<∞ such that, for all 0<t≤1 and x∈S,
(27)Ψt(x+λt(x)ht(x))-Ψt(x)≤-αKs∥∇Ψt(x)∥2t,
where λt(x) is the stepsize of Algorithm 6 (see (24)).Lemma 9.
Suppose that Assumption1 holds and that {x1,i},{x2,i},…,{xk,i},… is a bounded sequence generated by Algorithm 6. Then, for any k, the sequence {xk,i} is finite; that is, there exists a ik∈N such that (25) holds for i=ik.The proofs of Lemmas7–9 are similar to that in [23] and omitted here.Theorem 10.
Suppose that Assumption1 holds and that {xk,ik}k=0∞ is a bounded sequence generated by Algorithm 6. Then, there exists an infinite subset K⊂N and a x^∈Rn such that xk,ik→Kx^ and 0∈∂ψ(x^).Proof.
Suppose that{xk,ik}k=0∞ is a bounded sequence generated by Algorithm 6. For the sake of a contradiction, suppose that there exists an ε>0 such that
(28)limk→∞inf∥∇Ψtk(xk,ik)∥≥ε.
Since {xk,ik}k=0∞ is a bounded sequence, it has at least one accumulation point. Hence, by Lemma 7, tk→0 as k→∞. Next, by Lemma 8, there exists an M<∞ such that
(29)Ψtk(xk,i+1)-Ψtk(xk,i)≤-αM∥∇Ψtk(xk,i)∥2tk
for all k∈N. Hence, for all k∈N,
(30)Ψtk+1(xk+1,ik+1)-Ψtk(xk,ik)=(Ψtk+1(xk+1,ik+1)-Ψtk+1(xk+1,ik+1-1))+(Ψtk+1(xk+1,ik+1-1)-Ψtk+1(xk+1,ik+1-2))+⋯+(Ψtk+1(xk+1,1)-Ψtk+1(xk+1,0))+(Ψtk+1(xk+1,0)-Ψtk(xk,ik))≤-αM∑j=0ik+1-1‍∥∇Ψtk+1(xk+1,j)∥2tk+1,
where we have used the fact from Proposition 4 that Ψtk+1(xk+1,0)-Ψtk(xk,ik)≤0 for all k∈N. In view of Algorithm 6, we know xk+1,0=xk,ik. Then,
(31)Ψtk+1(xk+1,ik+1)=(Ψtk+1(xk+1,ik+1)-Ψtk(xk,ik))+(Ψtk(xk,ik)-Ψtk-1(xk-1,ik-1))+⋯+(Ψt1(x1,i1)-Ψt0(x0,i0))+Ψt0(x0,i0)≤-αM∑j=0ik+1-1‍∥∇Ψtk+1(xk+1,j)∥2tk+1-αM∑j=0ik-1‍∥∇Ψtk(xk,j)∥2tk-⋯-αM∑j=0i0-1‍∥∇Ψt0(x0,j)∥2t0+Ψt0(x0,i0).
By Lemma 7, ∑k=0∞tk=+∞. It follows from (28) and (31) that Ψtk(xk,ik)→-∞, as k→∞. But for every accumulation point x* of {xk,ik}, that is, there exists an infinite subset K*⊂N such that xk,ik→K*x*, we have by continuity Ψtk(xk,ik)→K*Ψ(x*), which is a contradiction. Hence,
(32)limk→∞inf∥∇Ψtk(xk,ik)∥=0.
Now by (16) and (18) and Proposition  3.2 in [25], we have that, for all x∈Rn,
(33)limk→∞∇Ψt(xk,ik)=∑j∈s(x^)‍λ^j(x^)∇fj(x^)+C∑i∈q(x^)‍ξ^i(x^)∇ϕi(x^)=0.
According to Proposition  3.2 in [25], we have λ^j(x^),ξ^i(x^)≥0, j∈p(x^), and i∈q(x^) such that ∑j∈p(x^)λ^j(x^)=1, ∑i∈q(x^)ξ^i(x^)=1. Hence by (32) and (33), there has to exist an infinite subset K⊂N and x^∈Rn such that xk,ik→Kx^ and 0∈∂ψ(x^). This completes the proof.


## 3. Numerical Experiment

In this section, we consider the rotated cone fitting problem in [27]. In this example, the error δj associated with uj is zero, and uj=(uj1,uj2). Let (uj1,uj2,andwj)∈R3, i=1,…,m, be the measured data. Define the orthogonal matrix
(34)T(x4,x5)=[1000cos(x4)sin(x4)0-sin(x4)cos(x4)]×[cos(x5)0-sin(x5)010sin(x5)0cos(x5)].
Then (uj1,uj2,andwj) is transformed to (u-j1,u-j2,andw-j)=((uj1+x1),(uj2+x2),(wj+x3))T(x4,x5). The fitting cone V is defined as w-=x6(u-j1)2+(u-j2)2, x6≥0. Denote the unknown variables to be x=(x1,x2,x3,x4,x5,andx6); then the fitting problem is equivalent to
(35)minf(x)=max1≤j≤2m{fj(x)},s.t.x∈Λ={x∣-π≤x4≤π,-π≤x5≤π,0≤x6},
where
(36)fj(x)={w-j(x)-x6(u-j1(x))2+(u-j2(x))2j=1,…,m-fj-m(x)j=m+1,…,2m.Sincefj(x) in (35) is nonsmooth in ℜj={x∣(u-j1(x))2+((u-j2)(x))2=0}, we try to smooth it by the following function:
(37)fj,t(x)={w-j(x)-x6(u-j1(x))2+(u-j2(x))2+t+x6t,aaaaaaaaaaaaaaaaaaaaaaaaaaj=1,…,m,x6(u-j1(x))2+(u-j2(x))2+t-w-j(x),aaaaaaaaaaaaaaaaj=m+1,…,2m.Letϕ(x)=max{ϕ1(x),…,ϕ6(x)}, where ϕ1(x)=x4-π, ϕ2(x)=-x4-π, ϕ3(x)=x5-π, ϕ4(x)=-x5-π, ϕ5(x)=-x6, ϕ6(x)=0.Use penalty function with penalty parameterC>0 to transform (35) into the following unconstrained programming problem:
(38)min{ψ(x)={max1≤j≤2mfj,t(x)+Cϕ(x)}}.
According to the definition of fj,t(x), it is easy to obtain the following proposition.Proposition 11.
For anyt>0, fj,t(x) defined as (37) is twice continuously differentiable, and for any given x∈Λ, fj,t(x) is monotonically increasing with respect to t>0, and fj,t(x)→fj(x) as t→0.Remark 12.
Under Proposition11, Ft(x) is monotonically increasing with respect to t>0.Remark 13.
Suppose that{xk,ik}k=0∞ is a bounded sequence generated by Algorithm 6. Then, there exists an infinite subset K⊂N and a x^∈R6 such that xk,ik→Kx^ and
(39)0∈convj∈s(x^){∂fj(x^)}+Cconvi∈q(x^){∂ϕi(x^)}.
Moreover, if x^∉ℜ=⋃1≤j≤2m‍ℜj, then 0∈∂ψ(x^).We have implemented the SN algorithm using the MATLAB for problem (38). In order to show the efficiency of the algorithm, we also have implemented TSN algorithms using similar procedures and an SQP algorithm that is implemented by calling MATLAB function fminimax directly. Algorithm TSN was proposed by Xiao et al. in [27].The test results were obtained by running MATLAB R2011a on a desktop with Windows XP Professional operation system, Intel(R) Core(TM) i3-370 2.40 GHz processor and 2.92 GB of memory. The default parameters are chosen as follows:(40)α=0.8,β=0.77,t^=105ln(2m),κ1=10-7,κ2=1030,κ3=1000t^,τ(t)=10-3,t0=1,(ϵa,ϵb)=(0.01,0.2),δ=0.1,C=100.The results are listed in Table1, where x* denotes the final approximate solution point and f(x*) is the value of the objective function at x*. Time is the CPU time in seconds.Table 1
Test result for the example,x0=(0,0,-2,2,-2,1), t=10-5.
m
Method
x
*
f
(
x
*
)
Time
SN
(
-
0.156746
,
-
0.124595
,
…
,
1.732075
)
0.744533
1.891536
20000
TSN
(
-
0.156611
,
-
0.123928
,
…
,
1.736362
)
0.743643
3.878818
SQP
(
-
0.155824
,
-
0.126803
,
…
,
1.760751
)
0.927184
31.302830
SN
(
-
0.156984
,
-
0.125677
,
…
,
1.732449
)
0.851957
4.166439
50000
TSN
(
-
0.155253
,
-
0.123374
,
…
,
1.733202
)
0.851803
10.761201
SQP
(
-
0.155824
,
-
0.126194
,
…
,
1.760880
)
0.971305
33.716271
SN
(
-
0.155158
,
-
0.124954
,
…
,
1.736343
)
0.856780
16.437595
100000
TSN
(
-
0.155159
,
-
0.124956
,
…
,
1.736346
)
0.856767
37.235130
SQP
(
-
0.157418
,
-
0.125016
,
…
,
1.761221
)
1.025103
101.817872We test our algorithm for the artificial rotated cone data points which are generated as that in [28]. At first, produce points {(z~1j,z~2j,andz~3j)}j=1m on an unrotated cone by defining the z~1j=rjtan(π/6)cos(γj), z~2j=rjtan(π/6)sin(rj), z~3j=rj, where rj and γj are equally distributed pseudorandom numbers in [1;10] and [0;2π], respectively. Then, perturb z~3j by adding error item which follows the Gaussian distribution N(0;0.3), and make rotations and translation to obtain the final data ((uj1+2.1),(uj2-1.4),(wj+1.3))=(z~1j,z~2j,z~3j)T-1(π/20,π/25).



---

*Source: 393482-2013-04-16.xml*