Re: Re: about ConstrainedMin
- To: mathgroup at smc.vnet.net
- Subject: [mg29888] Re: [mg29855] Re: [mg29806] about ConstrainedMin
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Tue, 17 Jul 2001 01:00:33 -0400 (EDT)
- References: <200107140536.BAA18344@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
qing.cheng at icos.be wrote:
>
> Thank you, Mark.
>
> But I would like to know what ConstrainedMin inside does, not only the
> usage.
> The problem was rised from one of our applications. There we need to
> measure the position of each pins in a leads component(electronic chip).
> Based on these individual points, we need to calculate a plane, which
> should reflect a physical plane where the component can 'sit' stably, so
> call it seating plane. One way we try to achieve this is to convert this
> problem to a LP problem. The objective function is to Minimize the Sum
> distance between the measurement points and the plane.
> We have implemented a Simplex method besed on the algorithm in Numerical
> Recipes to solve this linear problem, and found it worked not very well for
> ">=" type constraints. I also brought the same problem to Mathematica, and
> found the situation that BasicSimplex failed as same as our C
> implementation, while ConstrainedMin found good solution. Now, we have done
> a data transformation before pass them to Simplex algorithm to ensure that
> all the constraints are "<=". It works in that way. But still I would like
> to know how ConstrainedMin improved BasicSimplex. (In Mathematica hand book
> from Stephen Wolfram, page 1061 says that ConstrainedMax and related
> function use an enhanced version of the simplex algorithm).
>
> Could you give me some more information or suggestions about it?
>
> Best Regards.
>
> /Qing
> [...]
I do not know the details of the Mathematica implementation of the
simplex algorithm. I will make a few comments about other aspects of
this particular problem, in case any are of relevance e.g. for obtaining
preferable or faster results.
As worded, the problem above appears to be one of total least squares.
That is, the task is to minimize distances from a plane (ordinary least
squares would minimize squares of vertical discrepancies; we'll address
that later).
First note that a plane may be represented by an equation of the form
a*x+b*y+c*z==d. Provided the x-coefficient is not zero we can normalize
so that it is one. We'll work with an example. We take the plane below,
find some points on it, and perturb them. We then hope to recover a
plane with approximately the same parameters.
plane = 3*x + 2*y + z - 10;
npoints = 10;
SeedRandom[1111];
xyvals = 2*Table[Random[], {npoints}, {2}];
zvals = Map[10 - 3*#[[1]] - 2*#[[2]] &, xyvals];
points = Map[Flatten, Transpose[{xyvals, zvals}]];
offsets = Table[.4*Random[] - .2, {npoints}, {3}];
fuzzedpoints = points + offsets;
Map[10 - 3*#[[1]] - 2*#[[2]] - #[[3]] &, fuzzedpoints];
Given the plane a*x+b*y+c*z==d and the point {x0,y0,z0} the square of
the distance from point to plane is given by
(a*x+b*y+c*z-d)^2/(a^2+b^2+c^2). If we want to minimize the sum of these
square distances we may discard the common denominators. We form the
resulting sum as below.
newpoints = Map[Append[#, -1] &, fuzzedpoints];
vec = {1, b, c, d};
dotprods = newpoints.vec;
sumsquares = dotprods.dotprods;
Now we can use FindMinimum to get the best values for the parameters (at
least assuming we can give it reasonable starting values).
In[16]:= bestparams = FindMinimum[sumsquares, {b,1}, {c,1}, {d,1}]
Out[16]= {0.117926, {b->0.654815, c->0.357271, d->3.50226}}
We'll check the result.
In[17]:= bestplane = x + b*y + c*z - d /. bestparams[[2]]
Out[17]= -3.50226 + x + 0.654815 y + 0.357271 z
In[18]:= Expand[3*bestplane]
Out[18]= -10.5068 + 3 x + 1.96445 y + 1.07181 z
We could instead minimize the sum of distances rather than square
distances.
sumdistances = Apply[Plus, Map[Abs, dotprods]];
For this we do not have a differentiable function so I used a secant
method by giving two initial values.
In[26]:= bestparamsl1 = FindMinimum[sumdistances, {b,.5,1}, {c,.5,1},
{d,1,2}]
Out[26]= {0.879003, {b -> 0.636147, c -> 0.34802, d -> 3.42449}}
In[27]:= bestplanel1 = x + b*y + c*z - d /. bestparamsl1[[2]]
Out[27]= -3.42449 + x + 0.636147 y + 0.34802 z
In[28]:= Expand[3*bestplanel1]
Out[28]= -10.2735 + 3 x + 1.90844 y + 1.04406 z
I found that it took some experimentation to make this work; "bad"
initial values gave bad results. Likewise for avoidung Abs and using
instead:
sumdistances = Apply[Plus, Map[Sqrt[#^2]&, dotprods]]
This latter does not have differnetiability problems, hence again does
not need two initial values. But I did not get good results unless I
started reasonably close to the actual solution.
Another scenario is that you might want to treat x and y values as
"exact" and z values as subject to experimental error, and then
determine a best-fitting plane of the form z = a*x+b*y+d. This is
readily cast as an ordinary linear least-squares optimization problem.
First we separate z values from the rest in the data.
vec2 = {a, b, d};
pointsnoz = newpoints[[All, {1, 2, 4}]];
zvals = newpoints[[All, 3]];
pointsnoz can be regarded as a matrix by which we multiply vec2; were
the data to lie exactly on a plane we would have zvals == pointsnoz.vec2
for appropriate values of the parameters in vec2. In the exactly
determined case one can multiply by an inverse matrix to solve for
{a,b,d}. For a least-squares solution, one instead may multiply by the
(Moore-Penrose) pseudo-inverse.
In[43]:= bestparams2 = Thread[vec2 -> PseudoInverse[pointsnoz] . zvals]
Out[43]= {a -> -2.68791, b -> -1.81179, d -> -9.63997}
In[44]:= bestplane2 = z - (a*x + b*y + d) /. bestparams2
Out[44]= 9.63997 + 2.68791 x + 1.81179 y + z
Note that there are ways to do this that avoid explicit computation of
PseudoInverse. In version 4.1 of Mathematica the "further notes" for
QRDecomposition and SingularValues demonstrate some such methods.
Moreover one can formulate the problem as a single call to Fit.
One reason to prefer this second approach is that it is alot faster. If
you have 1000 data points the first method will take several seconds and
the second a split second. This is probably not applicable for the
stated problem (can one work with 1000 pins in a single chip? What would
the connection be like?)
I was initially at a loss to see how this might be formulated as a
linear programming problem. Then one possibility came to mind, that you
constrain the plane so that it lies "under" the data. In this way we get
distances without recourse to Abs and we have a legitimate LP problem.
constraints = Map[# >= 0 &, newpoints.vec];
In[46]:= objfunc = Apply[Plus, newpoints.vec]
Out[46]= 12.9605 + 8.96492 b + 45.3205 c - 10 d
In[47]:= cbestparams = ConstrainedMin[objfunc, constraints, {b, c, d}]
Out[47]= {1.10822, {b -> 0.602984, c -> 0.353536, d -> 3.32804}}
In[48]:= cbestplane = x + b*y + c*z - d /. cbestparams[[2]]
Out[48]= -3.32804 + x + 0.602984 y + 0.353536 z
In[49]:= Expand[3*cbestplane]
Out[49]= -9.98412 + 3 x + 1.80895 y + 1.06061 z
I believe there is a way to formulate the total least squares problem so
that one can use linear algebra techniques e.g. SingularValues. But I do
not know how to do that myself. One advantage that may make such an
approach attractive (other than speed, if you have a truly monstrous
chip) is that these methods are immune to the problems FindMinimum can
have with respect to starting values.
Daniel Lichtblau
Wolfram Research
- Follow-Ups:
- non-commuting objects with dummy indices
- From: K G Zloshchastiev <zlosh@yahoo.com>
- non-commuting objects with dummy indices
- References:
- Re: about ConstrainedMin
- From: qing.cheng@icos.be
- Re: about ConstrainedMin