Le 11/10/2013 09:48, Luc Maisonobe a écrit :
> Hi Li Li,
>
> First of all, please note that this mailing list is shared among many
> Apache Commons components, so the name of the component (here [math])
> must be included in the subject line as I did when answering your
> question), so it helps the numerous subscriber to filter.
>
> Le 11/10/2013 08:31, Li Li a écrit :
>> hi all
>> I am using commons math to solve a linear equation. When I use
>> LUDecomposition, it throws a SingularMatrixException. If I switch it
>> to QRDecomposition, it's ok.
>> The problem is a simple polynomial regression. it is taken from
>> Bishop's book "pattern recognition and machine learning", chapter 1.
>> the data is generated by sin*2*pi*x + gaussian noise.
>> x=[0,1/9,2/9,...,1]
>> t=sin*2*pi*x + gaussian noise with starndard deviation of 0.3
>>
>> y(w,x)=w0 + w1 * x +w2 * x^2 + .. + wm * x^m
>> loss function: E(w)=1/2 {t1-y(w,x1)}^2 + ....
>> take differentiating with respect to wi
>> we get coefficients w = {wi} that minimize this error function are
>> given by the solution to the following
>> set of linear equations:
>> Aw=T
>> Aij=(x1)^(i+j) + (x2)^(i+j) + .. + (xn)^(i+j), A is a (M+1)*(M+1)
>> symmetric matrix
>> Ti=(x1)^i * t1 + (x2)^i * t2 + ... + (xn)^i * tn, T is a (M+1) vector
>> w=(w0,w1,...wM)' is a (M+1) vector
>>
>> by solve w, we can find best w to minimize loss function.
>>
>> for n=10 and M=0 to 9, it's correct.
>> But for n=10 and M=10, it throws an exception.
>>
>> what's wrong with it?
>
> It is a normal behaviour.
> You are building a Vandermonde matrix, which may have a large condition
> number and is notoriously difficult to solve with Gaussian elimination
> (which is what LUDecomposition does).
>
> You can check that in order to solve your case with this method, you
> have to reduce the singularity threshold from the default 1.0e-11 to
> something like 1.0e-2 by changing the line
the line above should refer to 1.0e-12, of course, not 1.0e-2
Luc
>
> DecompositionSolver solver =
> new LUDecomposition(coefficients).getSolver();
>
> into
>
> DecompositionSolver solver =
> new LUDecomposition(coefficients, 1.0e-12).getSolver();
>
>
> However, I would not suggest it. LU decomposition is simply not the way
> to handle such matrices. QR decomposition on the other hand is much more
> robust and should be privileged in many cases.
>
> Also note that Apache Commons Math provides a dedicated polynomial curve
> fitter (see the fitting package).
>
> best regards,
> Luc
>
>
>>
>> Following is my codes:
>>
>> PolynomialCurveFit .java
>>
>> public class PolynomialCurveFit {
>>
>> /**
>> * @param args
>> */
>> public static void main(String[] args) {
>> SyntheticRegressionDataGenerator srd=new SyntheticRegressionDataGenerator();
>> Point2D[] trainingData=srd.generateDataXWithFixedInterval(10, 0.3);
>> Point2D[] testData=srd.generateData(10, 0.3);
>> PolynomialCurveFit pcf=new PolynomialCurveFit();
>> System.out.println("\nLU Decomposition\n");
>> for(int order=0;order<=12;order++){
>> System.out.println("order: "+order);
>> try{
>> double[] weights=pcf.solveByCommonMathLUDecomposition(trainingData, order);
>> double trainErr=pcf.calcRootMeanSquareError(trainingData, weights, true);
>> double testErr=pcf.calcRootMeanSquareError(testData, weights, false);
>> System.out.println("training error: "+trainErr);
>> System.out.println("test error: "+testErr);
>> System.out.println();
>> }catch(Exception e){
>> e.printStackTrace();
>> }
>> }
>> System.out.println("\nQR Decomposition\n");
>> for(int order=0;order<=12;order++){
>> System.out.println("order: "+order);
>> try{
>> double[] weights=pcf.solveByCommonMathQRDecomposition(trainingData, order);
>> // System.out.print("Weights:");
>> // for(double weight:weights){
>> // System.out.print("\t"+weight);
>> // }
>> // System.out.println();
>> double trainErr=pcf.calcRootMeanSquareError(trainingData, weights, true);
>> double testErr=pcf.calcRootMeanSquareError(testData, weights, false);
>> System.out.println("training error: "+trainErr);
>> System.out.println("test error: "+testErr);
>> System.out.println();
>> }catch(Exception e){
>> System.out.println("Fail");
>> }
>> }
>> }
>> public double calcRootMeanSquareError(Point2D[] array,double[]
>> weights,boolean isTraining){
>> double error=0;
>> for(Point2D point:array){
>> double x=point.getX();
>> double y=this.predictByPolynomial(x, weights);
>> double realY;
>> if(!isTraining){
>> realY=Math.sin(2*Math.PI*x);
>> }else{
>> realY=point.getY();
>> }
>> error+=(y-realY)*(y-realY);
>> }
>> error=Math.sqrt(error/array.length);
>> return error;
>> }
>> private double predictByPolynomial(double x,double[] weights){
>> double result=0;
>> for(int i=0;i<weights.length;i++){
>> result+=weights[i]*Math.pow(x, i);
>> }
>> return result;
>> }
>> public double[] solveByCommonMathLUDecomposition(Point2D[] array,int order){
>> Array2DRowRealMatrix coefficients =new Array2DRowRealMatrix(order+1,
>> order+1);
>> ArrayRealVector constants=new ArrayRealVector(order+1);
>> for(int i=0;i<=order;i++){
>> double Ti=0;
>> for(Point2D point:array){
>> Ti+=point.getY()*Math.pow(point.getX(), i);
>> }
>> constants.setEntry(i, Ti);
>> }
>> for(int i=0;i<=order;i++){
>> for(int j=0;j<=order;j++){
>> double Aij=0;
>> for(Point2D point:array){
>> Aij+=Math.pow(point.getX(), i+j);
>> }
>> coefficients.setEntry(i, j, Aij);
>> }
>> }
>>
>> DecompositionSolver solver = new LUDecomposition(coefficients).getSolver();
>> RealVector solution = solver.solve(constants);
>> return solution.toArray();
>> }
>> public double[] solveByCommonMathQRDecomposition(Point2D[] array,int order){
>> Array2DRowRealMatrix coefficients =new Array2DRowRealMatrix(order+1,
>> order+1);
>> ArrayRealVector constants=new ArrayRealVector(order+1);
>> for(int i=0;i<=order;i++){
>> double Ti=0;
>> for(Point2D point:array){
>> Ti+=point.getY()*Math.pow(point.getX(), i);
>> }
>> constants.setEntry(i, Ti);
>> }
>> for(int i=0;i<=order;i++){
>> for(int j=0;j<=order;j++){
>> double Aij=0;
>> for(Point2D point:array){
>> Aij+=Math.pow(point.getX(), i+j);
>> }
>> coefficients.setEntry(i, j, Aij);
>> }
>> }
>> DecompositionSolver solver = new QRDecomposition(coefficients).getSolver();
>> RealVector solution = solver.solve(constants);
>> return solution.toArray();
>> }
>>
>> }
>>
>> SyntheticRegressionDataGenerator.java
>>
>> public class SyntheticRegressionDataGenerator {
>> public Point2D[] generateDataXWithFixedInterval(int N,double
>> noiseStandardDeviation){
>> if(N<1){
>> throw new IllegalArgumentException("N must great than 0. N="+N);
>> }
>>
>> RandomDataGenerator rdg=new RandomDataGenerator(new MersenneTwister());
>> Point2D[] array=new Point2D[N];
>> double interval=1.0/(N-1);
>> double x=0;
>> double y=0;
>> for(int i=0;i<N;i++){
>> y=Math.sin(2*Math.PI*x)+rdg.nextGaussian(0, noiseStandardDeviation);
>> array[i]=new Point2D(x, y);
>> x+=interval;
>> }
>> return array;
>> }
>> public Point2D[] generateData(int N,double noiseStandardDeviation){
>> if(N<1){
>> throw new IllegalArgumentException("N must great than 0. N="+N);
>> }
>> RandomDataGenerator rdg=new RandomDataGenerator(new MersenneTwister());
>> Point2D[] array=new Point2D[N];
>>
>> double x=0;
>> double y=0;
>> for(int i=0;i<N;i++){
>> x=rdg.nextUniform(0, 1.0, true);
>> y=Math.sin(2*Math.PI*x)+rdg.nextGaussian(0, noiseStandardDeviation);
>> array[i]=new Point2D(x, y);
>>
>> }
>> return array;
>> }
>>
>> }
>>
>> ---------------------------------------------------------------------
>> 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]
