matrix is singular when doing LU Decomposition

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

matrix is singular when doing LU Decomposition

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

Reply | Threaded
Open this post in threaded view
|

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

Luc Maisonobe
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: [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

Luc Maisonobe
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: [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]