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
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]
