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 re-written as
(Q R)' (Q R) u = (Q R)' z
R' R u = R' Q' z
This allows back-substitution 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/commons-math/apidocs/org/apache/commons/math3/linear/DecompositionSolver.html>solver
= new
RRQRDecomposition<http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/linear/RRQRDecomposition.html>
(A).getSolver<http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/linear/RRQRDecomposition.html#getSolver()>
();
RealVector u =
solver.solve<http://commons.apache.org/proper/commons-math/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 <[email protected]
> 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 <[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]
> >>
> >>
> >>
> >
> >
> > ---------------------------------------------------------------------
> > 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]
>
>
