Algorithms for Constrained and Weighted Nonlinear Least Squares.
A hybrid algorithm consisting of a Gauss--Newton method and a second-order method for solving constrained and weighted nonlinear least squares problems is developed, analyzed, and tested. One of the advantages of the algorithm is that arbitrarily large weights can be handled and that the weights in the merit function do not get unnecessarily large when the iterates diverge from a saddle point. The local convergence properties for the Gauss--Newton method are thoroughly analyzed and simple ways of estimating and calculating the local convergence rate for the Gauss--Newton method are given. Under the assumption that the constrained and weighted linear least squares subproblems attained in the Gauss--Newton method are not too ill conditioned, global convergence towards a first-order KKT point is proved.
Introduction
. Assume that f : R n continuously differentiable
function and that diagonal matrix with weights
We will discuss the Gauss-Newton method and a second order method for
solving the problem
min
x2R
denotes the 2-norm. For simplicity, and without loss of generality, we
assume that the weights are normalized and sorted such that
The normalization is easily done by first sorting out the zero weights, reducing the
problem and then dividing the remaining nonzero weights with the smallest positive
weight.
To our knowledge all existing algorithms for solving (1.1) are based on the un-weighted
problem
min
Assume that the ordinary Gauss-Newton method
is used to solve (1.3). The search direction, p, is then got by solving
min p2
(1.
rg. Note that (1.4) is solved as an unweighted problem and thus the
condition of this problem is determined by kKk kK y k where K y is the pseudo inverse
of K.
If we on the other hand linearize (1.1), without explicitly multiplying
with the weights, we solve the weighted linear least squares problem
to obtain the search direction p. The condition for the problem (1.5) is mainly determined
by kBk kJk where . For a more detailed discussion on condition
numbers for (1.5) see [11]. The problem (1.4) may be very ill conditioned (regarded
as an unweighted linear least squares problem) despite the fact that (1.5) is well conditioned
(regarded as a weighted linear least squares problem). Obviously it is very
important to look at (1.1) as the class of weighted nonlinear least squares problem.
Another important advantage of using (1.1) instead of (1.3) is that the former
defines a more general problem class than the latter. This is evident if we allow the
weights to be infinitely large. To be more precise, we define the vector -
the equations
weights correspond to zero elements in M . Note that
is the Lagrange multiplier corresponding to the ith constraint and
consequently - i is not defined by (1.6). We will return to the proper way of calculating
these Lagrange multipliers. Problem (1.1) is rewritten, using (1.6), as
Hence, by allowing infinite weights, our original problem formulation (1.1) defines the
class of weighted nonlinear least squares problems with nonlinear equality constraints.
To be even more specific we assume that we have p infinite weights such that
Problem (1.7) can now be stated as
\Theta f T
An equivalent formulation of problem (1.8)
using - is
min
2 . Of course, we could have started by defining our problem as
the one in (1.9) instead of (1.5) (without the need of (1.7) and (1.8)) but then the
notations would get unnecessarily complicated.
In the next section we describe the Gauss-Newton method for solving (1.1). The
local convergence properties of the Gauss-Newton method is analyzed in Section 3
and in Section 4 we show that, under certain assumptions on nondegeneracy, global
convergence is achieved. If the Gauss-Newton method is too slow or does not converge
and second derivatives are available at a reasonable cost then the Newton method
may be used to solve (1.1). However, when there are large and possibly infinite
weights a pure Newton method based on forming the Hessian of g(x) may not work
or, with infinite weights, is not even defined. The natural approach is then to use the
Perturbation method [9] that we will call the generalized Newton-Raphson method
(the gNR method). In Section 5 we construct and analyze an algorithm for solving
(1.1) based on the gNR method. Computational experiments are presented in Section
6 and finally we discuss our results and give hints of possible future work.
2. The Gauss-Newton Method Using the System Equations. In the Gauss-Newton
method, the nonlinear least squares problem (1.1) is linearized around the
current iteration point, x k , and the search direction, p k , is computed as the solution
to
The next iterate is x
the steplength. In the presence of large weights, possibly infinite, it is adequate to
reformulate (2.1) as
\Gammaf#
where we for simplicity have dropped the iteration index k. There are several names
to the linear system of equations in (2.2) such as the equilibrium equations, the system
equations or the augmented system equations. We call (2.2) the system equations and
the matrix in (2.2) is called the system matrix. A less obvious reason for using (2.2)
is that the elements in - corresponding to infinite weights are approximations to the
Lagrange multipliers and that - can be used in a second order method as described
in Section 5.
The following lemma gives the relevant conditions for the system matrix to be
nonsingular.
Lemma 2.1. The system matrix in (2.2) is nonsingular if and only if the rows in
J that correspond to infinite weights are linearly independent and J has full column
rank.
There exist several stable algorithms that solve (2.2), see e.g. [5] for further
references. We have chosen to use the modified QR decomposition, see [5], and the
reasons are the following. The modified QR decomposition is simple and easy to
compute and it is identical to the ordinary QR decomposition when the weights are
equal. The modified QR decomposition is also easily reused in the second order gNR
method, see Section 5.
The modified QR decomposition of J 2 R m\Thetan is defined as
R#
n\Thetan is an upper triangular matrix and \Pi is a permutation
matrix. The decomposition in (2.3), with Q and R nonsingular, exists if and only if
the system matrix in (2.2) is nonsingular (see Lemma 2.1).
The system equations are solved with the modified QR decomposition in the
following way. Using the decomposition (2.3) in (2.2) we get6 6 4
R#
If we make the partition Mm\Gamman ) the solution to (2.4) is
m\Gamman
3. The local rate of convergence for the Gauss-Newton method. In this
section we will describe the local convergence properties of the Gauss-Newton method
described in the previous section. Our analysis depends much upon the perturbation
analysis of the constrained and weighted linear least squares problem done in [11, 12].
After having defined the inverse of the system matrix, using the same notation as in
[11], we state and prove two important theorems on the local convergence rate for
projected residual). In fact, the local convergence properties of
these two quantities are, as we shall see, very similar. Finally we show that J k p k and
the local convergence rate for x
x and J k p k are independent of the parametrization
in R n .
Assuming that b x is a solution of (1.1), we define b
and the corresponding
notation for other quantities evaluated at b
x.
A necessary condition for our algorithm to converge without regularization is that
the system matrix in (2.2) has full rank and it is convenient to make the following
definition.
Definition 3.1. If the system matrix in (2.2) is nonsingular at x we say that x
is a nondegenerate point.
At a nondegenerate point the inverse of the system matrix in (2.2) is given by
is a generalized inverse of J , see [11]. From (2.5) we immediately
get
The following theorem describes the local behaviour of x
Theorem 3.1. Assume that fp k g are generated by solving (2.2) and that all
points x are nondegenerate. If b x is the solution of (1.1) and b
- is the
vector - from (2.2) at b x then
R 1f 00
Proof. From x
f)
Using the Taylor expansion
the first term in (3.4) can be expressed as
To express the second term, \GammaB k
f , in (3.4) we use the perturbation identity (2.2) p.
in [11] which says that
Using the identity
Z 1[f 00
equation (3.6) becomes
The equations (3.5) and (3.7) inserted into (3.4) gives the theorem.
The Gauss-Newton method can be written as
and with b
From Theorem 3.1 we conclude that
and from [8] we get the following theorem.
Theorem 3.2. Define
and - i as the eigenvalues of H x . Then
lim sup
xk
xk
It is easy to get an estimation of the local convergence rate if we use the matrix
defined by (3.2), because then b
R
R \GammaT \Pi T .
A useful quantity for estimating how close x k is to the solution, b x, is the projected
residual is the oblique projection of f k onto R(J k ).
The following theorem shows that J k p k locally has the same convergence behaviour
as
Theorem 3.3. Assume that fp k g are generated by solving (2.2) and that all
points x are nondegenerate. If -(x k ) is the vector - from (2.2) at x k
then
R 1f 00
Proof. Denote the projection J k B k by P k . Then we have
Using the Taylor expansion
multiplying with P k+1
we obtain
I , the equality B k s holds, and we can identify the last term in
equation (3.14) as
From (3.1) we get
and hence
From the perturbation identity (2.1) p. 16 in [11] we get
Using (3.18) and the fact that Y k+1 ffiJ k B k
the equation (3.17) becomes
where the last equality follows from (3.16). The identities -
\Gammap k together with a Taylor expansion of (ffiJ k ) T give
The theorem follows by inserting (3.15) and (3.19) into (3.14).
The matrix corresponding to H x for the projected residual, s k , is
B:
and it is easy to show that H x and H s has the same nonzero eigenvalues. Hence, we
have from Theorem 3.3 the following corollary.
Corollary 3.1. Define B k from the inverse of the system matrix in (3.1). If
lim sup
ks
are the eigenvalues of the matrix H x defined in (3.11).
The relation (3.12) can also be used to determine when ks k+1 k=ks k k reflects
the linear convergence rate and if a second order method should be used. If the
convergence of the Gauss-Newton method is slow we use a higher order method if2
see also Algorithm 6.1.
Several of the above quantities are invariant under a change of parametrization
and as an example we have the following theorem.
Theorem 3.4. The matrix
B:
is independent of the parametrization in R n .
Proof. Assume that
f(x(')). We want to show that
y
B y is the generalized inverse of ry( b
'). Now, consider the Taylor expansion
where
\Delta' T x 00
\Delta' T x 00
By comparing the Taylor expansion (3.21) above with the Taylor expansion
using that J T
that
From (3.23) we finally get
y
which proves the theorem.
A consequence of Theorem 3.4 is that the local convergence for is the same
as for x
The main argument for choosing Jp as a measure of the closeness to the solution
is the following theorem which is a direct consequence of (3.23).
Theorem 3.5. The projection of f on R(J), is independent of the
parametrization in R n .
4. Global convergence. In this section we assume that x k , where k is the
iteration index, is nondegenerate and that p k is the solution of (2.2) at x k . If nothing
else is stated we assume that all limits denoted by ! are when k ! 1, and that all
sums with no explicitly stated upper or lower limit are from one to infinity.
4.1. The merit function. As a merit function we have chosen
The goal is to find a matrix D k of merit weights and a step length ff k , at each
iteration, such that global convergence towards a first order Kuhn-Tucker point can
be proved. To compute D k we will use the approximation
of D). For a fixed matrix D, we define Obviously
a sufficient condition on p k to be a descent direction to \Phi(x; D) at x k is that OE 0
We realize that we can determine a good matrix, \Upsilon(x k ), of merit weights by
solving
min
k\Upsilonk
s:t:
where ffi is a small positive constant and - i is a lower limit for the weights determined
by some previously computed weights, see below. There is always a solution to (4.2)
because
lim
argf min ff
Note that keeping the weights not too large is important in practice but for the global
convergence it is only the constraints in (4.2) that must be satsified. We will now
describe the algorithm for computing the merit weights D k , using \Upsilon(x k ), such that
does not become unnecessarily large. We first describe a method for solving (4.2)
and then an algorithm for computing the actual merit weights D k .
When solving (4.2) we have chosen to use the max-norm since this gives a simple
algorithm. The problem (4.2) can be rewritten as
min kuk1
where u is the diagonal in \Upsilon(x k ), ! is the diagonal in W , and y
with that when Jp is given the problem (4.3) consists of only vectors
and no matrices. The first step in our algorithm is to reduce (4.3) such that u
We then get a new problem
y are the corresponding parts of u; -; ! and y left after the reduction
ae we are ready with the solution -
-. Otherwise
we choose -
y, where e is a vector of ones and thus attain equality in the
constraints. If -
respectively. Again we
can reduce the problem to a copy of (4.4) but where the vectors are shorter and ae
is smaller. The procedure is then repeated until the whole of u is found. It is easily
realized that the infinite weights in ! do not change the algorithm and the algorithm
will terminate with a solution of (4.4).
We determine the actual merit weights D k from the solution \Upsilon(x k ) of (4.2). The
weights may get large close to a saddle point and when the iterates diverge from this
saddle point (that is always the case with the Gauss-Newton method) we would like the
weights to decrease. This is accomplished by saving, say t, older versions,
the merit weight matrices. Initially, at iteration
and at the kth iteration we update
m ), as in Algorithm 4.1.
Algorithm 4.1.
Solve (4.2) for the vector u(x k ).
If d (k)
i be the new
sequence - (1)
In Algorithm 4.2 our Gauss-Newton algorithm is described with line search and
quadratic merit function.
Algorithm 4.2.
Initiate the start vector x k .
while not convergence
Compute
Compute using the modified QR decomposition of J k .
Determine D k from Algorithm 4.1.
Determine the step length ff k such that
4.2. Proving global convergence. We will need the following two technical
lemmas to prove that our algorithm is globally convergent. In the lemmas we use d k
as an arbitrary diagonal element in D k .
Lemma 4.1. Assume that d k - 0; and that fd k g is bounded. Let
be the subsequence of fd k g such that d k j+1 ? d k j . Then the positive series
converges if and only if
converges.
Proof. Take
a
a
where a
converges too. Now assume that b +
N converges to b
N . Hence,
k g is a bounded
sequence that increases to a limit b \Gamma and
converges to b
Lemma 4.2. Assume that an arbitrary component, d k , in the diagonal of D k
stays bounded as k ! 1 and let v k be the corresponding diagonal element in V t .
Then lim and the series
converges.
Proof. Let us first exclude the trivial case that v k becomes equal to the upper
bound ! for a finite k.
The sequence fv k g is an increasing infinite sequence. Hence, lim v k exists and is
denoted v. Take
ffl be an arbitrary small but
fixed positive number. Then d k ? -
k-values. Hence, v ? -
and since ffl ? 0 was arbitrary this implies that v -
d. From d k - v k it follows that
thus we have v -
d - d - v and consequently d k ! v.
Let fd i k g be the subsequence of fd k g with d i k+1 ? d i k . From lemma 4.1 we know
that the series
converges if and only if
converges. Let
us now prove that the latter series converges. From v i k - d i k it follows that
and hence
Since
a subserie of
increases to v, the series
converges. Since
positive series it is sufficient to prove
that it is bounded. Hence, it only remains to prove that the series
converges. Since d the saved older weights are updated in step i 1 . When
we reach d i 1+t
there have been t updates and v i 1+t
equals one of the earlier d
t. In this way we can eliminate both this v i 1+t and the corresponding d i j . In
the same way it is seen that v i 1+t+1
equals one of the d i j . That pair can also be
eliminated from the series. We go on and eliminate elements in this way to get
t. Thus the positive series
bounded and so converges.
That completes the proof.
Our main global convergence theorem covers both bounded and unbounded sequences
of merit weights.
Theorem 4.3. Let fx k g and fD k g be generated by Algorithm 4.2. Assume that
is bounded and that the system matrix in (2.2) is nonsingular in the closure of
g. Then the sequence fx k g has either finite termination at a KKT point or an
accumulation point that is a KKT point of (1.1).
Proof. It is trivial that there is finite termination just at KKT points. Let us now
assume that we have an infinite sequence. Algorithm 4.2 implies that it is sufficient
to consider the following two cases :
These cases will now be treated separately.
There exist a subsequence fx i k g of fx k g such that kD i
is bounded it its possible to choose a subsequence fx j k g of fx i k g such
that x
x for some e
x. From Algorithm 4.2 it follows that kD k k ! 1 only when
is continuous for all points in the closure of fx k g except
KKT points, e
x is both an accumulation point of fx k g and a KKT point.
ii) From the inequality
one can prove that a point e x cannot be an accumulation point of fx k g if there exist
(The proof of (4.5) is a trivial extension of a similar proof in [7] pp. 21-22.)
From Lemma 4.2 we know that
converges and from the Goldstein-
Armijo condition in Algorithm 4.2 then, for a given D k , it follows that for every point
e
x in the closure of fx k g, that is not a KKT point, there exists constants ffl ? 0
that (4.5) is satisfied. Hence, only KKT points remain as possible
accumulation points. That proves the theorem in case ii).
4.3. Line search. We have chosen to keep things simple and therefore we use
a standard cubic interpolation from [3] to approximate the minimum of our merit
function OE(ff). Another, more efficient, line search algorithm can be found in [6].
4.4. Regularization. We use a simple form of subspace minimization described
for the unweighted and constrained case in [7]. We have not been able to prove a
general global convergence result as the one in Theorem 4.3 but as we shall see in the
computational experiments our regularization seems to work appropriately.
5. The generalized Newton-Raphson method. A constrained Newton method
for solving (1.9) can be based on the quadratic subproblem
G)p
are first order approximations
of the Lagrange multipliers. The solution, -
p, to (5.1) is given by the linear
system of equations
G
\Gammaf#
The main disadvantage with using (5.2) is that for very large weights in W 2 the
quadratic subproblem (5.1) and the matrix in (5.2) may be very ill conditioned.
To avoid the ill conditioning due to large weights in W 2 we solve
\Gammaf#
and - is from (2.2). This method is the generalized Newton-Raphson
method [9], or just the gNR method.
The gNR method has an interesting theoretical motivation. Assume that we have
reached a point x k . From the first order approximation (1.5) it is known that x k
solves the perturbed problem
min
is a projection onto R(J k ). Hence, we know
the solution x k of (5.4) and want to compute the solution of the perturbed problem
min
Then we can use the quadratic approximation of z(x) at x k to compute a solution of
problem (5.5) whose error is O(kP k f k k 2 ). If we change back to the original notations
in f(x), this perturbed solution is found by solving problem (5.3) for
From (5.3) it is seen that there exists a matrix N k such that
x as the solution to (1.1) and
Take x . Then from the quadratic approximation in (5.6) we get
From (5.6) it is also seen that J k N k only depends on the surface and not on the
parameterization in x and consequently J k p k is independent of the parameterization
in R n . The generalized Newton-Raphson method is in fact the only quadratically
convergent method with J k p k independent of the parametrization. To see this we
assume that there exists another method which computes e
e
. The series expansion (5.6) is unique and we have J k
e
which implies that e
If we define Z 1 as a matrix whose columns span the null space of J 1 we call p a
descent direction if p T Z T
drawback with both the constrained Newton
method based on (5.2) and the gNR method is that a nonsingular matrix in (5.2)
or (5.3) is not sufficient for p to be a descent direction. However, we use the gNR
method only when we are close to the solution, see (3.20) and Algorithm 6.1, and
therefore we use the gNR method undamped. From (2.2) we get -, needed for G, and
the Gauss-Newton search direction and if the matrix in (5.3) is singular we use the
already available Gauss-Newton direction.
If we use the modified QR decomposition to solve (2.2) it is possible to reduce the
size of the system in (5.3). Ignoring the permutation matrix it is possible to rewrite
R#
Now implies that Q and we can reduce (5.7) to
R T \GammaG
\Gamma-
are the first n elements in Q T - and Q
respectively.
The matrix in (5.8) may be indefinite and we must either use a stable method
for indefinite systems, see e.g. [4], or add some condition on the submatrices in (5.8).
One possibility of the latter kind is to assume that R is well conditioned and use R T
to reduce (5.8) to
R T \GammaG
"\Gamma-
The solution is the matrix R+MnR \GammaT G is nonsingular,
if not we take a Gauss-Newton step.
6. Computational experiments. The algorithm we use in our tests is shown
below.
Algorithm 6.1.
Initialize
while
Determine the Jacobian J and the vector f .
Compute the GN direction p and - by solving (2.2).
If regularization was needed then Second := false.
If Close and Second and Rate ? 0:5
Compute the gNR direction, p gNR , by solving (5.3).
If the matrix in (5.3) is nonsingular then
If GN
Compute the merit weights by Algorithm 4.1.
Determine the step length ff using the line search described
in Section (4.3) with the merit function OE(ff).
x
To use a pure GN method then the variable Second has a fixed value of false.
We have tested our algorithm on three different problems described in the Ap-
pendix; Schittkowski 308 [10], Boggs 2 and Boggs 8 [2]. The intention with the tests is
not to show that the algorithms are faster than other existing algorithms but to show
how our algorithms handles large weights and inadequate models (ill conditioning in
the linear problems). Another important aim with the tests is to verify our theoretical
results on the local convergence rate. Therefore it has been natural to use small and
simple test problems.
We define as two different
measures of the convergence rate for the Gauss-Newton method. We emphasize that
k is an excellent way of estimating the convergence rate when regularization is not
needed and when b x is not known.
The first problem, Schittkowski 308, is first solved with the Gauss-Newton method
and the result is in Tab. 1. The largest weight is 10 20 and if the weights are multiplied
explicitly with f , forming then the algorithm breaks down because of
numerical instability. Note the slow growth of the merit weights. The first problem
Schittkowski 308 with the Gauss-Newton method.
5 7.3e-3 6.6e-5 2.7e-3 4.6e-2 4.4e-2 3.0 99 1.0
6 3.4e-4 3.0e-6 1.2e-4 4.6e-2 4.6e-2 2.9 3.9 1.0
8 7.2e-7 6.5e-9 2.7e-7 4.6e-2 4.7e-2 3.0 1.0e+2 1.0
9 3.3e-8 3.0e-10 1.2e-8 4.6e-2 4.6e-2 3.0 1.0e+2 1.0

Table
Schittkowski 308 with the gNR method.
5   5.9e-5 5.6e-12 2.2e-5 1.0
6   2.9e-11 2.3e-16 1.1e-11 1.0
solved with the gNR method is showed in Tab. 2. The asterisk indicates that the
gNR method was used in that step. The second problem, Boggs 2, is a constrained
problem and it has been solved with the Gauss-Newton method, Tab. 3, and the gNR
method, Tab. 4. All the merit weights for the Gauss-Newton method were equal to
one and are not shown in the Tab. 3. The remaining two test problems illustrate the
regularization. The rank of the problem is shown under the headline Rank. In Tab.
5 the second test problem, Boggs 2, is solved with the Gauss-Newton method when
the Jacobian is rank deficient at the starting point. In the third problem, Boggs 8,
the Jacobian at the solution is rank deficient and the result is shown in Tab. 6.
7. Discussion. We claim that we have developed an efficient and fairly robust
algorithm for solving (1.1) (with possibly infinite weights as discussed in the intro-
duction). However, it is difficult for us to measure the effectiveness of the algorithm
Boggs 2 with the Gauss-Newton method.
3 1.2e-2 2.8e-2 2.3 0.13 0.50 2.6 1.0
5 2.4e-4 2.5e-4 7.3e-3 0.10 1.9e-2 1.6e-2 1.0
6 6.4e-5 6.6e-5 3.0e-4 0.27 4.2e-2 1.6e-2 1.0

Table
Boggs 2 with the gNR method.
3 1.2e-2 2.8e-2 2.3 1.0
5 2.4e-4 2.5e-4 7.3e-3 1.0
6 6.4e-5 6.6e-5 3.0e-4 1.0
9   8.2e-16 1.0e-15 1.9e-15 1.0
because there are, to our knowledge, no other algorithm that can solve such a general
problem as (1.1).
The local convergence properties are well understood for the Gauss-Newton algo-
rithm. It is especially interesting that the local convergence results are valid for the
whole problem class defined by (1.1) and that they are independent of the parametrization
in R n .
The merit function is especially suited for our weighted and constrained problem
and our technique for choosing the merit weights is effective and do not lead to
unnecessary large weights.
Boggs 2, Gauss-Newton and rank deficient at the starting point.
9 0.47 0.71 0.14 1.0 0.12 3
As for robustness, we have shown that our algorithm is globally convergent when
the iteration points are nondegenerate. It remains to find a way to regularize when
the rows in J corresponding to very large weights become (almost) linearly dependent.
We believe that this is a difficult and challenging problem to solve.


Appendix

. Test problems. In this appendix we define our three test problems
and the weight sequences. We also give the starting points, x start ; solutions, b
the residuals f(bx). The examples are from [10] and [2] and includes unconstrained as
well as constrained problems.
Schittkowski 308 [10]
An unconstrained problem which we have modified by incorporation of weights.
A constrained problem where the Jacobian is rank deficient at the second
Boggs 8, Gauss-Newton and rank deficient at the solution.
28 1.1 1.0 0.65 1.0 3.4e-10 5
29 1.1 0.99 1.6 3.5 9.5e-7 4
53 4.2e-9 0.25 0.50 3.5 1.0 4
54 1.0e-9 0.25 0.50 3.5 1.0 3
starting point, x start2 .
Boggs 8 [2]
A constrained problem where the Jacobian is rank deficient at the solution.



--R


A Strategy for Global Convergence in a Sequential Quadratic Programming Algorithm
Numerical methods for unconstrained optimization and nonlinear equations
Matrix Computations
Modifying the QR decomposition to weighted and constrained linear least squares


Iterative solution of nonlinear equations in several variables
A comparision of some algorithms for the nonlinear least squares problem
More Test Examples for Nonlinear Programming Codes
Perturbation theory and condition numbers for generalized and constrained linear least squares problems

--TR

--CTR
Hiroshi Hosobe, Hierarchical nonlinear constraint satisfaction, Proceedings of the 2004 ACM symposium on Applied computing, March 14-17, 2004, Nicosia, Cyprus
nonlinear least squares;parameter estimation;optimization;weights