Return-Path: X-Original-To: apmail-commons-user-archive@www.apache.org Delivered-To: apmail-commons-user-archive@www.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 0C9F110875 for ; Mon, 9 Dec 2013 12:23:05 +0000 (UTC) Received: (qmail 18990 invoked by uid 500); 9 Dec 2013 12:23:02 -0000 Delivered-To: apmail-commons-user-archive@commons.apache.org Received: (qmail 18517 invoked by uid 500); 9 Dec 2013 12:22:56 -0000 Mailing-List: contact user-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: "Commons Users List" Delivered-To: mailing list user@commons.apache.org Received: (qmail 18503 invoked by uid 99); 9 Dec 2013 12:22:55 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 09 Dec 2013 12:22:55 +0000 X-ASF-Spam-Status: No, hits=-0.7 required=5.0 tests=RCVD_IN_DNSWL_LOW,SPF_PASS X-Spam-Check-By: apache.org Received-SPF: pass (athena.apache.org: domain of andrea.antonello@gmail.com designates 74.125.82.175 as permitted sender) Received: from [74.125.82.175] (HELO mail-we0-f175.google.com) (74.125.82.175) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 09 Dec 2013 12:22:51 +0000 Received: by mail-we0-f175.google.com with SMTP id t60so3296334wes.6 for ; Mon, 09 Dec 2013 04:22:30 -0800 (PST) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20120113; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; bh=rlwC45JTmZ++n5Dl1jb6I+JvqVFfc0mTvwb4fpXsaJU=; b=r9JlmEHfgpkvXYqSQLP1UUvTPi2cczleHLv7CDsgZL/CxCGv5aoGO0wfGE2CEtxb5/ dU9xOr1+Xbqc5c2fjN8aZ/blWHG5S0Nxbt/SfeVm0/ht1c4Lijl0KGPsKJeJmZcYRUyE NlVYzk/pOeTXWBxFmJLTAcYA9N9gE9Aj903fd1ygzXQDll60uF17G69MGLiZVQ/FL04q IA8xYg8vh300U2zqnL+B3YJ0dF9BNKV55GgGl/t1sflVDG97NcBFzP+4b6E41yqiS1x5 LuSo2JRkKH2ZgeN7Udq+SJxzMI22iQdMPSgHHoFaK6BssG11Yi/6iJbHjnI5LObalAlR Pnyg== X-Received: by 10.180.108.42 with SMTP id hh10mr13991306wib.15.1386591749998; Mon, 09 Dec 2013 04:22:29 -0800 (PST) MIME-Version: 1.0 Received: by 10.217.0.142 with HTTP; Mon, 9 Dec 2013 04:22:09 -0800 (PST) In-Reply-To: References: <52A19FB9.3040904@spaceroots.org> <52A1F8BE.4060509@spaceroots.org> From: andrea antonello Date: Mon, 9 Dec 2013 13:22:09 +0100 Message-ID: Subject: Re: [math] solving bivariate quadratic equations To: Commons Users List Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable X-Virus-Checked: Checked by ClamAV on apache.org 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 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 colum= n > should be the constant 1. We want to solve > > A u =3D z > > where u =3D [a, b, c, d, e, f]' and z =3D [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 =3D 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 =3D I and R is upper triangular. The normal equation can be re-written= as > > (Q R)' (Q R) u =3D (Q R)' z > R' R u =3D 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 th= at > 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: > > DecompositionSolversolver > =3D new > RRQRDecomposition > (A).getSolver > (); > RealVector u =3D > solver.solve > (z); > > > > > On Sat, Dec 7, 2013 at 7:44 AM, andrea antonello > 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 point= s >> > are considered to be more accurate than others. Typically, each point >> > should have a weight related to the variance of the error for the poin= t, >> > 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 sam= e >> > 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 weig= ht >> > 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 >> wrote: >> >>> Hi Andrea, >> >>> >> >>> Le 06/12/2013 10:02, andrea antonello a =E9crit : >> >>>> Hi Ted, >> >>>> thanks for the reply. >> >>>> >> >>>>> How would you like to handle the fact that you may have an infinit= e >> 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 (to= o >> >>>> few it seems). >> >>>> In my usecase I am trying to define a terrain model from sparse poi= nts >> >>>> (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 elevati= on >> >>>> (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 suc= h >> >>> 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 y= our >> >>> 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 simpl= e >> >>> (linear Z), both should give almost immediate convergence, pick up o= ne >> >>> 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 =3D ...; >> >>> final double[] yArray =3D ...; >> >>> final double[] zArray =3D ...; >> >>> >> >>> // start values for a, b, c, d, e, f, all set to 0.0 >> >>> double[] initialGuess =3D new double[6]; >> >>> >> >>> // model: this computes the theoretical z values, given a current >> >>> // estimate for the parameters a, b, c ... >> >>> MultivariateVectorFunction model =3D >> >>> new MultivariateVectorFunction() { >> >>> double[] value(double[] parameters) { >> >>> double[] computedZ =3D new double[zArray.length]; >> >>> for (int i =3D 0; i < computedZ.length; ++i) { >> >>> double x =3D xArray[i]; >> >>> double y =3D yArray[i]; >> >>> computedZ[i] =3D 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 =3D >> >>> new MultivariateMatrixFunction() { >> >>> double[][] value(double[] parameters) { >> >>> double[][] jacobian =3D >> >>> new double[zArray.length][parameters.length]; >> >>> for (int i =3D 0; i < computedZ.length; ++i) { >> >>> double x =3D xArray[i]; >> >>> double y =3D yArray[i]; >> >>> jacobian[i][0] =3D x * x; >> >>> jacobian[i][1] =3D y * y; >> >>> jacobian[i][2] =3D x * y; >> >>> jacobian[i][3] =3D x; >> >>> jacobian[i][4] =3D y; >> >>> jacobian[i][5] =3D 1; >> >>> } >> >>> return jacobian; >> >>> } >> >>> }; >> >>> >> >>> // convergence checker >> >>> ConvergenceChecker checker =3D >> >>> new SimpleVectorValueChecker(1.0e-6, 1.0e-10); >> >>> >> >>> // configure optimizer >> >>> LevenbergMarquardtOptimizer optimizer =3D >> >>> new LevenbergMarquardtOptimizer(). >> >>> withModelAndJacobian(model, jacobian). >> >>> withTarget(zArray). >> >>> withStartPoint(initialGuess). >> >>> withConvergenceChecker(checker); >> >>> >> >>> // perform fitting >> >>> PointVectorValuePair result =3D optimizer.optimize(); >> >>> >> >>> // extract parameters >> >>> double[] parameters =3D 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 skeleto= n. >> >>> >> >>> 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 =3D 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 >> > >> >> --------------------------------------------------------------------- >> 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