# matrix is singular when doing LU Decomposition

 Classic List Threaded
3 messages
Reply | Threaded
Open this post in threaded view
|

## matrix is singular when doing LU Decomposition

 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?    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
Reply | Threaded
Open this post in threaded view
|

## [math] Re: matrix is singular when doing LU Decomposition

 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 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 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 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: [hidden email] > For additional commands, e-mail: [hidden email] > > --------------------------------------------------------------------- To unsubscribe, e-mail: [hidden email] For additional commands, e-mail: [hidden email]
Reply | Threaded
Open this post in threaded view
|

## Re: [math] Re: matrix is singular when doing LU Decomposition

 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> 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> 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> 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: [hidden email] >> For additional commands, e-mail: [hidden email] >> >> > > > --------------------------------------------------------------------- > To unsubscribe, e-mail: [hidden email] > For additional commands, e-mail: [hidden email] > > --------------------------------------------------------------------- To unsubscribe, e-mail: [hidden email] For additional commands, e-mail: [hidden email]