Abstract On the basis of the ideas, proposed by one of the authors (I.N. Silin), a suitable software has been developed for constrained data &quot;tting. Constraints may be of arbitrary type; i.e. equalities and inequalities. The simplest possible way has been used. The widely known program FUMILI was re-written in the C# # language. Constraints in the form of inequalities ( )5a were taken into account by changing them into equalities ( )&quot;t and simple inequalities of type t5a. The equalities were taken into account by means of quadratic penalty functions. The suitable software was tested on the model data for the ANKE setup (COSY accelerator, Forschungszentrum JuK lich, Germany)

I N Silin

S N Dymov

S V Yaschenko

V S Kurbatov

CiteSeerX

Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437
Constrained minimization in the C## environment
S.N. Dymov, V.S. Kurbatov, I.N. Silin, S.V. Yaschenko*
Laboratory of Nuclear Problems, Joint Institute for Nuclear Research, 141980 Dubna, Russia
Received 11 December 1998; received in revised form 24 June 1999 accepted 6 July 1999
Abstract
On the basis of the ideas, proposed by one of the authors (I.N. Silin), a suitable software has been developed for
constrained data "tting. Constraints may be of arbitrary type; i.e. equalities and inequalities. The simplest possible way
has been used. The widely known program FUMILI was re-written in the C## language. Constraints in the form of
inequalities /(h)5a were taken into account by changing them into equalities /(h)"t and simple inequalities of type
t5a. The equalities were taken into account by means of quadratic penalty functions. The suitable software was tested
on the model data for the ANKE setup (COSY accelerator, Forschungszentrum JuK lich, Germany). ( 2000 Elsevier
Science B.V. All rights reserved.
Keywords: FUMILI; ANKE; Codes
1. Introduction
In the present paper we describe two realizations
(in C## language) of constrained minimization
for s2-like functionals. One of them is the algorithm
of the FUMILI code, which was available for users
as a part of CERN library [1]. The description of
this algorithm was published in Russian [2] at the
end of the 1960s. Due to the fact that the access to
this publication is not easy for an English reader,
we give a short description of the FUMILI algo-
rithm. This algorithm is now coded in the C##
language.
The second part is the realization of the idea
proposed by one of the authors (I.N. Silin) for
solving the constrained minimization problem in
a general case, where constraints are of arbitrary
*Corresponding author.
type (arbitrary equalities and inequalities) [3].
Technically, here constraints are taken into ac-
count by the method of penalty functions (though
there are other ways of doing this [3]). The algo-
rithm described below was tested on the model
data for the calibration process ppPdp` under
the conditions of the ANKE setup [4].
2. Algorithm of FUMILI
For simplicity, let us assume that the function to
be minimized has the form1
s2"
1
2
n
+
j/1
A
f
j
(x
j
, h)!F
j
p
j
B
2
(1)
1What follows can be easily generalized to the case where the
covariance matrix of the data F
j
has non-diagonal terms.
0168-9002/00/$ - see front matter ( 2000 Elsevier Science B.V. All rights reserved.
PII: S 0 1 6 8 - 9 0 0 2 ( 9 9 ) 0 0 7 5 8 - 5
where f
j
(x
j
, h) are the measured functions at the
points x
j
, F
j
are the measured values, p
j
are their
errors, and h are parameters to be estimated.
The minimum condition is
Ls2
Lh
i
"
n
+
j/1
1
p2
j
)
Lf
j
Lh
i
[ f
j
(x
j
, h)!F
j
]"0, i"1,2, m
(2)
where m is the number of parameters.
Expanding the left hand side of Eq. (2) in para-
meter increments and retaining only linear terms
we get
A
Ls2
Lh
i
Bh/h0#+k A
L2s2
Lh
i
Lh
k
Bh/h0 ) (hk!h0k )"0.
Here h0 is some initial value of parameters. In
a general case
L2s2
Lh
i
Lh
k
"
n
+
j/1
1
p2
j
)
Lf
j
Lh
i
)
Lf
j
Lh
k
#
n
+
j/1
( f
j
!F
j
)
p2
j
)
L2f
j
Lh
i
Lh
k
.
(3)
In the FUMILI algorithm an approximate ex-
pression (3) for L2s2/Lh
i
Lh
k
is used in which the last
term is discarded (it is often done, not always wit-
tingly, and sometimes causes trouble), i.e.:
L2s2
Lh
i
Lh
k
+Z
ik
"
n
+
j/1
1
p2
j
)
Lf
j
Lh
i
)
Lf
j
Lh
k
.
As a result the equations for parameter increments
have the following form:
A
Ls2
Lh
i
Bh/h0#+k Zik ) (hk!h
0
k
)"0, i"1,2, m.
A remarkable feature of the algorithm is the
technique used for step restriction. For the current
approximation of the optimal parameters h0 a par-
allelepiped P
0
is built with a centre at h0 and axes
parallel to the coordinate axes h
i
. The lengths of the
parallelepiped sides along the i-th axis are 2 ) b
i
,
where b
i
have such values that the functions f
j
(h)
are quasi-linear all over the parallelepiped. If the
step *h gives a new point h1"h0#*h outside P
0
,
the crossing h1 of the vector *h with the surface of
P
0
is found and taken as a new value for the
parameters. After selection of a new value for
the parameters, it is checked whether the function
reduction is big enough compared with that ex-
pected in the quadratic approximation. If it is not,
the step reduction is performed. Some parallel-
epiped lengths can be increased too.
In addition, FUMILI takes into account simple
linear inequalities such as
h.*/
i
4h
i
4h.!9
i
. (4)
They form a parallelepiped P (P
0
may be deformed
by P). If the value of the parameter coincides with
its restriction (the approximation lies on the surface
of P) and the gradient component is such that s2 is
not going to increase beyond P, the corresponding
parameter is "xed.
Then the step is calculated for all non-"xed para-
meters and if some parameters, lying on the surface
of P, go beyond P, one of them is temporarily "xed
too (the parameter, for which the ratio
D*h
i
D/J(Z~1)
ii
is maximal), and so on.
The criterion to be ful"lled for the end of the
iteration process is that all parameters are "xed due
to only the gradient component signs and step
increments for non-"xed parameters
D*h
i
D(e )J(Z~1)
ii
where e is a small number &0.01. As the number of
"xation combinations is "nite, the number of steps
will also be "nite, at least in the convex quadratic
case.
Very similar step formulae are used in FUMILI
for the negative logarithm of the likelihood func-
tion with the same idea of linearizing the functional
argument.
3. Minimization of s2 functionals with arbitrary
constraints
3.1. Formulation of the problem
Again, let us assume that the function to be
minimized has the same form (1), but in addition to
simple linear constraints (4) there are two more
types of constraints, i.e. non-linear inequalities and
equalities:
a
r
4/
r
(h)4b
r
, r"1,2,md
(5)
t
s
(h)"c
s
, s"1,2, me
. (6)
432 S.N. Dymov et al. / Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437
Here /
r
(h), t
s
(h) are the regular functions of the
parameter h; a
r
, b
r
, m
d
are the low and upper
boundaries of the inequalities and their number;
c
s
, m
e
are any constant and number of equations.
The regularity is taken to mean continuous sec-
ond-order derivatives. The problem of taking into
account the constraints in the form of the equalities
of type (6) was solved before [5}7]. As for the
constraints in the form of inequalities (5), the
authors did not know a simple solution until one of
them (I.N. Silin) proposed a method for taking
them into account [3]. According to Ref. [3], any
constraint of the form a
r
4/
r
(h)4b
r
can be re-
placed by a simple inequality and equality, namely
a
r
4t
r
4b
r
(7)
/
r
(h)"t
r
. (8)
Here t
r
is an additional variable constrained by two
boundaries a
r
, b
r
, (8) is a constraint in the form of
the equation. You can see that constraints (7) have
the same form and structure as those of Eq. (4), so
we can combine them and introduce just one type
of simple constraint:
h.*/
i
4h
i
4h.!9
i
where index i changes in a wider range:
i"1,2, m,2, m#md
and for i'm
h
i
"t
i~m
, h.*/
i
"a
i~m
, h.!9
i
"b
i~m
.
Then the problem of constrained minimization
in a general case can be reformulated as follows. We
"nd a minimum of function (1) under the con-
straints:
h.*/
i
4h
i
4h.!9
i
, i"1,2, m,2, m#md
(9)
m
u
(h)"d
u
, u"1,2, md
,2, md
#m
e
(10)
where for 14u4m
d
, m
u
"/
u
, d
u
"t
u
and for
m
d
(u4m
d
#m
e
, m
u
"t
u~md
, d
u
"c
u~md
.
After such reformulation the number of the para-
meters to be "tted and of the constraints in the form
of simple inequalities of type (9) becomes m#m
d
;
and the number of constraints in the form of equa-
tions of type (10) becomes m
d
#m
e
.
When non-simple constraints are only equations
they can be taken into account by using either the
method proposed in Ref. [7] or the penalty func-
tion method. Here we use the latter. In the penalty
function method a minimum of function (11) is
searched for as ¹PR.
U"
1
2
n
+
j/1
A
f
j
(x
j
, h)!F
j
p
i
B
2
#
1
2
¹ A
md
+
r/1
(/
r
!h
r`m
)2
p2
r
#
me
+
s/1
(t
s
!c
s
)2
p2
s
B.
(11)
Here ¹ is the penalty factor (normally it is a su$-
ciently big number), and p
r
, p
s
are the formally
calculated errors of constraints.
3.2. Iteration scheme
Let us rewrite Eq. (11) in the form
U"U
1
#
1
2
md
+
r/1
(/
r
!h
r`m
)2
w
r
where
U
1
"
1
2
n
+
j/1
A
f
j
(x
j
, h)!F
j
p
j
B
2
#
1
2
¹
me
+
s/1
(t
s
!c
s
)2
p2
s
and w
r
"p2
r
/¹. With a chosen ¹ the minimum
condition is
LU
Lh
k
"
LU
1
Lh
k
#
md
+
r/1
L/
r
Lh
k
)
(/
r
!h
r`m
)
w
r
"0,
(k"1,2, m) (12)
LU
Lh
r`m
"!
(/
r
!h
r`m
)
w
r
"0, (r"1,2, md
). (13)
In both Eqs. (12) and (13) derivatives are taken
only for those parameters which are not "xed, i.e.
kOi
&
, r#mOi
&
, where i
&
is the index of a "xed
parameters. The functions on the left-hand sides of
Eqs. (12) and (13) depend on the m#m
d
para-
meters. Near the minimum we can expand the
left-hand sides of the equations in parameter in-
crements retaining only linear terms. For Eq. (12)
we have
C
LU
1
Lh
k
#
md
+
r/1
L/
r
Lh
k
)
(/
r
!h
r`m
)
w
r
D
#
m
+
l/1
C
L2U
1
Lh
k
Lh
1
#
md
+
r/1
L/
r
Lh
k
)
L/
r
Lh
1
)
1
w
r
D ) dh1
!
md
+
r/1
L/
r
Lh
k
)
dh
r`m
w
r
"0. (14)
S.N. Dymov et al. / Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437 433
We wrote Eq. (14) in the approximation of the
functional argument linearization method [8], in
which the derivatives L2//Lh
k
Lh
l
are discarded. All
values of functions and derivatives are taken for the
current values of the parameters. Let us also note
that index l (lOi
&
) in the second term runs over
indices of non-"xed parameters. By analogy, for Eq.
(13).
[/
r
!h
r`m
]#
m
+
l/1
L/
r
Lh
l
) dh
l
!dh
r`m
"0. (15)
From Eq. (15), for the non-"xed parameter
h
r`m
(r"1,2,md
) we have
dh
r`m
"[/
r
!h
r`m
]#
m
+
l/1
L/
r
Lh
l
dh
l
. (16)
Substituting Eq. (16) into Eq. (14) and after some
algebra we obtain
G
k
#
m
+
l/1
Z
kl
) dh
l
"0 (17)
where
G
k
"
LU
1
Lh
k
#
md
+
r/1
C
L/
r
Lh
k
)
(/
r
!h
r`m
)
w
r
D
Z
kl
"
L2U
1
Lh
k
Lh
l
#
md
+
r/1
L/
r
Lh
k
)
L/
r
Lh
l
)
1
w
r
.
A remarkable feature of the last expressions is
that the index l runs only over non-"xed para-
meters l"1,2,m, and the index r runs only over
those inequalities for which additional parameters
t
r
(7) are "xed!
Finally, the solution of Eq. (17) is
dh"!(Z~1 )G).
The increments of the additional parameters
dt
r
"dh
r`m
are calculated according to formula
(16).
The advantage of this iteration scheme is that the
matrix inversion only of order m]m is done irre-
spective of the number of constraints.
The scheme can be improved in order to avoid
the in#uence of non-linear e!ects of permanently
valid inequalities. To do this h
r`m
which is neither
on low nor upper boundary is substituted by the
quantity t
r
(h) (or its boundary value if the former
lies outside the boundary) after every iteration step.
4. Test
Both realizations described above are coded in
C## and tested on the model data for the calib-
ration reaction ppPdp` under the conditions of
the ANKE setup [4]. According to the plans,
ANKE will consist of three detectors: a side de-
tector, forward and backward ones. At the moment
the side detector is fully assembled, only a scintilla-
tion hodoscope is ready for the forward detector.
The side detector consists of two scintillation hodo-
scopes (START, STOP), and two proportional
chambers each consisting of three sensitive planes.
It allows one to reconstruct all the kinematic para-
meters of the ejectiles passing through the side
detector. The scintillation hodoscope incorporated
in the forward detector is capable of measuring the
coordinates of the particle and its time of #ight.
The "rst data were obtained in May and July
1998, the accuracies being studied. As the main
calibration process required for the analysis of de-
tector performance is the reaction ppPdp`, we
took this process for the tests. A number of events
was simulated for the beam kinetic energy
¹
"%!.
"425 MeV with the p` meson passing
through the side detector and the deuterons passing
through the scintillation hodoscope of the forward
detector. Simulation was done by the GEANT code
with all physical processes switched on except the
decay of p` mesons. In the case where the kin-
ematic parameters of the beam proton and second-
ary p` are known, there is one constraint having
the form of an equality, namely, the missing mass of
the process should be equal to the mass of the
deuteron:
(E
"%!.
#M
1
!E
p
` )2!( p
"%!.
!p
p
`)2"M2
$
(18)
where E
"%!.
, E
p`
are the energies of the beam pro-
ton and the secondary p` meson, p
"%!.
, p
p
` are
their 3-momenta, M
1
, M
$
are the masses of the
proton and the deuteron, respectively.
As noted above, for the deuterons detected by the
forward hodoscope, their coordinates and times of
#ight (t
$
) will be measured hopefully with the
434 S.N. Dymov et al. / Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437
Fig. 1. Accuracies of determining the particle kinematic parameter.
accuracies permitting 4c "t (using all 4 conserva-
tion laws). As at the moment not all accuracies are
known, we assume that their coordinates and times
of #ight lie between some boundaries and put for-
ward requirements having the form of the following
three inequalities:
y
.*/
4y
d
4y
.!9
, z
.*/
4z
d
4z
.!9
,
t
.*/
4t
$
4t
.!9
. (19)
The "rst two requirements come from the geo-
metrical dimensions of the scintillation hodoscope,
whereas the last one does from the simulation data.
Three functions y
d
, z
d
, t
d
were expressed as the ones
(in the form of polynomials to the third-order in-
clusive) of two angles of the pion in the laboratory
system of coordinates.
The total number of "tted parameters was six,
the "rst three are angles h
xz
, h
yz
of the pion relative
to the beam proton and the pion momentum in the
laboratory system. The last three parameters were
additional parameters t
r
, corresponding to three
inequalities (19). Initial pion angles were always 0,
and initial momenta were calculated as a function
of these angles. The coordinates of the pion detec-
ted in the side detector were expressed as functions
of three pion variables, i.e., two angles and mo-
mentum. The total number of events was &3000,
the maximum number of iterations was 40.
Two "ts corresponding to two di!erent realiz-
ations, described above, were performed. In the "rst
"t the constraint of the form of non-linear equation
(18) was disabled, while in the second one it was
enabled. In Fig. 1 the accuracies for both realiz-
ations are shown. Figs. 1(a)}(c) are for the "rst "t,
Figs. 1(d)}(f) are for the second one. It is necessary
to stress a drastic improvement of accuracy in *p/p
in the second case, which is the result of additional
constraint.
S.N. Dymov et al. / Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437 435
Fig. 2. Illustration of the Richardson approximation.
Each event was "tted to three values of the pen-
alty factor ¹. The initial value of this factor was
selected by using the formula
¹"100 )
n
%91
n
#0/
where n
%91
is the number of experimental points,
and n
#0/
is the number of constraints. In our case
n
%91
"6, n
#0/
"4.
Each successive value of ¹ was ten times larger
than the previous one. According to Ref. [3], in this
case, according to Richardson, we should have the
convergence, i.e. the parameter improvements
*
32
"p
3
!p
2
should be 10 times smaller than
*
21
"p
2
!p
1
. Here p
1
, p
2
and p
3
mean the values
of the "tted parameter for ¹ equal to ¹
1
, ¹
2
and
¹
3
, respectively.
In Fig. 2 the ratios *h32
xz
/*h21
xz
, *h32
yz
/*h21
yz
,
*p32/*p21, *(s2)32/*(s2)21 are shown. It is seen
that they are close to 10, and this indicates that the
statements made in Ref. [3] are correct.
5. Conclusion
Two codes are developed for the minimization of
s2-like functionals in the C## language. One of
them is realization of the FUMILI code with con-
straints of the form of simple boundaries. The sec-
ond one is the minimization with constraints of any
type. With FUMILI as a starting point, the C##
code is developed and tested on model data. The
results of the test show a high performance of the
algorithms developed. In conclusion, the authors
express their gratitude to their colleagues from the
436 S.N. Dymov et al. / Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437
ANKE collaboration for providing the necessary
details and to L. Pashkevich for her help in prepar-
ing an English version of the paper.
References
[1] I.N. Silin, CERN Program Library D 510, FUMILI, 1983.
[2] I.N. Silin, Appendix III, Statisticheski metodi v eksperimen-
talnoi "zike, Atomizdat 1976 (translated into Russian from
W.T. Eadie et al., Statistical Methods in Experimental Phys-
ics, CERN, Geneva, 1971, North-Holland, Amsterdam).
[3] I.N. Silin, Resolute progress in a constrained minimization
problem, JINR Rapid Communications, No 3 [89]-98, pp.
25}30.
[4] Forschungszentrum JuK lich, Germany, Annual Report 1996
(all the details about ANKE setup and related bibliography
can be found in this report).
[5] J.P. Berge, F.T. Solmitz, H.D. Taft, Rev. Sci. Instr. 32 (1961)
538.
[6] R. Bock, CERN 60-30 (1960).
[7] V.S. Kurbatov, I.N. Silin, Nucl. Instr. and Meth. A 345
(1994) 346.
[8] S.N. Sokolov, I.N. Silin, Preprint JINR D-810, Dubna,
1961.
S.N. Dymov et al. / Nuclear Instruments and Methods in Physics Research A 440 (2000) 431}437 437


Constrained minimization in the C# # environment

http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=12CA83CFADB81C4705C2212F413365E3?doi=10.1.1.1074.6101&rep=rep1&type=pdf

Constrained minimization in the C# # environment

Abstract

Similar works

Full text

Available Versions

CiteSeerX