Hi Ted,
thanks a lot for this comment.
I implemented it and it works fine.
One thing that for sure is positive of this method is the fact that it
can be performed with the current available 3.2 release, which is
available in maven, whereas for the Optimization solution the
development version is needed.
Anyways I really need to say that I am amazed by the degree of quality
help one gets from this mailinglist.
Thanks a ton,
Best regards,
Andrea
On Sun, Dec 8, 2013 at 12:48 AM, Ted Dunning <ted.dunning@gmail.com> wrote:
> I think that this can actually be solved using linear methods (for the
> quadratic error).
>
> The values of x, x^2, xy, y, y^2 and 1 can be considered coefficients of
> the terms a, b, c, d, e and f. Linear least squares systems of this form
> can be solved directly without use of an optimizer by using QR
> decomposition of the normal equation.
>
> To do this, form the matrix A where each row i is an observation and
> columns are values of x_i, x_i^2, x_i y_i and so on. The rightmost column
> should be the constant 1. We want to solve
>
> A u = z
>
> where u = [a, b, c, d, e, f]' and z = [z_1 ... z_n]' for the solution u
> that results in the smallest error in z. One simple way to do this is to
> note that setting the derivative of this squared error to zero is
> equivalent to solving this equation
>
> A' A u = A' z
>
> exactly. You can go ahead and do this, but if you use the QR
> decomposition, you get more stable results. The QR decomposition breaks A
> into two matrices, Q and R. Q is a orthornormal matrix which means that Q'
> Q = I and R is upper triangular. The normal equation can be rewritten as
>
> (Q R)' (Q R) u = (Q R)' z
> R' R u = R' Q' z
>
> This allows backsubstitution to be used to solve for u.
>
> The particularly nice thing about QR decomposition used in this way is that
> all of the QR implementations that I know of have a solve method.
>
> If you form the matrix A and vector z, then the code that you want in
> commons math is something kind of like this:
>
> DecompositionSolver<http://commons.apache.org/proper/commonsmath/apidocs/org/apache/commons/math3/linear/DecompositionSolver.html>solver
> = new
> RRQRDecomposition<http://commons.apache.org/proper/commonsmath/apidocs/org/apache/commons/math3/linear/RRQRDecomposition.html>
> (A).getSolver<http://commons.apache.org/proper/commonsmath/apidocs/org/apache/commons/math3/linear/RRQRDecomposition.html#getSolver()>
> ();
> RealVector u =
> solver.solve<http://commons.apache.org/proper/commonsmath/apidocs/org/apache/commons/math3/linear/DecompositionSolver.html#solve(org.apache.commons.math3.linear.RealMatrix)>
> (z);
>
>
>
>
> On Sat, Dec 7, 2013 at 7:44 AM, andrea antonello <andrea.antonello@gmail.com
>> wrote:
>
>> >> I didn't expect to have to use weights for this calculation, is there
>> >> a way to define the weights properly (for dummies like me)?
>> >
>> > The weights are mainly used as relatively to each other, if some points
>> > are considered to be more accurate than others. Typically, each point
>> > should have a weight related to the variance of the error for the point,
>> > or to the unit of the measurements (for example when mixing distances
>> > and angles in the residuals, which is not your case).
>> >
>> > If all your points are measured in the same way, they all have the same
>> > variance and hence they should share the same weight. Then, as the
>> > meaning of the weight is only relative, if they all have the same weight
>> > you may well use simply 1.0 as the global weight for everyone.
>>
>> That all makes perfectly sense, thank you very much!
>>
>> Best regards,
>> Andrea
>>
>>
>> >
>> > best regards,
>> > Luc
>> >
>> >>
>> >> Thanks,
>> >> Andrea
>> >>
>> >>
>> >>
>> >>
>> >> On Fri, Dec 6, 2013 at 10:58 AM, Luc Maisonobe <luc@spaceroots.org>
>> wrote:
>> >>> Hi Andrea,
>> >>>
>> >>> Le 06/12/2013 10:02, andrea antonello a écrit :
>> >>>> Hi Ted,
>> >>>> thanks for the reply.
>> >>>>
>> >>>>> How would you like to handle the fact that you may have an infinite
>> number
>> >>>>> of solutions? Will you be happy with any of them? Or do you
>> somehow want
>> >>>>> to find all of them?
>> >>>>
>> >>>> I am afraid my math knowledge goes only to a certain point here
(too
>> >>>> few it seems).
>> >>>> In my usecase I am trying to define a terrain model from sparse
points
>> >>>> (lidar elevation points).
>> >>>> So my idea was to get a set of points (given X, Y, Z) to define
the
>> >>>> parameters of the equation and then be able to calculate an elevation
>> >>>> (Z) for any position I need interpolated.
>> >>>>
>> >>>> I would then expect to have an exact value once the parameters are
>> defined.
>> >>>>
>> >>>> I am not sure how to find the parameters though...
>> >>>
>> >>> For this use case, I would suggest to use any kind of multivariate
>> >>> optimizer. Your Z variable is linear in 6 unknown parameters (a, b,
c,
>> >>> d, e and f). What you want is find values for these 6 parameters such
>> >>> that the quadratic sum of residuals is minimum, i.e. you want to
>> minimize:
>> >>>
>> >>> sum (Z_i  Z(X_i, Y_i))^2
>> >>>
>> >>> Where X_i, Y_i, Z_i are the coordinates for point i and Z(X_i, Y_i)
is
>> >>> the theoretical value from you model equation, which depends on the
6
>> >>> parameters.
>> >>>
>> >>> As Z is linear in the 6 parameters, the sum above is quadratic, so your
>> >>> problem is a least squares problem. Take a look at the
>> >>> fitting.leastsquares package. You have the choice among two solvers
>> >>> there: GaussNewton and LevenbergMarquardt. As you problem is simple
>> >>> (linear Z), both should give almost immediate convergence, pick up one
>> >>> as you want. Start with something roughly along this line (of course,
>> >>> you should change the convergence threshold to something meaningful
for
>> >>> your data):
>> >>>
>> >>> // your points
>> >>> final double[] xArray = ...;
>> >>> final double[] yArray = ...;
>> >>> final double[] zArray = ...;
>> >>>
>> >>> // start values for a, b, c, d, e, f, all set to 0.0
>> >>> double[] initialGuess = new double[6];
>> >>>
>> >>> // model: this computes the theoretical z values, given a current
>> >>> // estimate for the parameters a, b, c ...
>> >>> MultivariateVectorFunction model =
>> >>> new MultivariateVectorFunction() {
>> >>> double[] value(double[] parameters) {
>> >>> double[] computedZ = new double[zArray.length];
>> >>> for (int i = 0; i < computedZ.length; ++i) {
>> >>> double x = xArray[i];
>> >>> double y = yArray[i];
>> >>> computedZ[i] = parameter[0] * x * x +
>> >>> parameter[1] * y * y +
>> >>> parameter[2] * x * y +
>> >>> parameter[3] * x +
>> >>> parameter[4] * y +
>> >>> parameter[5];
>> >>> }
>> >>> return computedZ;
>> >>> }
>> >>> };
>> >>>
>> >>> // jacobian: this computes the Jacobian of the model
>> >>> MultivariateMatrixFunction jacobian =
>> >>> new MultivariateMatrixFunction() {
>> >>> double[][] value(double[] parameters) {
>> >>> double[][] jacobian =
>> >>> new double[zArray.length][parameters.length];
>> >>> for (int i = 0; i < computedZ.length; ++i) {
>> >>> double x = xArray[i];
>> >>> double y = yArray[i];
>> >>> jacobian[i][0] = x * x;
>> >>> jacobian[i][1] = y * y;
>> >>> jacobian[i][2] = x * y;
>> >>> jacobian[i][3] = x;
>> >>> jacobian[i][4] = y;
>> >>> jacobian[i][5] = 1;
>> >>> }
>> >>> return jacobian;
>> >>> }
>> >>> };
>> >>>
>> >>> // convergence checker
>> >>> ConvergenceChecker<PointVectorValuePair> checker =
>> >>> new SimpleVectorValueChecker(1.0e6, 1.0e10);
>> >>>
>> >>> // configure optimizer
>> >>> LevenbergMarquardtOptimizer optimizer =
>> >>> new LevenbergMarquardtOptimizer().
>> >>> withModelAndJacobian(model, jacobian).
>> >>> withTarget(zArray).
>> >>> withStartPoint(initialGuess).
>> >>> withConvergenceChecker(checker);
>> >>>
>> >>> // perform fitting
>> >>> PointVectorValuePair result = optimizer.optimize();
>> >>>
>> >>> // extract parameters
>> >>> double[] parameters = result.getPoint();
>> >>>
>> >>>
>> >>> BEware I just wrote the previous code in the mail, it certainly has
>> >>> syntax errors and missing parts, it should only be seen as a skeleton.
>> >>>
>> >>> best regards,
>> >>> Luc
>> >>>
>> >>>
>> >>>>
>> >>>> Cheers,
>> >>>> Andrea
>> >>>>
>> >>>>
>> >>>>
>> >>>>>
>> >>>>>
>> >>>>>
>> >>>>> On Thu, Dec 5, 2013 at 5:31 AM, andrea antonello <
>> andrea.antonello@gmail.com
>> >>>>>> wrote:
>> >>>>>
>> >>>>>> Hi, I need a hint about how to solve an equation of the
type:
>> >>>>>> Z = aX^2 + bY^2 + cXY + dX + eY + f
>> >>>>>>
>> >>>>>> Is an example that handles such a case available? A testcase
code
>> maybe?
>> >>>>>>
>> >>>>>> Thanks for any hint,
>> >>>>>> Andrea
>> >>>>>>
>> >>>>>> PS: is there a way to search the mailinglist for topics?
I wasn't
>> able
>> >>>>>> to check if this had already been asked.
>> >>>>>>
>> >>>>>>
>> 
>> >>>>>> To unsubscribe, email: userunsubscribe@commons.apache.org
>> >>>>>> For additional commands, email: userhelp@commons.apache.org
>> >>>>>>
>> >>>>>>
>> >>>>
>> >>>> 
>> >>>> To unsubscribe, email: userunsubscribe@commons.apache.org
>> >>>> For additional commands, email: userhelp@commons.apache.org
>> >>>>
>> >>>>
>> >>>
>> >>>
>> >>> 
>> >>> To unsubscribe, email: userunsubscribe@commons.apache.org
>> >>> For additional commands, email: userhelp@commons.apache.org
>> >>>
>> >>
>> >> 
>> >> To unsubscribe, email: userunsubscribe@commons.apache.org
>> >> For additional commands, email: userhelp@commons.apache.org
>> >>
>> >>
>> >>
>> >
>> >
>> > 
>> > To unsubscribe, email: userunsubscribe@commons.apache.org
>> > For additional commands, email: userhelp@commons.apache.org
>> >
>>
>> 
>> To unsubscribe, email: userunsubscribe@commons.apache.org
>> For additional commands, email: userhelp@commons.apache.org
>>
>>

To unsubscribe, email: userunsubscribe@commons.apache.org
For additional commands, email: userhelp@commons.apache.org
