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 of norm by norm. This criterion may be a more appropriate one in the context of accept/reject decisions for manufacture parts. For 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.

1. Introduction

For fitting curves or surfaces to observed or measured data, a common criterion is orthogonal distance regression (ODR). Given the data pairs , , where is the independent variable and is the dependent variable; suppose where is a vector of parameters to be determined. We assumed that is the random error associated with and that is the random error associated with . To be precise, we relate the quantities , , , and to As shown in Boggs et al. [1] this gives rise to the ODR problem given by 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 where and are vectors of length that provide the lower and upper bounds on , 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 where ,   is an orthogonal matrix, , , , , and and are vectors of length .

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 norms replaced with other norms. We consider here 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 distance regression with bound constraints.

Let

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 [68]), bundle-type methods (see [912]), and smooth approximation methods (see [1321]). For (7), 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 approximate by the smooth splines introduced in [22].

Let us first recall the formulation of multivariate splines. Let be a polyhedral domain of which is partitioned with irreducible algebraic surfaces into cells . A function defined on is called a -spline function with th order smoothness, expressed for short as , if and , where is the set of all polynomials of degree or less in variables. Similar to the smooth splines which uniformly approximate given in [22], we can construct a spline function to uniformly approximate (as ), where is the homogenous Morgan-Scott partition of type two in [22], as follows: where , , , and the cell is the region defined by the following inequalities: 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 distance regression with bound constraints.

2. The SN Algorithm for Distance Regression with Bound Constraints

Firstly, some deformations for (7) are necessary. Due to , then (7) is equivalent to where

Assumption 1. We assume that the functions , , are twice continuously differentiable.

Let , where , and , , . Denote the unknown variables , to be ; then , where . Use penalty function with penalty parameter to transform (10) into the following unconstrained programming problem:

The following proposition concerning Theorem   in [24] gives the first-order optimality condition for (12).

Proposition 2. Suppose that Assumption 1 holds, and then if (12) attains the extremum at , then where , .

We use the following cubic spline to smooth and the aggregate function to smooth in (12).

Consider where

Remark 3. Under Assumption 1, and are twice continuously differentiable for arbitrary .

From Lemma  1.1 in [23], Proposition  3.3 in [25], and Proposition  2.4 in [26], we have the following proposition.

Proposition 4. For any , and are monotonically increasing with respect to .
Suppose that Assumption 1 holds. Then, for any , is twice continuously differentiable, and where and , where

Lemma 5. Suppose that Assumption 1 holds, Then, for every bounded set , there exists an such that for all , , and .

Proof. From Proposition 4, Lemma  1.2 in [23], and Lemma  2.2 in [17], we know that for every bounded set , there exist and such that and for all , , and . We also know . Let ; then .

Algorithm 6 (The smoothing Newton algorithm).
. ; , , and , , and ; , and ; functions , and , satisfying for all . Set , , , and .
Search the Cell

. Let ; is the cardinality of , and ; range according to . If , the cell is . Else, for every , if , we have . Let be the maximum element of , then the cell is .

The Stabilized Newton-Armijo Algorithm

. Go to Step 1, and compute . If , go to Step 3. Else go to Step 8.

. Compute where with denoting the minimum eigenvalue of ; then compute the Cholesky factor such that and the reciprocal condition number of . If and , go to Step 4. Else, if and the largest eigenvalue of satisfies , go to Step 4; else go to Step 5.

. Compute the search direction and go to Step 6.

. Compute the search direction

. Compute the step length , where is the smallest integer satisfying

. Set , . Go to Step 1, and compute . If go to Step 8; else go to Step 3.

Adaptive Parameter Decrease

. If , compute such that go to Step 9, else set , , and , and go to Step 2.

. If , set , , and , and go to Step 2; else set , , , and , and go to Step 2.

Lemma 7. Suppose that Assumption 1 holds and that sequences and , , are generated by Algorithm 6. Then the following properties hold: the sequences is decreasing; If has an accumulation point, then and .

Lemma 8. Suppose that Assumption 1 holds. Then, for every bounded set and parameters , there exists a such that, for all and , where is the stepsize of Algorithm 6 (see (24)).

Lemma 9. Suppose that Assumption 1 holds and that is a bounded sequence generated by Algorithm 6. Then, for any , the sequence is finite; that is, there exists a such that (25) holds for .

The proofs of Lemmas 79 are similar to that in [23] and omitted here.

Theorem 10. Suppose that Assumption 1 holds and that is a bounded sequence generated by Algorithm 6. Then, there exists an infinite subset and a such that and .

Proof. Suppose that is a bounded sequence generated by Algorithm 6. For the sake of a contradiction, suppose that there exists an such that Since is a bounded sequence, it has at least one accumulation point. Hence, by Lemma 7, as . Next, by Lemma 8, there exists an such that for all . Hence, for all , where we have used the fact from Proposition 4 that for all . In view of Algorithm 6, we know . Then, By Lemma 7, . It follows from (28) and (31) that , as . But for every accumulation point of , that is, there exists an infinite subset such that , we have by continuity , which is a contradiction. Hence, Now by (16) and (18) and Proposition  3.2 in [25], we have that, for all , According to Proposition  3.2 in [25], we have , , and such that , . Hence by (32) and (33), there has to exist an infinite subset and such that and . This completes the proof.

3. Numerical Experiment

In this section, we consider the rotated cone fitting problem in [27]. In this example, the error associated with is zero, and . Let , , be the measured data. Define the orthogonal matrix Then is transformed to . The fitting cone is defined as , . Denote the unknown variables to be ; then the fitting problem is equivalent to where

Since in (35) is nonsmooth in , we try to smooth it by the following function:

Let , where , , , , , .

Use penalty function with penalty parameter to transform (35) into the following unconstrained programming problem: According to the definition of , it is easy to obtain the following proposition.

Proposition 11. For any , defined as (37) is twice continuously differentiable, and for any given , is monotonically increasing with respect to , and as .

Remark 12. Under Proposition 11, is monotonically increasing with respect to .

Remark 13. Suppose that is a bounded sequence generated by Algorithm 6. Then, there exists an infinite subset and a such that and Moreover, if , then .

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:

The results are listed in Table 1, where denotes the final approximate solution point and is the value of the objective function at . Time is the CPU time in seconds.

We test our algorithm for the artificial rotated cone data points which are generated as that in [28]. At first, produce points on an unrotated cone by defining the , , , where and are equally distributed pseudorandom numbers in and , respectively. Then, perturb by adding error item which follows the Gaussian distribution , and make rotations and translation to obtain the final data .

Acknowledgments

This work was supported by the National Natural Science Foundation of China (11171051, 91230103, and 61153001) and the Fundamental Research Funds for the Central Universities of China (DC110305).