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)? Thanks, Andrea On Fri, Dec 6, 2013 at 10:58 AM, Luc Maisonobe <[email protected]> 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 <[email protected] >>>> 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: [email protected] >>>> For additional commands, e-mail: [email protected] >>>> >>>> >> >> --------------------------------------------------------------------- >> To unsubscribe, e-mail: [email protected] >> For additional commands, e-mail: [email protected] >> >> > > > --------------------------------------------------------------------- > To unsubscribe, e-mail: [email protected] > For additional commands, e-mail: [email protected] > --------------------------------------------------------------------- To unsubscribe, e-mail: [email protected] For additional commands, e-mail: [email protected]
