commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <...@spaceroots.org>
Subject Re: [math] solving bivariate quadratic equations
Date Fri, 06 Dec 2013 16:18:06 GMT
Le 06/12/2013 15:26, andrea antonello a écrit :
> Hi Luc,
> thanks a ton for your example code, much appreciated.
> 
> It took me a little while to understand that this works only with the
> development version :)
> 
> I was able to get it working on a small sample of data with results in
> the expected range.
> I didn't check yet really on the result validity because I got a NPE
> because of missing weights.
> In the testcaes I then saw:
> .withWeight(new DiagonalMatrix(weights))
> where the weights are calculated for each point as 1/Z_i.
> 
> 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.

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: Gauss-Newton and Levenberg-Marquardt. 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.0e-6, 1.0e-10);
>>
>>  // 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, e-mail: user-unsubscribe@commons.apache.org
>>>>> For additional commands, e-mail: user-help@commons.apache.org
>>>>>
>>>>>
>>>
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>> For additional commands, e-mail: user-help@commons.apache.org
>>>
>>>
>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>> For additional commands, e-mail: user-help@commons.apache.org
>>
> 
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
> For additional commands, e-mail: user-help@commons.apache.org
> 
> 
> 


---------------------------------------------------------------------
To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
For additional commands, e-mail: user-help@commons.apache.org


Mime
View raw message