Andreas Andreakis a écrit :
hi,
sorry to send you an email regarding commons math directly. I tried to
subscribe to the mailing list, but for some reason, I was considered as
spam.
Thus I chose two members randomly.
Well, anyway :) I want to find a solution to the following problem:
(explained by example)
1 *c1 + 4 * c2 + 4* c3 + 4 * c4 = 2.5
1 *c1 + 4 * c2 + 1* c3 + 4 * c4 = 3.0
here the constraints:
1) c1+c2+c3+c4 = 1 (of course this could be expressed as an additional
equation: c1 + c2 + c3 + c4 = 1.0)
2) c1,c2,c3 and c4 are all positive
Important is that I need only one and feasible solution, thus the
computation should be as fast as possible. So there is no need to search
for every possible result.
And it needs to work with doubles (not only integers, because some
libraries work only with integers)
Actually the solution needs to scale for several dozens c´s ( there will
not be only 4 like c1,c2,c3,c4), but more like 50 c1, c2,...,c49, c50
And there will not be only two equations (like in the example) but
several hundred.
So the questions are:
1) can I find a solution using commons math. There is an optimization
package, but I could not figure out how to use it?
In this case, the estimation package is more suited to your needs. There
are much better methods for such problems, but they are not available in
commons-math.
The Levenberg-Marquardt solver implemented in thie estimation package
should support both under- and over-determined problems (i.e. either
more equations than variables or more variables than equations). The
case with more variable than equations is still considered experimental,
though. I did use it for such non homogeneous linear problems without
problems myself.
This solver does not completely fit your needs as is because it
currently do not handle constraints. You will have to set up an
intermediate unconstrained problem that can be solved and from which the
solution of your constrained problem can be derived. Here are a few
tricks to do that, take care, they may induce numerical problems:
For equality constraints, either add the constraint to the equation set
or remove one parameter (say c1) from the parameter set and compute it
as needed by solving the constraint each time from the values of the
other ones, i.e. compute c1 = 1 - f(c2, c3, ..., cn).
For simple bounds constraints on parameters (such as ck > ak), replace
the real parameter ck by an unconstrained parameter pk using
ck = ak + pk² (or ak - pk² if the constraint was ck < ak)
For interval constraints on parameters (such as lowk < ck < highk),
replace the real parameter ck by an unconstrained parameter pk using
ck = [ lowkk + highk ] / 2 + atan(pk) * [ highk - lowk ] / pk
Once you have changed your problem using these pk, it is an
unconstrained problem that can be solved by our implementation of
Levenberg-Marquardt. Once you have the pk values, you compute the ck values.
2) how will the performance be ?
Since the problem is linear, the performance should be fair.
3) could you give me a code snipped which is able to compute the given
example ? (I will generalize it then myself)
Try something along these lines. Take care, I wrote this exemple
directly in the mail, I did not attempt to even compile it ;-)
Hope this helps,
Luc
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.EstimationException;
import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
import org.apache.commons.math.estimation.SimpleEstimationProblem;
import org.apache.commons.math.estimation.WeightedMeasurement;
public class LinearConstrainedProblem extends SimpleEstimationProblem {
public static void main(String[] args) {
try {
// create the uninitialized problem
LinearConstrainedProblem pb =
new LinearConstrainedProblem(0.25, 0.25, 0.25, 0.25)
// 1 *c1 + 4 * c2 + 4* c3 + 4 * c4 = 2.5
pb.addEquation(1.0, 4.0, 4.0, 4.0, 2.5);
// 1 *c1 + 4 * c2 + 1* c3 + 4 * c4 = 3.0
pb.addEquation(1.0, 4.0, 1.0, 4.0, 3.0);
// solve the problem, using a Levenberg-Marquardt algorithm
// with default settings
LevenbergMarquardtEstimator estimator =
new LevenbergMarquardtEstimator();
estimator.estimate(pb);
System.out.println("RMS after solve = "
+ estimator.getRMS(pb));
// display the results
System.out.println("c1 = " + pb.getC1()
+ ", c2 = " + pb.getC2()
+ ", c3 = " + pb.getC3()
+ ", c4 = " + pb.getC4());
} catch (EstimationException ee) {
System.err.println(ee.getMessage());
}
}
/** Fitter for a linear constrained model.
* @param c1Guess start value for c1 (in fact, ignored)
* @param c2Guess start value for c2
* @param c3Guess start value for c3
* @param c4Guess start value for c4
*/
public LinearConstrainedProblem() {
// one parameter (c1) is computed from the other ones
// and the other ck are derived from pk using ck = pk*pk
p2 = new EstimatedParameter("p2", Math.sqrt(c2Guess));
p3 = new EstimatedParameter("p3", Math.sqrt(c3Guess));
p4 = new EstimatedParameter("p4", Math.sqrt(c4Guess));
// provide the parameters to the base class which
// implements the getAllParameters and
// getUnboundParameters methods
addParameter(c2);
addParameter(c3);
addParameter(c4);
}
/** Add an equation.
* @param a1 coeff for c1
* @param a2 coeff for c2
* @param a3 coeff for c3
* @param a4 coeff for c4
* @param rhs right hand side of the equation
*/
public void addEquation(double a1, double a2, double a3,
double a4, double rhs) {
addMeasurement(new Equation(a1, a2, a3, a4, rhs));
}
public double getC1() {
// this solves the equality constraint
return 1.0 - getC2() - getC3() - getC4();
}
public double getC2() {
double p = p2.getEstimate();
return p * p;
}
public double getC3() {
double p = p3.getEstimate();
return p * p;
}
public double getC4() {
double p = p4.getEstimate();
return p * p;
}
/** Internal measurements class, representing one equation. */
private class Equation extends WeightedMeasurement {
public Equation(double a1, double a2, double a3,
double a4, double rhs) {
super(1.0, rhs);
this.a1 = a1;
this.a2 = a2;
this.a3 = a3;
this.a4 = a4;
}
public double getTheoreticalValue() {
return a1 * getC1() + a2 * getC2()
+ a3 * getC3() + a4 * getC4();
}
public double getPartial(EstimatedParameter parameter) {
if (parameter == p2) {
return 2.0 * p2.getEstimate() * (a2 - a1);
} else if (parameter == p3) {
return 2.0 * p3.getEstimate() * (a3 - a1);
} else {
return 2.0 * p4.getEstimate() * (a4 - a1);
}
}
private final double a1;
private final double a2;
private final double a3;
private final double a4;
}
private EstimatedParameter p2;
private EstimatedParameter p3;
private EstimatedParameter p4;
}
kind regards,
Andreas
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
